Non-stationary analysis of rogue wave probability over a shoal

Saulo Mendes, 2, ∗ Alberto Scotti, † Maura Brunetti, 2, ‡ and Jerome Kasparian 2, § Group of Applied Physics, University of Geneva, Chemin de Pinchat 22, 1211 Geneva 4, Switzerland Institute for Environmental Sciences, University of Geneva, Boulevard Carl-Vogt 66, 1211 Geneva, Switzerland Department of Marine Sciences, University of North Carolina at Chapel Hill, Chapel Hill, NC, 27599 USA (Dated: June 1, 2021) (Non-peer-reviewed manuscript submitted to EarthArXiv)


I. INTRODUCTION
Ocean statistics offers numerous applications, particularly in marine and offshore safety [1]. Models for short-term and long-term statistics of water waves are used to define the operating envelope for ocean vessels and fixed offshore structures, respectively. Furthermore, understanding the mechanisms of formation of rogue waves has received a considerable amount of attention in past decades [2,3]. Defined as waves at least twice as tall as the significant wave height [2,4], rogue waves present a looming danger to offshore operations [5,6]. The Benjamin and Feir instability [7] arising in surface gravity waves, described by the nonlinear Schrödinger equation [3,8], as well as the linear theory of quasi-determinism [9,10] are the two main theories for the generation of rogue waves. However, there is ongoing debate on which theory can provide the best description of rogue waves [11] and whether they can be unified [12]. From a statistical point of view, the success of either theory is seen in how well the tail of the wave amplitude probability distribution is reproduced.
Longuet-Higgins [13] was the first to attempt to demonstrate theoretically that ocean wave envelopes obey a Rayleigh distribution as a consequence of the central limit theorem. In essence, Longuet-Higgins applied methods and ideas from signal processing [14] to oceanography [15]. In particular, his approach took Rice [14] underlying assumptions of stationarity and ergodicity for granted. Therefore, the resulting distribution cannot account for the influence of sea state parameters [16][17][18]. For similar reasons, higher-order analytical distributions share this inability [19,20]. Though the standard approach has reached considerable success in explaining the observed directional spectrum and wave properties [21][22][23], according to Donelan et al. [24], the need for the ergodicity and stationarity assumptions essentially prevent the use of spectral analysis techniques in unsteady conditions or during isolated events, such as rogue waves. Furthermore, Haver and Andersen [25] points out that whether rogue waves are rare realisations of a Gaussian population or regular realisations of a non-Gaussian set can only be answered through the solution of the governing equations or by relaxing the assumption of stationarity. Nevertheless, we still lack a clear direction and framework to modify linear and nonlinear mechanisms into an unsteady setting.
Following Trulsen et al. [26] laboratory experiments, considerable attention has been given to the shoaling effect on rogue wave formation. Waves approaching a sudden change in bathymetry display substantial non-equilibrium features [see 27, and references therein], thus providing an ideal configuration to probe out-of-equilibrium theories. Over the past decade, additional experiments [28][29][30][31][32][33] and numerical studies [34][35][36][37] have attested heavier tails than expected by Longuet-Higgins [13]. Most studies conclude that nonlinearities and abrupt depth change lead waves out-of-equilibrium, deviating from Gaussian statistics [see 38-40, for the most up-to-date theoretical analysis]. However, analytical or numerical closed-form probability distributions able to describe these laboratory results are still lacking. Towards that end, the present work seeks to create a probability distribution for non-stationary (unsteady) conditions encountered by waves undergoing abrupt depth change. More precisely, we define a functional applied to the probability distribution that transforms the initial wave statistics into a transient one, guided by the reshaping of the energetics due to Stokes [41] perturbation, delineated in sections IV A-IV B. We shall demonstrate that modifying the water wave solution affects Khintchine [42] theorem in unsteady conditions typically up to 10%. Therefore, we lay down a procedure that brings the Trulsen [27] and Haver and Andersen [25] programs partially into effect.

II. BACKGROUND AND GENERAL FORMULATION
Following Massel [43], it is convenient to express distributions in terms of wave and crest heights normalised by the significant wave height H 1/3 , respectively, as: where H= Z c + Z t is the wave height, Z c is the crest height and Z t the trough depth. Assuming we can approximate H = 2Z c in narrow-banded seas, the Rayleigh exceeding probability distribution reads [13]: with f α denoting the probability density of wave heights. Whenever the narrow-band condition loses validity, the relation between wave crests and wave heights requires the knowledge of the spectrum bandwidth ν [44]. Alternatively, this conversion is described through the vertical asymmetry between crest and wave heights [45]. Hence, whenever it becomes necessary to convert normalised crests into wave heights in broad-banded seas, the measure of asymmetry (depending on both ν and skewness µ 3 ) is defined as follows [20]: with · r denoting the long-term wave record average, much exceeding the wave period.

III. THE ROLE OF ERGODICITY & STATIONARITY
The steady-state statistical theory of random waves depends on the relation between the energetics of the Airy [46] wave solution and the statistical treatment of the random noise [47]. Following the establishment of this connection [14], subsequent works applied it without taking into account the underlying assumptions on stationarity, ergodicity, uniform distribution of phases and narrow-band spectrum. The early successes of this paradigm with the narrowbanded hypothesis relaxed [21][22][23][24] made it reasonable to overlook the implication of the complete set of assumptions for the derivation of probability distributions (see Tayfun and Alkhalidi [48] for a review of methodologies). In order to understand how this connection arises, we start by considering the total energy per width contained within a wavelength [49]: where ρ is the density, g is the gravitational acceleration, λ the wavelength, h the water column depth, x the direction of motion and z the vertical axis so that g = −gẑ, ζ is the sea surface elevation and u i = ∂Φ/∂x i is the i-th velocity component derived from the velocity potential Φ. Assuming ζ/h 1 and ζ ∼ ζ = 0 for the second integral interval in eq. (4) [49], the Airy [46] energy reads: where a is the wave train amplitude. On the other hand, the autocorrelation of the sea surface elevation is ( · denotes the time average): such that one can find [47,50]: Then, using Khintchine [42] theorem, which relates the spectral density of a time series with its autocorrelation in stationary processes, one concludes: with S(ω) denoting the unidirectional ocean energy spectrum. However, by definition (see eq. (A1)), the distribution for the surface elevation ζ has to be extracted from its ensemble moments. Unless one assumes a uniform distribution of phases (see appendix A), the above rationale is unable to provide such moments. The above paradigm provides exact values for the time moments and complex physical systems have been typically treated as ergodic since Boltzmann [51]. Therefore, as demonstrated by appendix A, a sufficiently large time series with uniform distribution of phases leads to: thus narrowing down the possible solutions into the Gaussian distribution of the surface elevation because of the ensemble moments E[ζ 2n+1 ] = 0 for all n ∈ N and vanishing excess kurtosis. Despite the simplification brought by ergodic theory in physical systems of high complexity, proving its validity is always a challenging task [52]. Moreover, the stationarity of the ensemble average is not easily verifiable [47]. Hence, wave statistics models have to assume stationarity and ergodicity in order to find results similar to eq. (9). In the next section we shall challenge the assumption of stationarity and pave the way for the analysis of non-equilibrium wave statistics.

IV. NON-EQUILIBRIUM WAVE STATISTICS
Despite the usefulness of the Airy [46] formulation for the spectral analysis of water waves, the evolution of the ocean surface is neither stationary nor ergodic [47,53]. Therefore, if one still uses the Khintchine [42] theorem, higher-order corrections to the probability density will be overlooked. In the spirit of Das and Nason [54], we take into account the non-equilibrium condition by considering a numerical deviation from stationarity. Let us assume that a wave train travelling over a flat bottom has its statistical properties and stationarity well described by linear wave theory. If this group of waves is affected from a sudden change in its physical variables due to the bathymetry, we expect the Khintchine [42] theorem in eq. (8) to be slightly corrected by a factor Γ during the shoaling: hereby considering that the system can still be treated similarly to an ergodic one in spite of its non-stationarity. This approximation relies on the assumption that the distribution of phases is not too far from uniform, as in Tayfun [55]. It is furthermore supported by the fact that the first-moment correction due to Γ is much smaller than the second one (see appendix A). The consequences of including the Γ correction factor have not been fully considered to date, although studies such as Longuet-Higgins [44] relates his distribution to the validity of the Khintchine [42] theorem, and Tayfun [56] considers it a small correction, |Γ − 1| ≈ O(g −2 ). Let us for now focus on how the wave statistics will adapt to an unsteady correction parameter Γ in the ocean. To first order, the probability density of wave heights must fulfil the narrow-band identity: Due to the difficulty to convert surface elevation distributions into crest and height distributions in broad-banded seas [57], eq. (11) replaces the typical envelope approach to find the wave height distribution directly. In fact, analytical models obtained from physical principles typically provide divergent methodologies for the derivation for surface elevation, wave height and wave crest distributions. For instance, while Tayfun [56] derives closed-form expressions for crest and troughs based on a simple change of variables from the surface elevation, soon after Tayfun [58] derives a wave height distribution based on autocorrelation kernels for the surface elevation distribution and requires a Gauss-Chebyshev quadrature for its conversion to wave heights. In spite of that, we can alternatively start with the narrow-band integral of eq. (11) and include broad-banded effects into Γ through the apparatus developed in eq. (3) later in section IV C. Therefore, neglecting skewness µ 3 and kurtosis µ 4 , the change of variables m 0 → Γm 0 results in the narrow-banded correction: that is, the initial wave train Gaussian statistics will be affected by bathymetry and its exceeding probability R α will be transformed into R α,Γ during a sufficiently unsteady shoaling process. We shall proceed and probe the concrete effect of such correction in unsteady settings.

A. Stokes Model Statistics in Unsteady Conditions
Following the same energetics argument as in section III, let us generalise the velocity potential and surface elevation solutions without secular terms [59]: with the auxiliary variables ϕ = k(z + h) and mφ = mk(x − c m t + θ m ) where c m = c m (mk) is the phase velocity of the m-th spectral component. Subtracting the fixed depth h 2 0 /2 term from the energy and fulfilling the conditions ζ/h 1 and ζ ≈ 0 for the integral interval in eq.(4), one can prove for a depth variation ∂h(x)/∂x 1 (see appendix B): The coefficients (Ω m ,Ω m ) can be decomposed into a series of even powers of steepness coupled to factored out trigonometric functions (χ m ,χ m ) (see appendix B), leading to: Hence, by means of the change of the second moment in eq. (10) in comparison to eq. (8), the probability correction due to unsteady conditions reads: This expression demonstrates that the energetics is reduced to the coefficients of eq. (15). We assume that waves, initially propagating on a flat bottom, follow the Airy solution distribution. Afterwards, due to the bathymetry change, Stokes corrections become relevant since a much larger steepness is to be taken into account (see Eagleson [60]). Consequently, out-of-equilibrium dynamics will deform an initially Gaussian statistics due to higher-order effects in steepness. In other words, we are modelling a transition between Airy (before the shoal) and Stokes distributions (during the shoal). In order to generalise the initial distribution, appendix C shows that modelling any bathymetry change by a Airy -Stokes transition is very effective. Surprisingly, the Γ correction only depends on ε and kh, but not on the slope steepness, provided the latter is sufficiently high, i.e. ∂h(x)/∂x 1/20 [37], ensuring non-stationarity into the system in support of eq. (12). Indeed, the experiments in Trulsen et al. [33] had an approximate slope of tan 15 • = 1/3.8. Figures 6-8 of Gramstad et al. [35] showed that a 21% decrease in k p h led to a 320% increase in the amplification of the probability of rogue waves as compared with the pre-shoal conditions, whereas a 400% increase in the slope increased the amplification by 150%. This suggests that for the same magnitude of variation, i.e. 10% decrease (k p h) or increase (slope), will lead to a respective 100% and 5% increase in the amplification due to the shoal. Hence, we conclude that the slope is necessary to maintain the system out-of-equilibrium but the parameters in Γ control the probability amplification.

B. Stokes Second-Order Statistics
The energy per width for a second-order Stokes [41] perturbation in the narrow-banded case (S 0 = 1) reads (see eq. (B9)): with coefficients reading (see figure 1a): The condition for the validity of the model is governed by the Ursell number: Ur = Hλ 2 /h 3 = ε(2π/kh) 3 8π 2 /3 [49,61]. Hence, for small amplitudes (ζ h) we have: whereas the surface total variance reads: The previous equations have been derived for a single wave for simplicity. In order to consider the full wave distribution, we employ the wave number at the peak of the spectrum, k p , in the reduced water depth k p h, and the zero-crossing wavelength λ 2 in the significant steepness ε = H 1/3 /λ 2 [43]. Thus, in the limit of a large number of wave components, the previous expressions read The narrow-banded correction for a group of waves over a changing bathymetry is: Γ exceeds 1 for all depths, with a maximum of up to 1.13 in intermediate depths (k p h ∼ 0.5 − 1) and an asymptotic behaviour for deep water (see figure 1b). This shape defines three regimes, as marked on figure 1b. In shallow water (Regime I), a shoal reducing the depth will reduce Γ, hence the exceeding probability. Conversely, beyond the maximum of Γ (Regime II), the shoal will increase Γ and the exceeding probability. Finally, as long as the shoal  stays within the third regime, the depth variation will translate into a negligible change of Γ, hence will have no consequence on the exceeding probability. The latter regime is compatible with the second-order height distribution in deep water (see appendix D). This behaviour will allow to simplify investigations of shoals starting in deep water (k p h 2), well within Regime III. We can without loss of generality start the analysis of the wave statistics evolution at the point when it enters Regime II.
In order to generalise the derivation of eq. (22) to broad-band seas, we follow the alternative approach in eq. (3). Consequently, the steepness in eqs. (17)(18)(19)(20)(21)(22) will be corrected by the vertical asymmetry S 0 , which in turn modifies the correction parameter: The vertical asymmetry will increase the correction Γ S0 as compared to the narrow-band case by a few percent, as follows (see figure 1b): On the other hand, the parameter Γ must be corrected for wave breaking, leading to slightly smaller peaks than in figure 1b. We include a depth-dependent breaking limit of regular waves [62] by setting ε (ε 0 /7) tanh k p h with 0 ε 0 1: where even if the ratio ε 0 is constant, the actual (breaking-limited) steepness ε will drop considerably throughout the transition from deep to shallow waters (see figure 2a). The vertical asymmetry S 0 increases the value of Γ significantly (figure 2b), but less than the decrease in Γ due to wave breaking. Although eq. (24) typically increases the shoal correction by a few percent, it also shows that the correction has an upper bound Γ 1.20 when ε 0 = 1 and k p h = 0.5 and limiting asymmetry S 0 = 2 (skewed sea with η 1/3 ≈ 2.1). Therefore, a decrease in the depth k p h as reviewed in Trulsen et al. [33] will increase the significant steepness as well as Γ, hence R α,Γ given in eq. (12), i.e. the probability of rogue waves as compared to the pre-shoal, stationary Rayleigh distribution if the shoal reduces the depth toward values in Regime II.

C. The Normalised Height Diagram
Among many of Longuet-Higgins's seminal ideas, the demonstration that the rms surface elevation is connected to the spectral bandwidth ν is usually overlooked. While Longuet-Higgins [44] showed that a bandwidth correction to the standard measure H 1/3 = 4σ will affect the Rayleigh distribution, the converse can in principle be true as well. Goda [63] was among the first to point to the relevance of this ratio as a valuable test for the Rayleigh distribution.
Hence, the Rayleigh or any arbitrary distribution of wave heights contain previously unnoticed information on the ratio H 1/3 /σ, as follows: which shows that at fixed energy (fixed m 0 = σ 2 ) the change in the ratio H 1/3 /σ is balanced by the normalisation α = H/H 1/3 , therefore the overall exponent stays invariant, and thus the wave statistics. However, if and only if m 0 increases in comparison to Airy's solution, we expect the coefficient attached to α 2 to be smaller. In fact, this can be seen by the ratio of the energies at first-order (5) and second-order (17). To investigate the effect of broad-band on H 1/3 /σ by the direct measurement of the vertical asymmetry as in eq. (3), we combine eqs. (12) and (26) and find (see figure 3a): Asymmetry contributes to this expression via two different processes. On one hand, the 1/S 0 factor corresponds to its direct influence on the significant wave height. This factor does not depend on the occurrence of a shoal, although the asymmetry ultimately depends on the depth. Increasing the bandwidth and therefore S 0 will also lower H 1/3 /σ (see eq. (27)), as demonstrated by Vandever et al. [64]. On the second hand, the asymmetry influences the probability amplification by the shoal, via its impact on Γ S0 . Except in very shallow water (Regime I), both S 0 and Γ S0 decrease monotonically as a function of k p h, as we model the pre-shoal zone as a stationary Airy solution concurrently leading to S 0 → 1 and Γ S0 → 1. We can therefore define a parameterisation mapping the latter on the former, in the form S 0 = Γ κ S0 , or, equivalently, its linear version S 0 = 1 + κ(Γ S0 − 1). Note that this parameterisation aims at describing the simultaneous evolution of Γ and the asymmetry when the depth evolves, while eq. (24) only describes the direct impact of the asymmetry on Γ. Since Regime II is the one where we seek the evolution of the exceeding probability, the κ parameterisation appears suitable to our work. Furthermore, its use avoids the numerical issues that would arise if handling S 0 and Γ as independent. Consequently, the ratios √ m0 obtained from our model (figure 3) agree with the asymptotic values of 4 and 3.8, as reported by Goda [63], in the shallow-and deep-water limits, respectively. Moreover, as described in Vandever et al. [64], initial shallow regions (0.2 < k p h < 0.4) can still have a ratio slightly lower than 4 if the bandwidth is very large (1 < ν < 2). Since the typical Pierson and Moskowitz [23] spectrum features ν ≈ 0.45 [65], we conclude that our estimates are in match of both Goda [63] and Vandever et al. [64], as demonstrated by figure 3a.

D. Probability Evolution in the Broad-Banded Formulation
So far we have studied a pre-shoal Rayleigh distribution of the surface elevation. Following eq. (12), the Γ correction due to abrupt bathymetry reads in this case: However, the rationale of section IV B is general. It applies in fact to any initial surface elevation distribution P.
Taking into account the first-order models of appendix C, one can prove that the correction is the same regardless of P. In a narrow-banded sea: where P α denotes the exceeding probability at equilibrium prior to the shoal and P α , Γ is the distribution within the non-equilibrium zone. As a consequence, the correction term to the exceeding probability due to the bathymetry change is e 2α 2 (1− 1 Γ ) , regardless of the initial probability distribution function. Hence, the highest possible realistic values of Γ ∼ 1.08 (see figure 1) lead to a two-fold increase of rogue wave probability (α 2).
In the case of a broad-banded sea, the connection between crest and height statistics responsible for eq. (27) can be reinterpreted as a simple change of variables Γ → S 2 0 Γ S0 in eq. (12), so that eq. (29) becomes: To ensure the numerical stability of the distribution along the propagation, we use the κ parameterisation introduced in section IV C. For practical purposes, we estimate it as: where the maximum of Γ is taken over the propagation on the whole shoal, ε being the average of the pre-and postshoal steepness. It can be measured at both stages of the wave propagation, or, alternately, obtained from forecast or hindcast, from the mean period and significant wave height. Note that uncertainties in asymmetry marginally impact the value of κ 0 , as errors in the numerator and in the expression of Γ tend to partially compensate each other due to eq. (24). Consequently, the effect of the shoal on the exceeding probability evolves during the propagation along x as follows: When the evolution takes part in two different regimes, the above equation is necessary as a truncation to ensure stability and consistency. Otherwise, a continuous κ restricted to Regime II is applicable. Note that under adequate conditions, a more compact form can be derived, as detailed in appendix E.

V. COMPARISON WITH TRULSEN EXPERIMENTS
Raustøl [28] provides experiments of wave propagation over a shoal, later summarised in Trulsen et al. [33]. Remarkably, the probability distribution of rogue waves evolving with the distance from the wavemaker was recorded, offering a benchmark for our non-stationary correction to the wave height probability distribution. Since our model is expressed in terms of relative depth and significant steepness, we have extracted the raw data from Figures 5.4-5.5 of Raustøl [28], in accordance with the inversion k p h ≈ (k p a c /Ur) 1/3 of Trulsen et al. [33]. In order to facilitate the processing, we have smoothed the experimental data by fitting analytic functions on the data points, as described in detail in appendix F.

A. Results
The experiments feature a shoal starting at x = 0, rising up to 42 cm at x = 1.6 m, followed by a plateau till x = 3.2 m and a decay to zero until x = 4.8 m. The initial depth ranges between 50 and 60 cm depending on the runs, as detailed in table 1 of Trulsen et al. [33] and table 5.2 of Raustøl [28]. The values of κ 0 calculated according to eq. (31) amount to (5.9, 4.9, 4.3, 3.6, 2.9, 2.6, 2.3, 2.4, 2.3, 2.2) for the runs 1-2, 4-9 and 11-12, respectively. Based on the κ 0 parameterisation described in the previous section, we compute the evolution of the rogue-wave probability as a function of distance for each run, assuming a pre-shoal stationary Rayleigh distribution (figure 4). Our model that takes into account the evolution of the wave asymmetry S 0 over the shoal (cyan curves) due to the skewness of the surface elevation distribution (eq. (3)) reproduces well the experimental data, over the whole shoaling episode and for all runs. Disregarding the evolution of skewness reported in Trulsen et al. [33] while keeping a vertical fixed asymmetry S 0 = 1.2 (blue curves in figure 4) degrades the agreement only marginally, although the probability rises slightly earlier, decays slightly later, and the asymmetry between the up-and down-shoaling phases is reduced. Furthermore, we point out the remarkable difference of amplification between vertically asymmetrical (solid blue, cyan) and symmetrical seas (dashed red). This happens because eq. (28) leads to a maximal amplification between 75-100% for Γ ≈ 1.08 − 1.10. When we include the typical vertical asymmetry S 2 0 ∼ 1.5, the pre-shoal probability P α is transformed into P 2/3 α ∼ 10P α within Regime II, seemingly becoming an alternative to Gram-Charlier models. Finally, adjusting the stationary pre-shoal probability to the observed values instead of considering an initial Rayleigh distribution (dotted curve) improves the agreement (runs 6-12, panels e-j), demonstrating the amplification universality. This agreement over the whole range of experimental conditions reported by Raustøl [28] is remarkable, as it requires no specific parameter tweaking.

B. Discussion
In contrast, the early model of Tayfun [56] (see its eqs. (36)-(38)) prescribed no second-order effect for wave heights, also discussed in Tayfun [66], Tayfun and Fedele [67]. In fact, these models fall under the broad category of quasideterminism theories [50] of which Longuet-Higgins [44] and Naess [68] also take part and preclude wave heights from exceeding the Rayleigh distribution, typically being lower than the latter as the bandwidth broadens. Therefore, these formulations in stationary conditions would not be able to describe the Raustøl [28], Trulsen et al. [33] experiments. Also according to these models, any departure from a Rayleigh distribution of wave heights is due to third-order nonlinearities [67,69]. Nevertheless, it is important to remark that most of these theories were devised for deep waters, hence they fall within the Regime III of our model. Moreover, following Marthinsen [70], Dingemans [71] one can tentatively elaborate alternative finite-depth second-order mathematical structures as in appendix D. Comparing with our model, we observe a numerical equivalency in deep water, but not in transitional and shallow waters where the alternative second-order models display a sharp departure from Goda [63] observations. Unfortunately, the finite depth model of Tayfun and Alkhalidi [48] is not a distribution of wave heights, hence not applicable to our discussion. Our Γ-model of the probability evolution over a shoal is based on the Stokes perturbation, which is valid for Ur 8π 2 /3 [49]. This limit of validity is well beyond the Ursell number Ur 0.22 of the experiments considered in this work [28,33]. We can therefore expect that the Γ-model will still apply to more than 10 times larger waves and/or shallower waters.
The transient drop of P α , Γ /P α to almost zero in the rising region of the shoal (0 x ≤ 1.6 m) in runs 1, 2, 4, and 5, as well as its slow decay on the trailing side (3.2 ≤ x ≤ 4.8 m) could be due to higher-order nonlinear effects not captured in our model, including the evolution of the skewness and kurtosis of the surface elevation probability distribution when propagating on the shoal.
Recent developments in Zhang and Benoit [40] described the exceeding probability of Run 3 of Trulsen et al. [33] throughout the shoaling and de-shoaling stages by means of numerical simulations. However, a closed-form expression for the probability distribution was not provided, and cases with low or vanishing skewness and kurtosis like Run 12 were not addressed. Our model, possessing no free fitting parameter and directly obtained from Stokes perturbation, has shown to be valid for all runs in Raustøl [28], thus describing that the intensity of the anomalous wave statistics decreases from Run 1 to Run 12 almost monotonically. This evolution of the rogue wave probability could be related to the process proposed by Li et al. [38] and experimentally confirmed in Li et al. [39], in which the generation of additional wave packets that interact with the original pre-shoal wave packet under abrupt depth changes leads to a local peak some distance into the shallower region. Although we were not able to analyse Run 3 [72] due to the lack of wave height probability data, we are confident that our model can reproduce it well, as it does for the very similar Runs 1, 2, and 4. Furthermore, our model and in particular the three regimes discussed on eq. (22) and figure 1b also allow to understand the contradictory behaviours highlighted by Trulsen et al. [33] and Zhang and Benoit [40]. While the rogue wave probability is not affected by a shoal in initially deep water (Regime III), it does increase on a shoal in intermediate depth (Regime II). Interestingly, a third regime (Regime I of figure 1b) is consistent with observed surf zone statistics over less abrupt depth changes [73], as it ascribes likelihoods to rogue waves in shallow water to be smaller than in transitional waters. In addition, Barbariol et al. [74] shows an increase of maximum crest up to some cut-off in relative depth, upon which it starts to sharply decrease the higher the ratio H 1/3 /h becomes, thus providing support for Regime I in our model. However, numerical simulations and experiments with even shallower shoals are necessary for a conclusive assessment. Therefore, we have demonstrated that rogue wave statistics will be enhanced by non-equilibrium dynamics of abrupt depth change at a suitable interval of relative depth 0.5 < k p h < 1.5 but cannot be amplified endlessly. It is therefore less pronounced past the threshold k p h 0.3 until it starts to follow shallow water distributions such as Glukhovskii [73]. of Raustøl [28]. Blue: Model of eq.(29) for a pre-shoal (x < 0) Rayleigh distribution with constant S 0 = 1.2; Cyan: Same model considering the evolution of skewness over the shoal; Dotted: Model of eq.(29) for a pre-shoal probability matched to the experimental data. Dashed: same as blue but for symmetrical seas S 0 = 1.

VI. CONCLUSION
This work presented a connection between statistical distributions and the fluid mechanics of the second-order Stokes perturbation in non-stationary conditions, successfully providing a physical explanation for the rogue wave probability increase under a depth change [28,33]. To confirm the observation of maximal rogue wave enhancement atop the shoal [28,33], we have shown that our model reproduces very well the experiments of Raustøl [28] regarding the probability distribution as a function of the distance from the wavemaker. Moreover, we showed which physical variables could affect the validity and numerically assessed the extent of the deviation from the assumption of stationarity. Furthermore, instead of introducing new physics [25], our model has demonstrated that an effective theory arises by challenging the stationarity assumption. The latter distinguishes regimes of probability amplification and provides a new interpretation for the H 1/3 /σ ratio [63] as a guide to estimate probabilities in a transition from linear to second-order wave theory. Consequently, we have established that the deep water regime is expected to not produce any significant amplification of the height distribution, whereas the transitional water within 0.5 < k p h < 1.5 provides strong amplification, and the shallow water regime will tend to decrease this large amplification to a level smaller than the initial stage in deep water, as shown in figure 3.
Our model has been restricted to reformulating moments up to the variance only. The evolution of either skewness or kurtosis as a function of the distance from the wavemaker [33] shall be addressed in a subsequent work. On the other hand, we have shown that it is possible to fit the data without the application of either skewness nor kurtosis, as they are 'symptoms' of the dynamics and not the cause [18,75,76]. Since our model has relied on relations between parameters that might be affected by the abrupt depth change, namely the dependence of the probability amplification on steepness, slope and bandwidth, the empirical findings in Mendes et al. [20] regarding vertical asymmetry have to be extended to the current setting for an exact formulation. Therefore, future work has to investigate how the connection between vertical asymmetry, nonlinearity and skewness differ from steady to unsteady conditions towards a full treatment of this problem. Furthermore, the generalisation of this work to multidirectional spectra would be needed, since Ducrozet and Gouin [36] suggest that such configuration weakens the effect of a varying bathymetry. Whether rogue waves are enhanced in strong bathymetry changes throughout most oceans or regionally under suitable conditions is yet to be assessed. The exact computation of the ensemble average of the sea surface elevation at an instant of time t 0 is not trivial [47]. Without ergodicity, the probability density of the surface elevation is unknown unless one assumes its expected value: where dµ(ζ) is a measure on the space of possible surface elevations. For an oscillatory system, it is customary to write ζ = ξ cos(φ) and introduce a joint PDF on the space [0, ∞) × [0, 2π) that describes the statistical distribution of amplitudes and phases. Without loss of generality, we can always choose units in which E[ζ] = 1, and write: Longuet-Higgins [13] distribution is recovered if we assume that the phases are uniformly distributed and uncorrelated from the amplitudes, that is, f (ξ, φ) = 1 [3,14,77,78]. However, the uniform distribution of phases is only appropriate for narrow-banded signals [79,80]. For more realistic sea states, Tayfun [55] shows that the PDF introduced above should be corrected to account for correlations between phase and amplitude. As a first approximation, the correction proposed by Tayfun [55] is: where µ 3 is the skewness. Note that within this approximation, the expected value of even order powers of ζ are not modified relative to the uniform phase approximation. Odd order powers, which are zero when f = 1, are now non zero. On the other hand, it is known that the weakly nonlinear evolution of a sea state which at t = 0 has random and uniformly distributed independent phases, remains so over the nonlinear time interval t 2π/ω p [81,82]. Therefore, the source of correlations between phases and amplitudes cannot be attributed to the internal weakly nonlinear dynamics. However, this does not preclude that external factors (e.g., wind forcing) inducing nonequilibrium dynamics can nudge the phase distribution away from uniformity and/or impart a correlation between phases and amplitudes.
For the purpose of illustration, let us analyse the simplest effect of a uniform distribution. Through a change of variables [83] and the law of the unconscious statistician [84], we rewrite the ensemble average: Therefore, the uniform distribution of phases leads to ergodicity: If we assume that correlations develop between the phases, a weak departure from ergodicity will be observed. Below we show that, due to the ergodicity assumption, the accuracy of Gaussian statistics will deteriorate the narrower the superposition distribution is. For the sake of measuring appreciable deviations and without loss of generality, we use a Boltzmann-like distribution, such that the ensemble average of the sea surface reads: which is relatively small compared to the second moment of the surface elevation. For a tentative Gaussian shaped superposition, however, one finds: while for the square of the surface elevation we obtain: Lacking a closed-form, the Gaussian-like distribution of phases reads instead: We compare the deviations from the uniform superposition in eqs. (A7)-(A9): whereas for the Gaussian one we have: implying a decreasing gap between moments when the superposition distribution is narrower.
In these examples, we see that even the linear evolution of a sea state that, at some point in time, has a non uniform distribution of phases breaks ergodicity.
As an immediate corollary, we find 2I mm = cosh (2mϕ) + cos (2mφ). Using the algebra from eq. (B4) and periodic integration, following the expression for the energy in eq. (4) and subtracting the potential energy ρgh 2 0 /2 due to the water column, we find: As we assumed a small effect of ∂h/∂x in the previous integral, the bathymetry will appear in ϕ = k(z + h) as well as in Ω m andΩ m . Likewise, using the definition of eq. (13), we compute the time average of the squared sea surface elevation: To compare the generalised model with the specific case of Airy's solution we set m = 1. In this caseΩ 1 = a while Ω 1 = 1 · kf 1 = agk/ω cosh kh [71], the dispersion relation is expressed as ω 2 = gk tanh kh, so that the energy of a single wave is: thus recovering the energy in eq. (5). For the second-order we have instead: ; Ω 2 = 3ka 2 ω 4 sinh 4 kh ;Ω 1 = a ;Ω 2 = ka 2 cosh kh 4 sinh 3 kh 2 + cosh (2kh) .
Upon the steepness being expressed as ka= (2π/λ) · (H/2) = πε, the energy is computed for an individual wave, leading to eqs. (18)- (19): while the sub-Rayleigh second moment is expressed as, and the super-Rayleigh as well (see figure 5), Comparing the variances with the exponent in eq. (C1), we find µ ± 2 = Γ(1±L µ /Γ), and the corrections L µ can be readily available by isolating the quotients inside the brackets of the r.h.s in the above equations. These coefficients can be further approximated as follows (recalling that |Γ − 1| 1): Then, we showed a weak proof of the first part of eq. (29) by demonstrating that a change g ± µα in the pre-shoal Rayleigh probability will always be met by a change ±L µ in the variance, typically obeying L µ /Γ 1 (see figure 5b). We have also corroborated appendix A on the negligibility of the µ 1 amplification in comparison with the variance.

Generalised Proof
In this section we use the Γ-GC model to obtain exact closed-form wave height distributions, following the steps for the integration of the envelope in a two-dimensional random walk of Mori and Yasuda [78]. Thus we integrate the joint distribution of both surface elevation ζ and its Hilbert transformζ over the uniform distribution of phases, obtaining the marginal density of the surface envelope: where ζ 2 +ζ 2 = R is the envelope with ζ = R cos φ andζ = R sin φ. Performing this integration, changing variables to wave heights and later normalising by H 1/3 = 4m 0 = 4 and integrating again as in eq. (2), we find (see figures 6a,b): Then, in analogy with eq. (28) we are able to assess whether or not the proposition in eq. (C1) holds by verifying the following ratio: Accordingly, figures 6c,d display the magnitude of the correction 1 ± L µ /Γ in the variance. When Γ≈1.15 we see that super-Rayleigh distributions have a maximal 4% increase in the variance Γ, whereas sub-Rayleigh distributions exhibit the opposite but of smaller magnitude, confirming the estimates in the weak proof. As the validity of eq. (C1) has been demonstrated, one can prove the universality of the amplification regardless of the distribution, i.e. extend eq. (28) to an arbitrary distribution. Having in mind the order of magnitude of L µ /Γ, we can rewrite eq. (C1) as follows (defining | ln g ± µα | = ln g µα ): Furthermore, the relative probability becomes: As |Γ − 1| 1, it is straightforward to see that δ µΓ =L µ /Γ (Γ − 1) ∼ 10L µ ∼ δ µα = (ln g µα )/2α 2 , plotted for comparison in figure 7. Hence, we conclude that for the probability amplification the second correction term δ µΓ counters δ µα . Thus, we have proved the first-order amplification universality with µ 3 ∼ µ 4 < 1: Rogue waves in eq. (30) have 2α 2 1 − (S 2 0 Γ) −1 3, resulting in maximal correction e O(Γ−1) ≈ e 3×0.03 < 1.1, suggesting an upper bound of 8-9% variation from the universal amplification. On the other hand, figure 7 shows (1 ∓ δ µα )(1 ± δ µΓ )≈ 1.14 for ordinary waves (α 1) instead of (1 ∓ δ µα )(1 ± δ µΓ ) 1.04 for rogue waves, such that the main term of eq. (30) becomes 2α 2 1 − (S 2 0 Γ) −1 1, and the first-order correction reads e O(Γ−1) ≈ e 0.6×0.15 < 1.1. Hence, the bound is upheld by any normalised wave height. which normalised by the variance reduces to (with nomenclature σX := X): Applying the finite depth coefficients in eq. (2.382a) of Dingemans [71], one finds: leading to the mathematical structure (see Mendes et al. [20] for the F = 1 case): Note, however, that this distribution was not derived in the narrow-band model of Tayfun [56], it is rather an adaptation to extract the mathematical structure of a probability growing boundlessly with increasing steepness. Although Tayfun [87] puts limitations to the steepness as µ 3 kH 1/3 /(1 + ν 2 ), North Sea data [18,20] shows that the skewness can grow higher. Therefore, the suggested adaptation to wave heights is usefull since now ε can reach its breaking limit [62]. In figure 8a we compare our Γ model of eq. (24) with that of eq. (D4) through the H − σ diagram:  figure 4a for a higher reference Γ correction due to steepness ε + ≈ 1.16 ε (dashed) and its original description according to eq. (32) (solid). Notice that this 16% increase in the reference steepness decreases the blue curve model by 4% and the cyan one by 2%. (b) Plot of Γ correction for runs 1-4 of Raustøl [28] as a function of relative depth (dashed) and distance from the wavemaker (solid) with variables (ε, kh) modelled by appendix F, whereas the minimum threshold applicable (Γ 1.01) representative of Regime II is depicted by the thin horizontal line. The averages over these ranges approximately read Γ(x) = 1.036 and Γ(k p h) = 1.023.
The above relation holds because the typical difference between the Γ correction atop the shoal and the average moving under the shoal does not exceed δΓ 0 /Γ 0 2%. We corroborate this by comparing the equivalent of figure 4 (panel a) for a 0.5% higher Γ correction than its reference Γ 0 = 1.031 in figure 9a. Considering that Run 1 had maximal correction Γ ≈ 1.041, the choice for the reference κ 0 does not affect the validity of eq. (32). To obtain a compact formulation of eq. (32), we notice that at a fixed point in space, i.e. at a distance x = x i from the wavemaker, the asymmetry S 0 does not depend on how we plot Γ. Therefore, we can approximate As the two versions of the Γ correction differ, as shown in figure 9b, one concludes that κ(k p h) ≡ κ 0 = κ x . Thus we obtain: whose Regime II restriction is numerically translated to Γ 1 + O(g −2 ), leading to κ x ≈ 0.64 κ kh , as estimated from figure 9b. Therefore, we can conservatively estimate throughout the entire trajectory:  [28] experiments according to eq. (F1), corrected to the term ε = H 1/3 /λ 2 = 4 π k p a c . Dots represent data extracted from Figure 5.4 of Raustøl [28]. Vertical dashed lines depict the rising shoal end (x = 1.6) and beginning of the descending shoal (x = 3.2).

(E6)
Appendix F: Analytical Description of Steepness and Depth In order to smooth as well as to facilitate the handling of the experimental data on wave steepness and water depth, we fitted them against analytic functions, as follows: while the modelling for the depth is computed as, The values of these coefficients for ten runs of Raustøl [28], Trulsen et al. [33] are given in table I. As displayed in figure 10, these fits provide an accurate description of the actual data. Re-scaled by π/4, the first steepness coefficient ε 1 is equal to the pre-shoal steepness in table 1 of Trulsen et al. [33] while ε 1 + ε 2 is the shallower steepness of the shoal. Likewise, the coefficient D 1 , which was extracted from Raustøl [28] through the formula k p h ≈ (πε/4Ur) 1/3 (see Trulsen et al. [33]), approximately equals the deeper side relative depth while D 1 + D 2 recovers the shallower side. Note, however, that Trulsen et al. [33] displays the averages of each sides, while we model every value according to the 16 probes of figure 2 in Trulsen et al. [33].