Modelling of ocean waves with the Alber equation: application to non-parametric spectra and generalization to crossing seas

The Alber equation is a phase-averaged second-moment model for the statistics of a sea state, which recently has been attracting renewed attention. We extend it in two ways: firstly, we derive a generalized Alber system starting from a system of nonlinear Schr\"odinger equations, which contains the classical Alber equation as a special case but can also describe crossing seas, i.e. two wavesystems with different wavenumbers crossing. (These can be two completely independent wavenumbers, i.e. in general different directions and different moduli.) We also derive the associated 2-dimensional scalar instability condition. This is the first time that a modulation instability condition applicable to crossing seas has been systematically derived for general spectra. Secondly, we use the classical Alber equation and its associated instability condition to quantify how close a given non-parametric spectrum is to being modulationally unstable. We apply this to a dataset of 100 non-parametric spectra provided by the Norwegian Meteorological Institute, and find the vast majority of realistic spectra turn out to be stable, but three extreme sea states are found to be unstable (out of 20 sea states chosen for their severity). Moreover, we introduce a novel"proximity to instability"(PTI) metric, inspired by the stability analysis. This is seen to correlate strongly with the steepness and Benjamin-Feir Index (BFI) for the sea states in our dataset (>85% Spearman rank correlation). Furthermore, upon comparing with phase-resolved broadband Monte Carlo simulations, the kurtosis and probability of rogue waves for each sea state are also seen to correlate well with its PTI (>85% Spearman rank correlation).


Introduction
Ocean waves are a very active field of mathematical modelling and analysis. The first principles hydrodynamic models for gravity waves are by now well understood [1,31,32]. An array of approximate models is well established and widely used, including the nonlinear Schrödinger equation (NLS) and its variants [11,54], the Zakharov equation [57], the coupled-mode systems [9,10,45], the High Order Spectral Method (HOSM) [20,56,21] and others. (In shallow water an even larger collection of models is being used, but here we focus on deep water.) One reason for the wide use of approximate models in oceanography is the need to study large wavefields, with hundreds or thousands of individual wavelengths. Thus there is a trade-off between hydrodynamic fidelity and the ability to scale up models, as well as the accessibility of powerful qualitative insights.
In fact, actual ocean waves include multi-physics and non-differentiable phenomena, such as wind forcing, wave breaking etc. These are not accounted for in the classical "exact" hydrodynamic model anyway, and are still being understood on different levels [15,14,6,16,27].
In this paper we will use stochastic modelling of ocean waves and, furthermore, explore what phaseaveraged stochastic models may reveal about rogue waves in particular. This will lead us to novel mathematical results; it is also worth mentioning that this kind of overall approach has been identified as a priority within the broader marine research & industry community [24].
The most well-known stochastic models for ocean waves include the CSY equation [18,3,4], Hasselmann's equation [26] and the Alber equation [2] that we will focus on here. For a recent review of various stochastic models one can see [51]. Broadly speaking, they are moment equations, starting from phaseresolved equations for the sea surface (such as Zakharov's equation or the NLS) as an approximation for deterministic wave dynamics. One then takes stochastic moments of the deterministic equations; due to the nonlinearity of these equations, an infinite hierarchy of moments is produced. A gaussian second-order moment closure is then used, to produce a closed equation for the second stochastic moment. The resulting equation is phase-averaged, meaning that in no longer resolves individual wave peaks and troughs, but instead the evolution and propagation of the statistics of the wavefield.
A key approximation step is the gaussian moment closure. This is of course not exact, but in many cases the free surface is indeed close to being gaussian [29,39], making the gaussian closure plausible. It is important to keep in mind that fidelity to the deterministic, phase-resolved model is not the only consideration: for example, the exact infinite hierarchy of moments, is not just more complicated as a mathematical model: it is impossible to initialize meaningfully. The vast majority of synoptic data collected from the ocean is, or can be converted to, some kind of second moment [39]. There is some data involving moderately higher moments, but very little data or know-how exists for moments higher than 4 th order.
The question of using realistic data for initialisation is extremely important; all the more so in the study of extreme sea states and rogue waves. For example, every deterministic model can produce "rogue waves" on demand, by carefully preparing particular initial conditions; however, the real-life question is how often would these "initial conditions leading to rogue waves" realistically appear? A stochastic approach can directly attack this question e.g. by a phase-resolved Monte Carlo approach [52]; this has a number of advantages, but is clearly expensive to be applied indiscriminantly. Another possibility would be using phase averaged stochastic models, like the moment equations discussed above, to directly investigate whether a sea state is likely to support the rapid concentrations of energy. That way sea states of interest could be selected and computational resources focused on them. In fact, it appears that the sea states highlighted as more unstable by the Alber equation do turn out to exhibit a higher probability of extreme events in a phase-resolved Monte Carlo simulation (details in Section 3.2.3 and Figure 5).

Ocean wave modelling with the Alber equation
The Alber equation, its derivation and interpretation have been widely studied and explained. We refer the interested reader to [2,49,22,7] and the references therein for more details. Here we will briefly present its derivation and main features in order to make the paper self-contained.

Derivation
The cubic focusing nonlinear Schrödinger equation (NLS) is an approximate model for the envelope of a narrow-banded wavetrain with carrier wavenumber k 0 along its direction of propagation (unidirectional propagation), cf. e.g. [34]. Thus, in deep water the sea surface elevation ηpx, tq is related to the complex-valued envelope upx, tq through The coefficients of the equation also depend on k 0 " |k 0 |, This is an asymptotic model, where the order parameter is the steepness of the waves, and moreover assumes narrow-bandedness of the wavefield around the carrier wavenumber. The NLS and its variants are widely used as relevant to a wide array of realistic scenarios [52]. 1 A different kind of limitation of the NLS is that it is a deterministic phase-resolved model; i.e. any prediction with it will not be better or more accurate than the initial condition used. But realistic wave-systems are not widely available as phaseresolved initial conditions, as opposed to power spectra. In that context, Alber [2] proposed generating a second-order moment equation from the NLS, considered with stochastic initial data. Denoting by Rpx, y, tq " Erupx, tqupy, tqs the autocorrelation of the envelope u, one obtains iB t Rpx, y, tq`p p∆ x´∆y q Rpx, y, tq`qE " upx, tqupy, tqrupx, tqupx, tq´upy, tqupy, tqs By using the gaussian closure Er|upα, tq| 2 upα, tqupβ, tqs " 2Rpα, α, tqRpα, β, tq, the autocorrelation can now be seen to satisfy iB t Rpx, y, tq`p p∆ x´∆y q Rpx, y, tq`2qRpx, y, tqrRpx, x, tq´Rpy, y, tqs " 0.
A key restriction the Alber equation inherits from its starting point, the NLS equation, is narrowbandedness. This is arguably not too restrictive for unidirectional sea states [52]; however it completely fails for crossing seas, i.e. sea states where several different directions carry substantial wave energy. To address this limitation, we will derive a generalized Alber equation valid in a crossing seas scenario. This is made possible by starting from a system of NLS equations describing the crossing wavetrains which was derived in [25]. The resulting generalized Alber equation is reported in Section 3, and derived in detail in Section 4. Crossing seas are receiving increased attention as a possible incubator of rogue waves [42,23], so it is significant to have a moment equation applicable to such scenarios.

The stability-of-homogeneity question
It is empirically known that sea-states are typically homogeneous and stationary, at least for appropriate lengthscales and timescales [29,39,30,53]. This feature is reflected in the Alber equation; for example observe that Rpx, y, tq " Γpx´yq is a solution of equation (4) for any smooth function Γ. The next step is to investigate the stability of such homogeneous and stationary solutions: let us consider a weakly inhomogeneous initial sea-state, i.e. assume that the autocorrelation is initially of the form Rpx, y, 0q " Erupx, 0qupy, 0qs " Γpx´yq` ρpx, y, 0q, ! 1 for some nice functions Γ, ρ (e.g. Schwartz class test functions). By inserting eq. (5) in eq. (4) one can reformulate the problem in terms of the inhomogeneity ρ, iB t ρpx, y, tq`p p∆ x´∆y q ρpx, y, tq`2qrΓpx´yq` ρpx, y, tqsrρpx, x, tq´ρpy, y, tqs " 0.
So stability for homogeneous state-states is controlled by the boundedness (or lack thereof) of the inhomogeneity ρ in equation (6). Can ρ grow in time to the extent that ρpx, y, tq is no longer small? Or is there a guarantee that ρ stays bounded? In the latter case of stability, the autocorrelation will simply stay close to Γpx´yq for all times. In the unstable case however, even if initially very close to homogeneous, the autocorrelation could develop significant inhomogeneities.
In [2] a sufficient condition for linear instability was derived in terms of the spectrum Spkq, namely the Fourier transform of the autocorrelation function Γ, Spkq " F yÑk rΓpyqs.
Indeed, it was shown that the homogeneous sea-state with autocorrelation Γpx´yq is unstable if, for some X P R, there exists ΩpXq P C so that This was called an "eigenvalue relation" in [2]; we will call it an "instability condition". In [7] it was further shown that if the instability condition does not hold then linear stability follows. The instability condition (7) itself can be refined in two ways: one concerns a technical issue related to X " 0. However, the more serious one is that in (7) we are asked to guarantee the existence or not of solutions for a non-linear system of two equations (the real and imaginary part of eq. (7)) in three real unknowns (X, Re Ω, Im Ω). This is not straightforward in general, and historically it has been a challenge in using the Alber equation more widely [22]. In [7] this condition is reformulated so that a more constructive way to check it can be found: by dividing both sides of the fraction by X and setting This becomes very simple if we recognize it as the Hilbert transform of the divided difference of a rescaled spectrum. So, denoting P pkq :" k 3 0 Spkk 0 q, D X P pkq " # P pk`X 2 q´P pk´X 2 q X , X ‰ 0 P 1 pkq, and dropping the primes, the condition for instability finally becomes DX P R : DΩ " ΩpXq P C : HrD X P spΩq " 1 4π .
The benefit is that now the argument principle 2 can be used to reformulate this to a constructive condition that we can directly check [46,8,7]. To that end we will also need to introduce the signal transform, Sruspxq :" Hruspxq´iupxq; we are now ready to state the equivalent instability condition: . Consider a homogeneous sea state with autocorrelation function Γpx´yq and carrier wavenumber k 0 , and denote Spkq " F yÑk rΓpyqs, P pkq :" k 3 0 Spkk 0 q. Then, the sea state is Penrose-Alber unstable if Any X for which 1{4π PΓ X is called an unstable wavenumber.
This notion of instability is equivalent to condition (9) [7]. The closed curves Γ X are vizualized for a concrete example in the top of Figure 1.
We call this the Penrose-Alber instability condition after Alber's "eigenvalue relation" (7) [2] and Penrose's introduction of the argument principle in an analogous problem in plasma [46]. The point is that this formulation of the condition boils down to drawing closed curves on the complex plane, and looking whether the point 1{4π is inside or outside of them, cf. Figure 2. It was used in [7] to study parametric JONSWAP spectra; here we will apply it to non-parametric spectra as well.
This stability-of-homogeneity question that was first studied in [2] is in fact a variant of a much better known stability question. The modulation instability (MI) is well understood and widely studied in a phase-resolved setting, for a plane wave background [58]. It can easily be seen that a sequence of spectra approaching a plane wave Spkq " aδpk´k 0 q would become unstable in the Penrose-Alber sense. In fact, the Penrose-Alber instability is just the modulation instability for more general backgrounds than plane waves. A modulationally unstable sea state supports the rapid concentrations of energy -a possible mechanism for the formation of rogue waves [44,41,43,8].
Unlike the classical MI, where every plane-wave solution is always unstable, a spectrum can be stable or unstable in the Penrose-Alber sense. In the case of stability, the homogeneity of the sea-state is robust, and small perturbations will merely disperse. In these cases, despite having a focusing nonlinearity and infinite energy present, the dynamics are going to be dominated by the linear dispersion for all times. This is a familiar leading order approximation in ocean waves, but it doesn't necessarily have a name to itself in an oceanographic context. As the rigorous stability analysis of [7] highlighted, mathematically this stable regime looks exactly like what is called Landau damping for Vlasov equations (indeed this parallel was also drawn before [40]). Figure 1: Curves Γ X for a nonparametric spectrum on the complex plane. The smaller closed curves correspond to larger |X|; here X " p1`120¨nq5¨10´4, n " 0, . . . , 7. Taking |X| ă 5¨10´4 does not change the outermost curve noticeably, as D X P has effectively converged to P 1 . The real number 1{4π « 0.08 is not circumscribed by any of the curves, i.e. the spectrum does not exhibit modulation instability (cf. Definition 2.1). The spectrum used is a unimodal spectrum with a reference steepness of around 7%, taken out of the dataset of 100 nonparametric spectra.

Quantifying stability
As was seen just above, a key feature of the Alber equation is a classification of a given spectrum as either stable or unstable. A careful look at the asymptotics can offer more nuance. For example, consider two sea states: one with a barely stable spectrum Spkq, and the other with a slightly perturbed version, e.g. p1`εqSpkq, so that it becomes barely unstable. These sea states will behave similarly on physically realistic timescales -as one would intuitively expect. In particular, there is not a violent bifurcation from Landau damping to modulation instability, but rather a gradual transition [7,8]. On the other hand, as we will see in some detail here, a spectrum exhibiting Landau damping can be "more stable" than another spectrum which also exhibits Landau damping. In short, the effective stability of a spectrum is better thought of as belonging to a continuous range of values rather than a binary "stable / unstable" classification. One of the results of this paper is a non-dimensional index quantifying this "effective stability", cf. Figures 2,4. This is important because, even in the presence of Landau damping (i.e. for "stable" spectra) only small enough inhomogeneities are guaranteed to disperse if the nonlinearity is taken into account. Large inhomogeneities are not well understood, and could very well behave differently. This is a harder mathematical problem, because now by definition asymptotics are not sufficient. To better understand this question there is a numerical requirement (efficient and reliable solvers for the Alber equation both in the stable and unstable regime) as well as a data requirement (quantify how large inhomogeneities are in the ocean; this may be possible e.g. with X-band radar imaging [13]). It is quite possible that, if a spectrum is close enough to MI the likelihood of nonlinear events under realistic perturbations may substantially Figure 2: We can quantify how close each spectrum is to being modulationally unstable by measuring the distance by between Γ X and 1{4π (plotted here as a red star). In practice, it suffices to do this for X « 5¨10´4, as much smaller values of X yield similar results, and large values of X lead to smaller and smaller curves, cf. Figure 1. So now we can say that spectrum no. 11300 has come something like 42% of how far out it would need to be in order to be modulationally unstable.
increase, even if technically the spectrum still exhibits Landau damping.

Nonparametric spectra
The fundamental scaling of the problem shows that the vast majority of plausible sea states would be stable, in accordance with the well-known fact that linearized dynamics very often do a good job in describing phase-averaged energy propagation. On the other hand, it seems that instability is within reach, i.e. the scaling of the problem does not make instability so far removed as to be considered impossible. This much was established in [7,22,49], working in the context of fitted parametric JONSWAP spectra. While fitted spectra are widely used, realistic spectra from the field come in many different shapes and forms. In this paper we work with a set of nonparametric unidirectional spectra. We investigate whether they are stable or unstable in the Penrose-Alber sense, and proceed to examine how well the "proximity to instability" (defined more precisely in equation (19)) correlates with the probability of rogue waves appearing in the given sea state.
The dataset, methodology and results of our investigation are described in detail in Section 3.2.

Emergence of coherent structures
Another area with many open questions is what happens when instability arises. For modulationally unstable spectra with a small inhomogeneity, formal asymptotics can be used to describe the early evolution of the instability. A particular coherent structure emerges, determined by the unstable wavenumbers for the particular spectrum and their rate of growth. In that sense the coherent structure, at least in its early stages, is determined by the spectrum (and not by the inhomogeneity) and the Penrose-Alber instability analysis suffices to predict it [8]. Rather surprisingly, the same kind of universal coherent structure was reported in a fully numerical study by van den Eijden et al. [19], for the fully nonlinear stage of the instability, cf. Figure 3. This is a direction with many open questions, where more work is needed. Figure 3: Predicted profile of emergent localized extreme events for NLS, computed according to the methodology of [8]. A virtually identical universal profile of fully developed rogue waves is reported in Figure 2 of [19] for the largest extreme waves. This brings to mind the "three sisters" discussion [12,33,37,38,17,5], or the Greek τ ρικυµια.
3 Main results

The Alber equation for crossing seas
The original derivation of the Alber equation [2] allows for an oblique but small inhomogeneity on a unidirectional sea state. This does include some two dimensional aspects, however it does not allow for two large, different wavesystems, with different main directions of propagation, crossing. Such a situation is called "crossing seas", and is recently attracting a lot of interest. In particular, it is thought that modulation instability and rogue waves may be more prominent in crossing seas, but this is still very far from fully understood [23,50,25,41,42]. One way to study crossing seas is by deriving a coupled system of equations, each governing the evolution of one (quasi-uni-directional) wavefield. Let us denote the direction of propagation for the wavetrain A by k A " pk A 1 , k A 2 q and for the wavetrain B by k B " pk B 1 , k B 2 q. The corresponding frequencies are for ocean waves (i.e. in the limit of infinite depth). In [25] the following system of NLS equations for the envelopes of two wavetrains, v A , v B , in two spatial dimensions is derived: All the coefficients are completely determined in terms of k A , k B [25].
Theorem 3.1. Consider the two crossing wavefileds of equations (13), (14), and assume moreover that their autocorrelations at t " 0 can be written as for some " op1q. Then the system (13), (14) exhibits modulation instability if where h A , h B , r n 0 A , r n 0 B are defined in terms of the data of the problem in eq. (47) 3 . In that case, inhomogeneities are expected to grow in time.
On the other hand, if and sup Re ωą0 hold, then formally the problem exhibits linear Landau damping, i.e. inhomogeneities are expected to disperse and thus deviation from homogeneity to not grow noticeably.
The proof is found in Section 4.  (17) can also be resolved with the argument principle since, for each P P R 2 , it asks whether a holomorphic function attains the value 0 on the right half-plane, inf Re ωą0, PPR 2 |F P pωq| ą κ. Condition (18) apparently boils down to general regularity conditions for the Working these ideas out in full detail will require the extension of some technical results to a non-standard "two-dimensional Hilbert transform" which arises here.
3.2 Stability of unidirectional non-parametric spectra and Proximity to Instability (PTI)

The data
For the present analysis we have used 100 realistic nonparametric spectra taken from the Norwegian Meteorological Institute's operational spectral wave model [48], which is a third-generation wave model 3 Roughly speaking, h A and h B are transfer functions generated by the homogeneous backgrounds Γ A , Γ B , and r n 0 A , r n 0 B express the free space (i.e. ξ1 " ξ2 " ζ1 " ζ2 " 0) evolution of the initial inhomogeneities.
based on WAM. The model provides wave spectra every hour for many locations in the North Atlantic. For this study 100 spectra from one specific location south-west of the Norwegian coast (560 36 1 N, 30 12 1 E) were used. The selected spectra consisted of 80 spectra randomly selected from the full database (26 255 spectra covering the period from October 2016 to September 2019), as well as the 10 spectra having the largest mean wave steepness " H s k 0 {2 and the 10 spectra with the largest BFI values. Recall that the Benjamin-Feir Index (BFI) was defined from the frequency spectrum Epωq as where δ ω is a measure for the frequency-bandwidth here defined in terms of Godas' peakedness factor Q p , as suggested e.g. in [28]: and where m 0 " ş Epωqdω is the total energy of the spectrum. Note that for the following analysis the original frequency spectra Epωq were converted into wavenumber spectra Spkq using the linear dispersion relation ω " ? gk.

The algorithm for checking the instability condition
Here we will implement the stability criterion of Definition 2.1 to a number of non-parametric ocean spectra. First of all we will rescale them Spkq Þ Ñ P pkq :" k 3 0 Spk¨k 0 q. Actually the selection of k 0 is a nontrivial question. For example, one could plausibly use the mode, the mean or the median wavenumber. For narrow, unimodal spectra the choice would make very little difference, but for more irregular shapes, including e.g. bimodal spectra, the resulting differences might be noticeable. In this case, we use the k 0 provided in the dataset by the Norwegian Meteorological Institute, and which has been used for the computation of the BFI, steepness contained in the dataset etc.
Once k 0 is determined and the rescaling is complete, we will need to interpolate the discrete data to a finer grid. (The original discrete spectra come sampled in 36 non-uniformly spaced wavenumbers.) This will be crucial in using effective quadrature methods. Splines are well suited to this task; the additional requirement is to minimize overshooting at maxima, as this could make a big difference with regard to our investigation. To that end we use pchip, a piecewise polynomial interpolation routine in MATLAB that minimizes overshoot.
The non-trivial part of the computation is the Hilbert transform, which is a singular integral. To that end we will use the Sokhotski-Plemelj formula, @u P CpRq X L 1 pRq lim ηÑ0`1 π ż s upsq t´s´p ?´1 qη ds " Hrusptq´iuptq " Srusptq and truncate the limit by taking an appropriate η " compl tol ! 1. After extensive testing it is found that the result doesn't change noticeably once η « 10´4 That way we avoid the singular integral. A detailed pseudocode for how the curve Γ X is generated is presented in Algorithm 1. It is a moderately heavy computation if a good approximation for all of Γ X is required, as a few million points are typically required in order to achieve stringent error tolerances (" 10´2 relative error tolerance or " 10´6 absolute error tolerance). However, one can check stability more quickly, by checking only the points t˚where D X P pt˚q " 0; these are the points where the curve Γ X crosses the real axis. If it never crosses the real axis to the right of 1{4π, then topologically 1{4π cannot possibly be in the interior of the curve. Some more details involve the selection of X; it is found that X ! 1 will produce the curves that have the most chance of coming closer to 1{4π, as when X increases the curves Γ X shrink to zero, see Figure  1. Numerical testing shows that X « 10´4 gives a good picture of what happens as X Ñ 0. So, if for X « 10´4 Γ X is not winding around 1{4π, nor coming very close to it, we can accept the spectrum as stable.
Algortithm 1 Pseudo-code for the computation of Γ X Input S j , k j (Sampled values of wavenumber-resolved spectrum, k j P r0.00479, 3.78s.) Rescale pS j , k j q Þ Ñ pk 3 0 S j , k j {k 0 q " pP j , ξ j q (k 0 is simply taken to be the peak wavenumber.) Interpolate P pξq Set compl tol" 10´4, rel tol" 10´2, abs tol" 10´6 For t i :" t min to t max step δt While rel err ą rel tol AND abs err ą abs tol Integrate I " 1 π ż s P ps`X 2 q´P ps´X 2 q X t i´s´p ?´1 q compl tol ds using composite Simpson on two quadrature grids: a finer one and a coarser one, generating two approximations, I fine and I coarse.
(Fine grid has 3ˆthe number of points compared to coarse grid.) Set rel err " |I fine´I coarse| I fine , abs err " |I fine´I coarse|.

End While
Set Γpt i q " I fine

End For
Plot the line´Re Γ X pt i q, Im Γ X pt i q¯, i " 1, 2, ... and the point p 1 4π , 0q Check whether p 1 4π , 0q is inside Γ X (This can be done with the MATLAB inpolygon function.)

Summary of the results
Most of the non-parametric spectra examined were found to be stable, i.e. exhibit Landau damping. However, three spectra were found to be modulationally unstable, and a handful more were extremely close to being unstable. We also define as proximity to instability (PTI) the quantity PTI " 1´d pΓ, 1 4π q 1 4π ; (19) this is 1 for any modulationally unstable spectrum, and 0 for the zero spectrum. It provides a nondimensional way to quantify how close a spectrum comes to being modulationally unstable in the sense of Definition 2.1. We will compare this with quantities of interest for the same spectrum (cf. Figure 4), along with a Monte Carlo estimation of the likelihood of rogue waves and nonlinear events (cf. Figure 5) Figure 4: The new metric of PTI offers a quantitative and nondimensional way to assess how far a spectrum is from being modulationally unstable. The Benjamin-Feir Index (BFI), essentially a rescaled version of steepness, was introduced with a similar task in mind; so it is only natural to ask how well the two are correlated. In the left graph above we see a log-log scatterplot of the BFI v. PTI. In the right graph we see a similar scatterplot of a representative wave steepness, " H s k 0 {2, v. PTI.
Indeed, for each of the 100 spectra selected for this study, we have run numerical simulations with the Higher Order Spectral Method (HOSM) [20,56] in a Monte-Carlo approach where each spectrum were simulated 100 times with different initial random phases and random amplitudes each run. That is, for a given spectrum Spkq the initial surface elevation is in the form where ∆k j " k j`1´kj is the grid spacing between the discrete wavenumbers and where Z j are independent complex standard normal variables. That is, the real and imaginary parts of Z j are independent normally distributed random variables with zero mean and variance 1{2, meaning that |Z j | are Rayleigh distributed with parameter σ " 1{ ? 2 so that Er|Z j | 2 s " 1 and Er|A j | 2 s " 2Spk j q∆k j , and the phases ArgpZ j q are uniformly distributed on r0, 2πq.
HOSM applies a regular discretization of the wavenumbers so that ∆k " 2π{x max " 2π{n∆x, leading to a periodic domain of length x max in space. In the present simulations we have used n " 1024, representing wavenumbers up to k max " 8k 0 .. This means that ∆k " 2k max {n, corresponding to ∆x " λ 0 {16 and x max " 64λ p , where λ 0 " 2π{k 0 is the reference wavelength. Each simulation was run for 30 minutes, from which time series of surface elevation were extracted from four locations distributed over the simulation domain. Thus, 200 hours of surface elevation time series were obtained for each sea state. Since here we are interested in relations between the instability (or proximity of such) obtained from the stochastic approach (i.e. Alber equation) and the occurrence of rogue and extreme events in random realizations of the sea states (i.e. phase-resolved Monte-Carlo simulations), we will consider the following parameters related to the occurrence of extreme wave events: the sea surface kurtosis and the probability of extreme wave crests.
The sea surface kurtosis is a measure for how much the tail of the distribution deviates from Gaussian statistics. For a Gaussian distribution the kurtosis is equal to zero, while a positive value for kurtosis indicates more large events than in a Gaussian population. Hence, the kurtosis is often used as an indicator for the probability of extreme and rogue waves.
Secondly we consider the probability that a wave crest, defined as the maximum between each zerocrossing of the surface elevation, exceeds the significant wave height H s ; P pC ą H s q 4 . This is a common definition of a rogue wave, which under linear & gaussian assumptions (Rayleigh distributed crests) has the probability P pC ą H s q " exp p´8q. The probability of crest exceedance is estimated as the relative frequency of individual crests in the HOSM simulations that exceed the threshold H s . For brevity we will refer to this probability simply as "rogue wave probability".

Proof of Theorem 3.1
First of all observe that, by virtue of a simple gauge transform, we can eliminate the terms with the first derivatives 5 ; without loss of generality we will work with the simplified system Special symmetric cases of this system which simplify the coefficients have been studied in [41,50], but we can in fact derive an Alber system for the general case.
To proceed we will use some shorthand notations, namely x " px 1 , x 2 q, y " py 1 , y 2 q, L j Now equation (22) can be compactified as and similarly for eq. (23) with L 2 x . At this point let us introduce the autocorrelation functions for A, B; R A px, y, tq " ErApx, tqApy, tqs, R B px, y, tq " ErBpx, tqBpy, tqs.
The last two lines consist of fourth order stochastic moments, and it is these terms that will have to be approximated by some closure scheme. For the fourth order moments involving A only, we will use the same idea as in the standard Alber equation; namely we will use the fact that for Gaussian processes the relationship Er|Apx, tq| 2 Apx, tqApy, tqs " 2R A px, x, tqR A px, y, tq 5 By plugging (21) in eq. (13) we see that eliminating the linear terms containing A and ∇A leads to The first two lines is a 2x2 system, and for the third, nonlinear equation we only need to substitute the obtained κ A , λ A values from above. The same steps lead to a symmetric expression for B.
holds (cf. Theorem A.1). This now provides the Gaussian closure for all terms of the same form, since Er|Apy, tq| 2 Apx, tqApy, tqs " 2R A py, y, tqR A px, y, tq (28) follows by the same argument. The new kinds of terms that come into play for the first time are the joint A´B moments. In physical terms, A and B represent two different wavetrains meeting in the ocean, each having been generated and propagated independently from each other. It is thus a reasonable assumption that they are stochastically independent. In that case, the joint moments simplify to Er|Bpx, tq| 2 Apx, tqApy, tqs " R B px, x, tqR A px, y, tq, Er|Bpy, tq| 2 Apy, tqApx, tqs " R B py, y, tqR A px, y, tq.
So, using these closures we finally end up with the following equation (whenever the independent variables are not shown explicitly, it is understood that R A :" R A px, y, tq): By the symmetry of equations (22), (23) one readily sees that, using the same kinds of closures, Now, we invoke the assumption of eq. (15), i.e.
Now taking the Laplace transform in time r f pP, Q, ωq " L tÑω rf pP, Q, tqs, r gpP, Q, ωq " L tÑω rgpP, Q, tqs, r n A pP, ωq " L tÑω rq n A pP, tqs, r n B pP, ωq " L tÑω rq n B pP, tqs we obtain the system iω r f´ P, Q 1 r f`2`x Γ A pQ´P 2 q´x Γ A pQ`P 2 q˘rξ 1 r n A pP, ωq`ζ 1 r n B pP, ωqs " f pP, Q, 0q, By dividing with the free-space part we obtain r f`2 x Γ A pQ´P 2 q´x Γ A pQ`P 2 q iω´ P,Q 1 rξ 1 r n A`ζ1 r n B s " f pP,Q,0q iω´ P,Q 1 , r g`2 y Γ B pQ´P 2 q´y Γ B pQ`P 2 q iω´4π 2 P,Q 2 rξ 2 r n B`ζ2 r n A s " gpP,Q,0q iω´ P,Q 2 .
The two key observations here are the following: • the rhs of both equations correspond to the linear free-space solutions (and thus can be treated as known and well-behaved functions), moreover • we can now integrate both equations in the Q variables and obtain a closed system for the position densities r n A , r n B ; achieving this is in fact the motivation for all the transforms and changes of variables.
We can summarize the result by denoting so that system (46) becomes , . - leading to r n A " as long as the determinant (itself a function of P, ω) is nonzero, So the problem is altogether linearly stable if eqs. (17), (18) hold. If both conditions (17), (18) hold, then q n A pP, tq " L´1rr n A pP, ωqs, q n B pP, tq " L´1rr n B pP, ωqs inherit an L 2´t ype decay in time, and one expects that Landau damping estimates similar to what was done in the scalar case in [7] are possible. Even without the rigorous estimates one readily sees that there is no exponential growth possible, i.e. no modulation instability.

Conclusions
Using the same approach as in the classical Alber equation, we derived a system applicable to crossing seas, i.e. two quasi-unidirectional wave systems meeting in the ocean. This is a genuinely 2-dimensional situation, and we produce for the first time the stability condition that controls whether Landau damping or modulation instability is present. This enables for the first time a systematic detection of modulation instability in crossing seas, including with realistic data. Such a study is now made possible, but involves several tedious steps that are beyond the scope of the current paper.
In this paper we also study a collection of non-parametric unidirectional spectra; this can be thought of as a first step towards the 2-dimensional study. The first subtlety that becomes apparent is that the choice of k 0 , the carrier wavenumber, is an important question that can affect the results in a non-negligible way. (For example, the mean, mode or median wavenumbers are all plausible choices, and we did see spectra where these were significantly different, and would even lead to different stability/instability classifications for a handful of spectra.) Here we used the k 0 provided by the Norwegian Meteorological Institute, and found that the vast majority of spectra in our collection are stable; however among the most extreme sea-states it is possible to find spectra that would be classified as (barely) unstable. This is in agreement with the main findings of [8,7], and indicates that the Penrose-Alber instability really is a limiting factor of how narrow a spectrum can be in the ocean. The fact that there exist spectra all the way up to instability places additional emphasis on the question of how would unstable (or even borderline stable) spectra behave in the field, something that is really not well understood at all.
Moreover, we find that the novel Proximity to Instability (PTI) metric, defined in equation (19), correlates well (ą 85% Spearman rank correlation) with the Benjamin-Feir Index (BFI) and the steepness of a sea state, as well as with Monte Carlo estimates for the kurtosis and the probability of rogue waves. This validates the PTI as a meaningful metric for the study of sea states in the unidirectional setting, where things are rather well understood.
Given that now we have a version of the Alber equation (and the corresponding stability condition) applicable to crossing seas, the natural next step is to investigate crossing seas situations and see how well a generalized version of the PTI metric would correlate e.g. with Monte Carlo estimates of the probability of rogue waves appearing there, as well as directional generalizations of the BFI [55,36]. Another option for this kind of analysis would be by working with the stability condition for the broadband Crawford-Saffman-Yuen equation (CSY) [18,4].