Non-homogeneous analysis of rogue wave probability evolution over a shoal

Abstract

Non-equilibrium evolution of wave fields, as occurring over sudden bathymetry variations, can produce rogue seas with anomalous wave statistics. We handle this process by modifying the Rayleigh distribution through the energetics of second-order theory and a non-homogeneous reformulation of the Khintchine theorem. The resulting probability model reproduces the enhanced tail of the probability distribution of unidirectional wave tank experiments. It also describes why the peak of rogue wave probability appears atop the shoal, and explains the influence of depth on variations in peak intensity. Furthermore, we interpret rogue wave likelihoods in finite depth through the $H$–$\sigma$ diagram, allowing a quick prediction for the effects of rapid depth change apart from the probability distribution.

1. Introduction

Ocean statistics offers numerous applications, particularly in marine and offshore safety (Toffoli et al. 2005). Models for short- 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 (Dysthe, Krogstad & Muller 2008; Onorato et al. 2013). Defined as waves at least twice as tall as the significant wave height (Haver 2000; Dysthe et al. 2008), rogue waves present a looming danger to offshore operations (Faulkner & Buckley 1997; Faukner 2002). Theories of rogue wave formation include the Benjamin and Feir instability (Benjamin & Feir 1967) arising in surface gravity waves, described by the nonlinear Schrödinger equation (Zakharov & Ostrovsky 2009; Onorato et al. 2013), and linear mechanisms such as constructive interference (Boccotti 2000; Fedele et al. 2016; Dematteis et al. 2019). From a statistical point of view, any particular theory relies on its ability to reproduce the tail of the wave amplitude probability distribution.

Longuet-Higgins (1952) applied methods and ideas from signal processing (Rice 1945) to oceanography (see St Denis & Pierson (1953) for a review). In particular, his approach took the underlying assumptions of Rice (1945) about homogeneity and ergodicity for granted. Therefore, the resulting non-dimensional Rayleigh distribution of wave heights cannot account for the varying sea state parameters, such as steepness (Stansell 2004) or depth (Glukhovskii 1966). The same limitation applies to higher-order analytical distributions (Karmpadakis, Swan & Christou 2020; Mendes, Scotti & Stansell 2021). Though these standard approaches have had considerable success in explaining the observed directional spectrum and wave properties (Phillips 1958; Hasselmann 1962; Pierson & Moskowitz 1964), the need for the ergodicity and spatial homogeneity assumptions essentially prevents the use of spectral analysis techniques in unsteady conditions or during isolated events, such as rogue waves (Donelan, Drennan & Magnusson 1996). Furthermore, Haver & Andersen (2000) were the first to suggest a link between rogue waves and non-stationarity. Nevertheless, we still lack a framework to account for unsteady conditions, such as shoaling. (For a review of the consequences of non-stationarity/homogeneity, see Appendix A.)

Following the laboratory experiments of Trulsen, Zeng & Gramstad (2012), considerable attention has been given to the shoaling effect on rogue wave formation. Waves approaching a sudden change in bathymetry provide an ideal configuration to probe out-of-equilibrium conditions (Trulsen 2018). Additional experiments (Raustøl 2014; Ma, Ma & Dong 2015; Bolles, Speer & Moore 2019; Zhang et al. 2019; Zou et al. 2019; Trulsen et al. 2020) and numerical studies (Zeng & Trulsen 2012; Gramstad et al. 2013; Ducrozet & Gouin 2017; Zheng et al. 2020; Zhang & Benoit 2021) have attested heavier tails than expected by Longuet-Higgins (1952). Nonlinearities and abrupt depth change lead waves out of equilibrium, deviating from Gaussian statistics. For a step, these bathymetry effects have been described by nonlinear evolution of interacting free modes with a truncation of the Korteweg–De Vries equation (Majda, Moore & Qi 2019; Moore et al. 2020) and as travelling wave packets subject to second-order effects in steepness (Li et al. 2021a,c). However, wave height probability distributions able to describe the laboratory results of Trulsen et al. (2020) are still lacking. In fact, Li et al. (2021b) have obtained a wave crest probability distribution for the step case based on Tayfun (1980), but as it is stated in the section following (36)–(38) of Tayfun (1980), this framework cannot lead to non-Gaussian wave height distributions. The present work seeks a probability distribution for non-homogeneous conditions encountered by waves undergoing rapid depth change. More precisely, we show that regardless of the pre-shoal probability distribution shape, the transformation of the sea surface elevation and spatial energy density through leading second-order effects in steepness will amplify the rogue wave probability. This amplification depends on the dimensionless depth ($k_{p}h$) and significant steepness ($\varepsilon$).

2. Background and general formulation

Following Massel (2017), we focus on wave and crest heights normalised by the significant wave height $H_{1/3}$, as, respectively,

(2.1a,b)\begin{equation} \alpha \equiv \frac{H}{H_{1/3}} = \frac{ \mathcal{Z}_{c} + \mathcal{Z}_{t} }{H_{1/3}},\quad \beta \equiv \frac{ \mathcal{Z}_{c} }{H_{1/3}}, \end{equation}

where $H$$= \mathcal {Z}_{c} + \mathcal {Z}_{t}$ is the crest-to-trough height (hereafter denoted as wave height), $\mathcal {Z}_{c}$ is the adjacent crest height, and $\mathcal {Z}_{t}$ is the trough depth. Assuming that we can approximate $H = 2\mathcal {Z}_{c}$ in narrow-banded seas, the Rayleigh exceedance probability distribution reads (Longuet-Higgins 1952):

(2.2)\begin{equation} \mathbb{P}_{\mathcal{R}}(H>\alpha H_{1/3}) \equiv \mathcal{R}_{\alpha} = \int_{\alpha}^{+\infty} f_{\alpha_{{0}}} \,{\rm d}\alpha_{{0}} = \int_{\alpha}^{+\infty} 4\alpha_{{0}} \, {\rm e}^{{-}2\alpha_{{0}}^{2}} \,{\rm d}\alpha_{{0}} = {\rm e}^{{-}2\alpha^{2}}, \end{equation}

with $f_{\alpha }$ denoting the probability density function of wave heights. For broader spectra, we define the vertical asymmetry (Kjeldsen 1984; Myrhaug & Kjeldsen 1986) between crest and trough as the ratio between crest and wave height (Linfoot, Stansell & Wolfram 2000). Hence whenever it becomes necessary to convert normalised crests into normalised wave heights in broad-banded seas, they are computed as follows (see § 6.1 of Mendes et al. 2021):

(2.3)\begin{align} {\mathfrak{S}_{0} {(\alpha)} \equiv } {2}\left\langle \frac{\beta}{\alpha } \right\rangle_{r} \approx \frac{{2}\eta_{{1/3}}}{1+\eta_{{1/3}}} \left[ 1 + \frac{ 2\eta_{{1/3}}\,\mathrm{Re} \left( \sqrt{\alpha-1} \right) }{7+ 2\,\mathrm{Re} \left( \sqrt{\alpha-1} \right)} \right],\quad {\eta_{1/3} \equiv \left( \frac{\langle \mathcal{Z}_{c} \rangle_{r}}{\langle \mathcal{Z}_{t} \rangle_{r}} \right)_{H>H_{1/3}}}, \end{align}

with $\langle \cdot \rangle _{r}$ denoting a wave record average. Experimental works typically use 20 min records. Remarkably, the vertical asymmetry is related to the skewness $\mu _{3}$ through the approximation $\eta _{1/3} \approx 1 + \langle \mu _{3} \rangle _{r}$ (see (14) and figure 8 of Mendes et al. 2021). Moreover, considering the correlation between asymmetry, skewness and significant steepness reported by Guedes Soares, Cherneva & Antao (2004), the asymmetry is weakly dependent on the bandwidth $\nu$ due to the bound $\mu _{3} (1+\nu ^{2}) \lesssim kH_{1/3}$ (Tayfun 2006).

2.1. Equilibrium wave statistics

In view of the equivalence between the spectral analysis of spatial and time domains (see Appendix A), we consider the average spatial energy density around $x$ (Dean & Dalrymple 1984), calculated over one spectral zero-crossing wavelength ${\bar {\lambda }}$ (Massel 2017):

(2.4)\begin{equation} \mathcal{E} = \frac{\rho }{2 {\bar{\lambda}} } \int_{x}^{x+{\bar{\lambda}}} \left[g \left( \zeta + h {(x)} \right)^{2} + \int_{{-}h {(x)}}^{{ \zeta }} ( {u^{2}_{1}} + {u^{2}_{3}} )\,{\rm d} z\right]{{\rm d}\kern0.06em x}, \end{equation}

where $\rho$ is the density, $g$ is the gravitational acceleration, $h$ is the water column depth, $x$ is the direction of motion, and $z$ is the vertical axis so that $\boldsymbol {g}=-g \hat {z}$, $\zeta$ is the sea surface elevation, and $u_{i} = \partial \varPhi /\partial x_{i}$ is the $i$th velocity component derived from the velocity potential. Indeed, the experiments of Trulsen et al. (2020) had length of slope 1.6 m while the typical peak wavelength ranged from 1 to 1.8 m, thus validating our calculation above. For an irregular wave field obeying the solution of linear theory (Airy 1845) with uncorrelated random phases $\theta _{i}$ and amplitudes $a_{i}$, one has the surface elevation (Tayfun 1980)

(2.5)\begin{equation} \zeta_{1} (x,t) = {\sum_{i}} a_{i} \cos{(k_{i}x - \omega_{i} t + \theta_{i})}, \end{equation}

where the $i$th components have wavenumber $k_{i}$ and frequency $\omega _{i}$. Assuming $\zeta / h \ll 1$ and $\langle \zeta \rangle \approx 0$ for the second integral interval in (2.4) (Dean & Dalrymple 1984), the energy reads

(2.6)\begin{equation} \mathcal{E}_{AIRY} = \tfrac{1}{8} \rho g H^{2} = \tfrac{1}{2} \rho g a^{2},\quad \therefore \quad \mathscr{E}_{AIRY} = \frac{\mathcal{E}_{AIRY}}{\rho g} = \tfrac{1}{2} \sum_{i} a_{i}^{2}, \end{equation}

where $a$ is the wave train amplitude. Then, using the spatial counterpart of the Khintchine (1934) theorem, which relates the spectral density of a spatial series to its autocorrelation in homogeneous processes, one concludes (see Appendix A) that

(2.7)\begin{equation} R_{x} (\xi=0) := \mathbb{E}[\zeta^{2}] =\langle \zeta^{2} \rangle_{{x}} = \int_{0}^{+\infty} S({k})\,{\rm d}{k} \equiv m_{0} = \mathscr{E}_{{AIRY}}, \end{equation}

with $S({k})$ denoting the unidirectional ocean energy spectrum based on the wavenumber, which is equivalent to computing the autocorrelation in time and reformulate it in terms of the spectrum $S(\omega )$ for stationarity in time when the system is both homogeneous and stationary (Massel 2017). This water wave solution features a Rayleigh distribution of wave heights in the form $\mathbb {P}(H > H_{0}) = {\rm e}^{-H_{0}^{2}/8m_{0}}$ in an irregular wave field with narrow-banded spectrum (Longuet-Higgins 1952). In the next section we will challenge the assumption of homogeneity implied by the use of the Khintchine (1934) theorem, paving the way for the analysis of non-equilibrium statistics.

3. Non-equilibrium wave statistics

Despite the usefulness of the Airy (1845) formulation for the spectral analysis of water waves, the evolution of the ocean surface is not stationary, spatially homogeneous or ergodic (Cherneva & Guedes Soares 2008; Goda 2010). Therefore, higher-order (unsteady) corrections to the Khintchine (1934) theorem should be considered. During the shoaling, the autocorrelation function is computed from (Here, the surface elevation $\zeta (x,t)$ has been denoted as $\zeta (x)$ to ease the notation.)

(3.1)\begin{equation} {R_{x} (\xi, x) = \mathbb{E}[\zeta (x )\,\zeta (x+\xi )] = \int_{-\infty}^{+\infty} \zeta (x)\,\zeta (x+\xi)\,f_{\varGamma {(x)}}(\zeta) \,{\rm d}\zeta,} \end{equation}

where $f_{\varGamma {(x)}}(\zeta )$ is the probability density of the surface elevation $\zeta (x,t)$ at a fixed point $x$, and is expected to depend on a correction $\varGamma {(x)}$ due to the effect of bathymetry. In the spirit of Longuet-Higgins (1980) and Das & Nason (2016), the correction is defined by comparing how the ensemble and spatial energy averages change at and past the shoal, as

(3.2)\begin{equation} {\varGamma (x) := \frac{\mathbb{E}[\zeta^{2} {(x,t)}] }{ \mathscr{E} } = \frac{\mathbb{E}[\zeta^{2} {(x,t)}]{(x)} }{ \mathscr{E}(x)}}, \end{equation}

which means that in a non-homogeneous setting, the surface elevation depends on $(x,t)$ but both ensemble and energy averages depend only on space. Clearly, due to (2.7), $\varGamma \rightarrow 1$ both prior to and several wavelengths after the shoal. Without the exact shape of the random phase distribution as in (A7), it is impossible to find the ensemble average ${\mathbb {E}[\zeta ^{2}]}$ because the probability density $f_{\varGamma }(\zeta )$ is unknown: the goal of the present work is to find its wave height counterpart. Since the targeted experiments are stationary in time but non-homogeneous in space, we have that ${\mathbb {E}_{x = x_{i}}[\zeta ^{2}]} = \langle \zeta ^{2} \rangle _{t, x = x_{i}}$, i.e. that the time average equals the ensemble average at any point $x_{i} \in \mathbb {R}$ of the spatial evolution. Because the evolution of a generalised spatiotemporal autocorrelation can be described as evolutionary spectrum $S(\omega, x)$, like those computed for the targeted experiments (Trulsen et al. 2020) whose statistics were generated considering the time series at different locations (Lawrence, Trulsen & Gramstad 2021), we may reduce (3.2) to the spatial evolution of the autocorrelation in time:

(3.3)\begin{equation} { \varGamma (x) = \frac{R_{t} ({\tau = 0, x})}{\mathscr{E}(x)} = \frac{\mathbb{E}[\zeta(x, t)\,\zeta(x, t + 0)] }{ \mathscr{E}(x) } \equiv \frac{\langle \zeta(x, t)\,\zeta(x, t + 0) \rangle_{t} }{ \mathscr{E}(x)}.} \end{equation}

Note that the spatial and temporal averages are functions not of time but of horizontal displacements, much like the ensemble average. Indeed, for homogeneous processes, it is customary to fix a free parameter $A$ such that the area under the spectral curve is equal to the mean power (Massel 2017), e.g. $R_{x}(\xi =0) = \int _{0}^{\infty } S({k}) \,{\rm d}{k}$ for a one-sided spectrum. Instead, we choose $A$ such that the area under the spectral curve matches the spatial energy density ($\mathscr {E}$) during homogeneous stages, and it reduces to the autocorrelation in the Gaussian case, i.e. in the strictly homogeneous case prior to the shoal. For non-homogeneous processes, this is not the case, and the methods for finding a spectrum produce anomalies and undesirable features (Loynes 1968; Cohen 1989; Flandrin 1989; Adak 1995; Bruscato & Toloi 2004). This implies that there is no canonical or unique way to define a non-homogeneous Khintchine theorem (Flandrin 1989) such that the ratio ${\mathbb {E}[\zeta ^{2}]} / \mathscr {E}$ computes the deviation from homogeneity. On the other hand, there are ergodic theorems for non-stationary (or non-homogeneous) processes that we will not discuss in detail; see, for instance, Nagabhushanam & Bhagavan (1969) and Salehi (1973), and references therein. In practice, we may speak of an ergodic approximation for non-homogeneous process in which the cumulative integral of the ensemble average in a given interval in $x$ is well approximated by the spatial average over the same interval, so that we could have had $\varGamma (x) = R_{x} ({\xi = 0, x}) /\mathscr {E}(x)$ as long as the spatial series is sufficiently long, and thus it converges to the ensemble average. For simplicity, we henceforth use the notation $\langle \zeta ^{2} \rangle$ for the temporal average, while in the appendices we specify whether we speak of temporal or spatial averaging.

Let us now focus on how the wave statistics will adapt to a non-homogeneous correction parameter $\varGamma$ in the ocean. To first order, the probability density of wave heights must fulfil the narrow-band identity

(3.4)\begin{equation} \int_{0}^{+\infty} f(H)\,H^{2} \,{\rm d} H = 8 \langle \zeta^{2} \rangle. \end{equation}

Due to the difficulty of converting surface elevation distributions into crest and height distributions in broad-banded seas (Janssen 2014), (3.4) replaces the typical envelope approach to find the wave height distribution directly. Therefore, neglecting skewness $\mu _{3}$ and kurtosis $\mu _{4}$, the change in the ratio $\langle \zeta ^{2} \rangle /m_{0} \rightarrow \varGamma \times (\langle \zeta ^{2} \rangle /m_{0})$ applied to (2.2) relates the $2\alpha ^{2}$ to $\langle \zeta ^{2} \rangle /m_{0}$, resulting in the narrow-banded correction

(3.5)\begin{equation} \mathcal{R}_{\alpha,\varGamma}(H>\alpha H_{1/3}) = \int_{\alpha}^{+\infty} \frac{4\alpha_{{0}}}{\varGamma}\,{\rm e}^{{-}2\alpha_{{0}}^{2}/\varGamma} \,{\rm d}\alpha_{{0}} = {\rm e}^{{-}2\alpha^{2}/\varGamma}; \end{equation}

that is, the initial wave train Gaussian statistics will be affected by bathymetry, and its exceedance probability $\mathcal {R}_{\alpha }$ will be transformed into $\mathcal {R}_{\alpha, \varGamma }$ by the shoaling process.

Following the same energetics argument as in § 2.1, let us generalise the velocity potential (Dingemans 1997) and surface elevation in (2.5) to an irregular wave field subject to second-order effects in steepness, with uncorrelated random amplitudes and phases (a brief discussion is provided in Appendix A) given by

(3.6)\begin{equation} \left.\begin{gathered} \varPhi (x,z,t) = {\sum_{i}} \sum_{m} {\frac{\varOmega_{m {, i} }(k{_{i}} h)}{mk{_{i}}}} \cosh{(m\varphi)} \sin{(m \phi)},\\ \zeta (x,t) = {\sum_{i}} \sum_{m} \tilde{\varOmega}_{m {, i}} (k{_{i}}h) \cos{(m\phi)}, \end{gathered}\right\} \end{equation}

with the auxiliary variables $\varphi = k_{i} (z+h)$ and $\phi = k_{i}(x-c_{m,i} t+\theta _{i})$, where $c_{m, i}=c_{m}(k_{i})$ is the phase velocity of the $i$th spectral component of $m$th order in steepness. To allow an analytical treatment, the above expression for surface elevation contains no directional effects or wave–wave interaction, and is restricted to the super-harmonic term of the second-order correction. This approximation is supported by Dingemans (1997) and Forristall (2000). In fact, sub-harmonics are at least one order of magnitude smaller than the super-harmonic (see figure 7 of Li et al. 2021c), and given that the super-harmonic contributes to only a 2–3 % change in the energy correction due to shoaling, one should expect the sub-harmonic term to not be fundamental to our analysis. Due to this approximation, however, the treatment is not equivalent to Stokes waves, whose propagation in deep water having a narrow-banded spectrum leads to the modulational instability (Dysthe et al. 2008; Zakharov & Ostrovsky 2009; Onorato et al. 2013).

Under the framework above, we subtract the fixed depth $h_{0}^{2}/2$ term from the energy and can prove for a depth variation $\partial h(x) / \partial x \ll 1$ (see Appendix B) that

(3.7a,b)\begin{equation} \langle \zeta^{2} \rangle = \frac{1}{2} {\sum_{i}} \sum_{m} \tilde{\varOmega}_{m {, i}}^{2},\quad \mathscr{E} = \frac{1}{4} {\sum_{i}} \sum_{m} \left[ \tilde{\varOmega}_{m{, i}}^{2} + \varOmega_{m{,i}}^{2}\, \frac{ \sinh{(2mkh)}}{2mgk} \right]. \end{equation}

In the limit $i \rightarrow \infty$, we can treat the leading-order coefficients $(\varOmega _{m}, \tilde {\varOmega }_{m})$ as being decomposed into a series of even powers of steepness coupled to factored out trigonometric functions $(\chi _{m}, \tilde {\chi }_{m})$ (see Appendix B), leading to

(3.8)\begin{equation} \left.\begin{gathered} \mathscr{E} = \frac{a^{2}}{4} \left[ 1 + \left( \frac{{\rm \pi} \varepsilon}{4} \right)^{2}\chi_{1} + \cdots \right] + \frac{a^{2}}{4} \left[ 1 + \left( \frac{{\rm \pi} \varepsilon}{4} \right)^{2}\tilde{\chi}_{1} + \cdots \right],\\ \langle \zeta^{2} \rangle = \frac{a^{2}}{2} \left[ 1 + \left( \frac{{\rm \pi} \varepsilon}{4} \right)^{2} \tilde{\chi}_{1} + \left( \frac{{\rm \pi} \varepsilon}{4} \right)^{4} \tilde{\chi}_{2} + \cdots \right]. \end{gathered}\right\} \end{equation}

Hence, due to the expressions in (3.3) we get

(3.9)\begin{equation} {\varGamma \approx \frac{\langle \zeta^{2} \rangle }{ \mathscr{E} } = \frac{\displaystyle 1+ \sum_{p} \left( \frac{{\rm \pi} \varepsilon}{4} \right)^{2p} \tilde{\chi}_{p} }{\displaystyle 1+ \frac{1}{2} \sum_{p} \left( \frac{{\rm \pi} \varepsilon}{4} \right)^{2p} \left( \tilde{\chi}_{p} + \chi_{p} \right)}.} \end{equation}

This expression demonstrates that the effect of energetics is reduced to the coefficients of (3.8). We assume that waves before the shoal propagate on a flat bottom and follow the Gaussian distribution associated with the linear wave theory. Afterwards, due to the bathymetry change, second-order corrections become relevant since a much larger steepness is to be taken into account (see Eagleson 1956). Consequently, out-of-equilibrium dynamics will deform an initially Gaussian distribution due to higher-order effects in steepness. In order to apply the same correction to an arbitrary initial distribution, Appendix C shows that modelling any bathymetry change by a transition from first- to second-order terms is very effective. The $\varGamma$ correction depends only on $\varepsilon$ and $kh$. Additionally, it is independent of slope of the shoaling, provided that the latter is sufficiently steep (Zheng et al. 2020). Thus it is applicable to the Trulsen et al. (2020) experiments that featured a slope of $1/3.8$. Indeed, results by Gramstad et al. (2013) suggest that the rogue wave probability is an order of magnitude more sensitive to a relative change in dimensionless depth $k_{p}h$ than in the slope of the shoal. Recently, Fu et al. (2021) has also demonstrated that for steep-shoaling slopes, the probability of rogue waves remains the same.

3.1. Second-order statistics

Under the validity of the above assumptions, the spatial energy density for a second-order perturbation in the narrow-banded case reads (see (B9))

(3.10a,b)\begin{equation} \mathcal{E} = \frac{1}{2} \rho g a^{2} \left[ 1+ {\left( \frac{{\rm \pi} \varepsilon}{4} \right)^{2} \left( \frac{\tilde{\chi}_{1} + \chi_{1}}{2} \right) } \right],\quad {ka = {\rm \pi}\varepsilon \times \mathfrak{S}_{0},} \end{equation}

with coefficients reading (see figure 1a)

(3.11a,b)\begin{equation} \tilde{\chi}_{1} = { \left[ \frac{\cosh kh \left[ 2 + \cosh(2kh) \right]}{\sinh^{3} kh} \right]^{2}},\quad \chi_{1} = \frac{9\cosh(2kh) }{\sinh^{6} kh}. \end{equation}

This model is valid provided that the Ursell number is $Ur =H \lambda ^{2} / h^{3} = \varepsilon (2{\rm \pi} / kh)^{3} \leqslant {8{\rm \pi} ^{2}/3}$ (Lé Méhaute 1976; Dean & Dalrymple 1984). Hence, for small amplitudes ($\zeta \ll h$),

(3.12)\begin{equation} \mathscr{E} = \frac{1}{2} \sum_{i} a_{i}^{2} + \frac{{\rm \pi}^{2}}{16} (\tilde{\chi}_{1}+\chi_{1}) \sum_{i} \frac{a_{i}^{4}}{\lambda_{i}^{2}} \equiv m_{0}, \end{equation}

whereas the surface total variance reads

(3.13)\begin{equation} \langle \zeta^{2} \rangle = \frac{1}{2} \sum_{i} a_{i}^{2} + \frac{\tilde{\chi}_{1} {\rm \pi}^{2}}{8} \sum_{i} \frac{a_{i}^{4}}{\lambda_{i}^{2}}. \end{equation}

In order to consider the whole ensemble of waves, we employ the wavenumber at the peak of the spectrum, $k_{p}$, in the dimensionless depth $k_{p}h$, and the zero-crossing wavelength ${\bar {\lambda }}$ in the significant steepness $\varepsilon = H_{1/3}/{\bar {\lambda }}$ (Massel 2017). Thus, in the limit of a large number of wave components, the previous expressions read

(3.14a,b)\begin{equation} \frac{2\langle \zeta^{2} \rangle}{a^{2}}= 1 + \left(\frac{{\rm \pi} \varepsilon}{4} \right)^{2} \tilde{\chi}_{1},\quad \frac{2\mathscr{E}}{a^{2}} = 1 + \left( \frac{{\rm \pi} \varepsilon}{4} \right)^{2} \frac{(\tilde{\chi}_{1}+\chi_{1})}{2}. \end{equation}

The narrow-banded correction for a group of waves over a changing bathymetry is

(3.15)\begin{equation} \varGamma = \frac{\langle \zeta^{2} \rangle}{\mathscr{E}} = \frac{32 + 2\tilde{\chi}_{1} {\rm \pi}^{2} \varepsilon^{2}}{32 +(\tilde{\chi}_{1}+\chi_{1}) {\rm \pi}^{2} \varepsilon^{2}} \approx 1 + \left( \frac{{\rm \pi} \varepsilon}{4} \right)^{2} \left(\frac{ \tilde{\chi}_{1} - \chi_{1} }{2}\right). \end{equation}

$\varGamma$ exceeds 1 for all depths, with a maximum of up to 1.13 in intermediate depths ($k_{p}h \sim 0.5\text{--} 1$) and an asymptotic behaviour for deep water (see figure 1b). This shape defines three regimes, as marked in figure 1(b). In shallow water (Regime I), a shoal reducing the depth will reduce $\varGamma$, hence the exceedance probability. Conversely, beyond the maximum of $\varGamma$ (Regime II), the shoal will increase $\varGamma$ and the exceedance probability. Finally, as long as the shoal stays within Regime III, the depth variation will translate into a negligible change in $\varGamma$, hence will have no consequence for the exceedance probability. Such absence in amplification is similar to the second-order height distribution in Tayfun (1980) in deep water (see § 3.5). This behaviour will allow us to simplify investigations of shoals starting in deep water ($k_{p}h \gtrsim 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.

Figure 1. (a) Trigonometric coefficients $(\chi _{m}, \tilde {\chi }_{m})$ of the second-order model. (b) Correction parameter $\varGamma$ as a function of steepness $\varepsilon = {H_{1/3}/{\bar {\lambda }}}$ and dimensionless depth $k_{p}h$ in both narrow-banded ($\mathfrak {S}_{0} = 1.0$) and broad-banded ($\mathfrak {S}_{0} = 1.2$) seas.

In order to generalise the derivation of (3.15) to broad-band seas, we use the definition of asymmetry from (2.3). Consequently, the steepness in (3.10a,b)–(3.15) will be corrected by the vertical asymmetry $\mathfrak {S}_{0}$, which in turn modifies the correction parameter:

(3.16)\begin{equation} \varGamma_{\mathfrak{S}{_0}} (\varepsilon, k_{p}h) = \frac{32 + 2\tilde{\chi}_{1} \mathfrak{S}_{0}^{2} {\rm \pi}^{2} \varepsilon^{2}}{32 +(\tilde{\chi}_{1}+\chi_{1}) \mathfrak{S}_{0}^{2} {\rm \pi}^{2} \varepsilon^{2}}. \end{equation}

The vertical asymmetry will increase the correction $\varGamma _{{\mathfrak {S}_0}}$ as compared to the narrow-band case by a few percent, as follows (see figure 1b):

(3.17)\begin{equation} \frac{\varGamma_{\mathfrak{S}{_0}} (\varepsilon, k_{p}h)}{\varGamma (\varepsilon, k_{p}h)} \approx 1 + \frac{1}{2} \left(\tilde{\chi}_{1} - \chi_{1}\right) ( \mathfrak{S}_{0}^{2} - 1 ) \left(\frac{{\rm \pi} \varepsilon}{4} \right)^{2} + {O}(\varepsilon^{4}). \end{equation}

On the other hand, the parameter $\varGamma$ must be corrected for wave breaking, leading to slightly smaller peaks (see figure 2a). We include a depth-dependent breaking limit of regular waves (Miche 1944) by setting $\varepsilon \leqslant (\varepsilon _{0}/7) \tanh {k_{p}h}$ with $0 \leqslant \varepsilon _{0} \leqslant 1$:

(3.18)\begin{equation} \varGamma_{\mathfrak{S}{_0},0} \approx \frac{1600 + 2 {\rm \pi}^{2} \mathfrak{S}_{0}^{2} \varepsilon_{0}^{2} \tilde{\chi}_{1} \tanh^{2}{k_{p}h}}{1600 + {\rm \pi}^{2} \mathfrak{S}_{0}^{2} \varepsilon_{0}^{2} (\tilde{\chi}_{1}+\chi_{1}) \tanh^{2}{k_{p}h}}, \end{equation}

where even if the ratio $\varepsilon _{0}$ is constant, the actual (breaking-limited) steepness $\varepsilon$ will drop considerably throughout the transition from deep to shallow waters (see figure 2a). The vertical asymmetry $\mathfrak {S}_{0}$ increases the value of $\varGamma$ significantly (figure 2b), but less than the decrease in $\varGamma$ due to wave breaking. Although (3.17) typically increases the shoal correction by a few percent, it also shows that the correction has an upper bound $\varGamma \leqslant 1.20$ when $\varepsilon _{0}=1$ and $k_{p}h = 0.5$, and maximum possible asymmetry $\mathfrak {S}_{0}=2$ (skewed sea with $\eta _{1/3} \approx 2.1$). Therefore, a decrease in the depth $k_{p}h$ as in the experiments of Trulsen et al. (2020) will increase the significant steepness as well as $\varGamma$, hence $\mathcal {R}_{\alpha,\varGamma }$ given in (3.5), i.e. the probability of rogue waves as compared to the pre-shoal, homogeneous Rayleigh distribution towards Regime II.

Figure 2. $\varGamma$ correction parameter with the same initial significant steepness in deep water (a) with (dashed) or without (solid) wave breaking in narrow-banded seas, and (b) accounting for wave breaking in narrow-banded (dashed) and broad-banded (solid) seas.

3.2. The normalised height diagram

Since a bandwidth correction to the standard measure $H_{1/3} = 4 \sigma$ will affect the Rayleigh distribution (Longuet-Higgins 1980), such a relation provides a valuable test for the Rayleigh distribution (Goda 1983). Hence the distribution of wave heights contains information on the ratio $H_{1/3}/\sigma$, as follows:

(3.19)\begin{equation} \exp\left(-\frac{H^{2}}{8m_{0}}\right) = \exp\left({- \frac{H^{2}}{8H_{1/3}^{2}} \left(\frac{H_{1/3}}{\sqrt{m_{0}}} \right)^{2}}\right) \equiv \exp\left({- \frac{\alpha^{2}}{8} \left(\frac{H_{1/3}}{\sqrt{m_{0}}} \right)^{2}}\right), \end{equation}

which shows that at fixed energy (fixed $m_{0} = \sigma ^{2}$), the change in the ratio $H_{1/3}/\sigma$ is balanced by the normalisation $\alpha = H/H_{1/3}$, therefore the overall exponent stays invariant, thus the wave statistics. However, if and only if $m_{0}$ increases in comparison to Airy's solution, we expect the coefficient $H_{1/3}^{2}/m_{0}$ attached to $\alpha ^{2}$ to be smaller. In fact, this can be seen by the ratio of the energies at first order (2.6) and second order (3.10a,b). To investigate the effect of broad-band seas on $H_{1/3}/\sigma$ by the direct measurement of the vertical asymmetry as in (2.3), we combine (3.5) and (3.19), and find (see figure 3a)

(3.20)\begin{equation} H_{1/3} = \frac{4}{\mathfrak{S}_{0}} \sqrt{\frac{m_{0}}{\varGamma_{\mathfrak{S}_{0}}}}. \end{equation}

Consequently, the ratios ${H_{1/3}}/{\sqrt {m_{0}}}$ obtained from our model agree with the asymptotic values of 4 (narrow-banded in deep water, solid lines) and 3.8 (broad-banded in deep water, dashed line), as reported by Goda (1983). Asymmetry contributes to this expression via two different processes. On one hand, the $1/\mathfrak {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 depends ultimately on the depth. Increasing the bandwidth and therefore $\mathfrak {S}_{0}$ will also lower $H_{1/3}/\sigma$ (see (3.20)), as demonstrated by Vandever et al. (2008). On the other hand, the asymmetry influences the probability amplification by the shoal, via its impact on $\varGamma _{\mathfrak {S}_{0}}$. Except in very shallow water (Regime I), both $\mathfrak {S}_{0}$ and $\varGamma _{\mathfrak {S}_{0}}$ decrease monotonically as a function of $k_{p}h$, as we model the pre-shoal zone as a homogeneous Airy solution concurrently leading to $\mathfrak {S}_{0} \rightarrow 1$ and $\varGamma _{\mathfrak {S}_{0}} \rightarrow 1$. We can therefore define a parameterisation mapping the latter to the former, in the form $\mathfrak {S}_{0} = \varGamma _{\mathfrak {S}_{0}}^\kappa$ (figure 3b), or, equivalently, its linear version $\mathfrak {S}_{0} = 1 + \kappa (\varGamma _{\mathfrak {S}_{0}}-1)$. Note that this parameterisation aims to describe the simultaneous evolution of $\varGamma$ and the asymmetry when the depth evolves, while (3.17) describes only the direct impact of the asymmetry on $\varGamma$. Since Regime II is the one where we seek the evolution of the exceedance probability, the $\kappa$ parameterisation appears suitable for our work. Furthermore, its use avoids the numerical issues that would arise if handling $\mathfrak {S}_0$ and $\varGamma$ as independent.

Figure 3. (a) $H$–$\sigma$ diagram for narrow-banded $\mathfrak {S}_{0}=\sqrt {\varGamma }$ (solid) and otherwise with $\mathfrak {S}_{0} \sim \varGamma ^{2}$ (dashed) for the specific case of $\varepsilon = 1/15$ for comparison. (b) Ratio $\kappa = \ln {\mathfrak {S}_{0}}/\ln {\varGamma }$ with steepness: $1/7$ (orange), $1/10$ (blue), $1/15$ (green) and $1/20$ (red).

3.3. Evolution of a narrow-banded arbitrary probability distribution

So far, we have studied a pre-shoal Rayleigh distribution of the surface elevation. Following (3.5), the $\varGamma$ correction due to a rapid bathymetry change reads

(3.21)$$\begin{gather} \mathcal{R}_{\alpha,\varGamma} = {\rm e}^{{-}2\alpha^{2}/\varGamma} = \left(\mathcal{R}_{\alpha}\right)^{1/\varGamma}\quad \therefore \quad \frac{\ln \mathcal{R}_{\alpha}}{\varGamma \ln \mathcal{R}_{\alpha,\varGamma}} = 1, \end{gather}$$

(3.22)$$\begin{gather}\frac{\mathcal{R}_{\alpha,\varGamma}}{\mathcal{R}_{\alpha}} = \left( \mathcal{R}_{\alpha} \right)^{({1}/{\varGamma}-1)} = \exp\left({2\alpha^{2} \left(1-\frac{1}{\varGamma}\right)}\right). \end{gather}$$

The same correction therefore applies to any initial surface elevation distribution $\mathbb {P}$, as detailed in Appendix C. For a narrow-banded sea,

(3.23)\begin{equation} {\mathbb{P}_{\alpha,\varGamma} \approx \left( \mathbb{P}_{\alpha} \right)^{{1}/{\varGamma}}\quad \therefore \quad \frac{\mathbb{P}_{\alpha,\varGamma}}{\mathbb{P}_{\alpha}} \approx \exp\left({2\alpha^{2} \left( 1 - \frac{1}{\varGamma}\right)}\right),} \end{equation}

where $\mathbb {P}_{\alpha }$ denotes the exceedance probability at equilibrium prior to the shoal, and $\mathbb {P}_{\alpha,\varGamma }$ is the distribution within the non-equilibrium zone. Hence the highest possible realistic values of $\varGamma \sim 1.08$ (see figure 1) lead to a two-fold increase of rogue wave probability.

3.4. Evolution of a broad-banded arbitrary probability distribution

In the case of a broad-banded sea, the connection between crest and height statistics responsible for (3.20) can be reinterpreted as a simple change of variables $\varGamma \rightarrow \mathfrak {S}_{0}^{2} \varGamma _{\mathfrak {S}_{0}}$ in (3.5), so that (3.23) becomes

(3.24)\begin{equation} \mathbb{P}_{\alpha,\varGamma_{\mathfrak{S}_{0}}} \approx \left(\mathbb{P}_{\alpha} \right)^{{1}/{\mathfrak{S}_{0}^{2} \varGamma_{\mathfrak{S}_{0}}}}\quad \therefore \quad \frac{\mathbb{P}_{\alpha,\varGamma_{\mathfrak{S}_{0}}}}{\mathbb{P}_{\alpha}} \approx \exp\left({2\alpha^{2} \left( 1 - \frac{1}{\mathfrak{S}_{0}^{2} \varGamma_{\mathfrak{S}_{0}}}\right)}\right). \end{equation}

To ensure the numerical stability of the distribution along the propagation, we use the $\kappa$ parameterisation introduced in § 3.2. For practical purposes, we estimate it as

(3.25)\begin{equation} \kappa_{0} = \frac{ \ln{ \left[ {\textrm{max}} \, \mathfrak{S}_{0} \right]} }{ \ln{\left[ \mathrm{max}\, \varGamma \left( \langle \varepsilon \rangle, k_{p}h, \mathfrak{S}_{0} \right) \right] } } \approx \frac{ \ln{ \mathfrak{S}_{0} } }{ \ln{ \varGamma \left( \langle \varepsilon \rangle, k_{p}h, \mathfrak{S}_{0} \right)}}, \end{equation}

where the maximum of $\varGamma$ is taken over the propagation on the whole shoal, $\langle \varepsilon \rangle$ being the average of the pre- and post-shoal steepness. It can be measured at both stages of the wave propagation, or, alternatively, obtained from forecast or hindcast, from the zero-crossing period and significant wave height. Note that uncertainties in asymmetry impact marginally the value of $\kappa _{0}$, as respective errors in the numerator and in $\varGamma$ tend to partially compensate for each other due to (3.17). Consequently, the effect of the shoal on the exceedance probability evolves during the propagation along $x$ as

(3.26)\begin{equation} \frac{\ln \mathbb{P}_{\alpha} }{\ln \mathbb{P}_{\alpha, \varGamma} } = \mathfrak{S}_{0}^{2} \times \varGamma = \left[\varGamma \left(\langle \varepsilon \rangle, k_{p}h, \mathfrak{S}_{0} \right)\right]^{2\kappa_{0}} \times \varGamma \left(\varepsilon (x), k_{p}h (x), \mathfrak{S}_{0} \right). \end{equation}

When the evolution spreads into Regimes I and III, the use of $\kappa _{0}$ is necessary to avoid divergence, though a continuous $\kappa$ restricted to Regime II is applicable. Note that under adequate conditions, a more compact form can be derived, as detailed in Appendix D.

3.5. Comparison with standard second-order models

In this subsection, we delineate similarities and differences between our model of second-order wave height probability and the typical treatment arising from Tayfun (1980). Following the second-order water surface, one finds (Tung & Huang 1985)

(3.27)\begin{equation} \zeta (\phi = 0) := \mathcal{Z}_{c} = \left[ a \cos{\phi} + \frac{ka^{2}}{2} \cos{(2\phi)} \right]_{\phi = 0} = a + \frac{ka^{2}}{2}, \end{equation}

which, normalised by the variance, reduces (with nomenclature $\sigma \tilde {X} := X$) to

(3.28)\begin{equation} \tilde{\mathcal{Z}}_{c} = \tilde{a} + \tilde{a}^{2} \times \frac{k\sigma}{2} \quad \therefore \quad \tilde{a} = \frac{\sqrt{1+2\tilde{\mathcal{Z}}_{c}\sigma k }-1}{\sigma k}. \end{equation}

Applying the finite-depth coefficients in (2.382a) of Dingemans (1997), one finds

(3.29a,b)\begin{equation} \tilde{\mathcal{Z}}_{c} = \tilde{a} + \tilde{a}^{2} \times \frac{k\sigma}{2} \times \mathcal{F},\quad \mathcal{F} =\mathcal{F}(kh) = \frac{3-\tanh^{2}{kh}}{2\tanh^{3}{kh}}, \end{equation}

leading to the mathematical structure (see Mendes et al. (2021) for the $\mathcal {F}=1$ case)

(3.30)\begin{equation} \mathbb{P}(H > \alpha H_{1/3}) = \exp \left\{ - \frac{8 }{k^{2}H_{1/3}^{2}\mathcal{F}^{2}} \left[ \sqrt{1+\alpha kH_{1/3}\mathcal{F}}-1 \right]^{2} \right\}. \end{equation}

Note, however, that this distribution was not derived in the narrow-band model of Tayfun (1980); it is rather an adaptation to extract the mathematical structure of a probability growing boundlessly with increasing steepness. Although Tayfun (2006) puts limitations on the steepness as $\mu _{3} \lesssim kH_{1/3} / (1+\nu ^{2})$, North Sea data (Stansell 2004; Mendes et al. 2021) shows that the skewness can grow higher. Therefore, the suggested adaptation to wave heights is useful since now $\varepsilon$ can reach its breaking limit (Miche 1944). In figure 4(a), we compare our $\varGamma$ model of (3.24) with that of (3.30) through the $H$–$\sigma$ diagram:

(3.31)\begin{equation} \frac{H_{1/3}}{\sqrt{m_{0}}} = \frac{8}{\alpha kH_{1/3} \mathcal{F}} \left[ \sqrt{1+ \alpha kH_{1/3}\mathcal{F} } - 1 \right]. \end{equation}

The finite-depth adapted height model yields 10–20 % lower ratios in deep water, and drops towards zero in shallow water, in contradiction with observations (Goda 1983). Nonetheless, our model can be consistent with the structure of (3.30) in both deep and transitional waters ($k_{p}h \geqslant 0.8$) for moderate values of $\alpha$ and within the narrow-band validity of Tayfun (1980), provided that (3.30) is corrected to the steepness $k_{p}H_{1/3} \rightarrow \varepsilon _{s}$, where $\varepsilon _{s} := \langle \varepsilon _{i} \rangle _{r} = \langle H_{1/3}/\lambda _{i} \rangle _{r}$; see figure 4(b). Indeed, we can use $\langle \lambda _{i} \rangle _{r} \approx 1.5 {\bar {\lambda }}$ (Mendes et al. 2021), which implies $\varepsilon = H_{1/3}/{\bar {\lambda }} \approx 1.5 \varepsilon _{s} \approx 3kH_{1/3}/4{\rm \pi}$. For transitional and especially shallower regions reaching the limit of the second-order approximation, the above adapted model would not reproduce the same probability amplification as our model, nor the original one (Tayfun 1980). Hence our model departs from Tayfun (1980) in finite depth and recovers it in deep water.

Figure 4. (a) $\varGamma$ model (solid) as compared to the finite-depth adjusted second-order structure (dashed), and (b) with corrected steepness.

4. Comparison with Trulsen experiments

Raustøl (2014) provides experiments of wave propagation over a shoal, later summarised in Trulsen et al. (2020). These experiments were carried out in a 24.6 m long and 0.5 m wide wave tank. Surface elevation measurements were made with ultrasound probes, and velocity measurements were made with an acoustic doppler velocimeter, with an array of 16 probes that was moved to four different locations such that the resolution before the shoal was 0.3 m and the resolution above the shoal was 0.1 m. The wavemaker generated a JONSWAP spectrum with peak parameter $\gamma = 3.3$, with typical peak periods $T_{p} \sim 1$ s. The probability distribution of rogue waves evolving with the distance from the wavemaker was recorded, offering a benchmark for our non-homogeneous correction to the wave height probability distribution. Since our model is expressed in terms of dimensionless depth and significant steepness, we have extracted the raw data from figures 5.4 and 5.5 of Raustøl (2014), in accordance with the inversion $k_{p}h \approx (k_{p}a_{c}/ {Ur})^{1/3}$ of Trulsen et al. (2020). 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 E. The experiments feature a shoal starting at $x = 0$, rising up to 42 cm at $x$ = 1.6 m, followed by a plateau until $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. (2020) and table 5.2 of Raustøl (2014).

4.1. Results

Based on the $\kappa _{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 homogeneous Rayleigh distribution (figure 5). The values of $\kappa _{0}$ calculated according to (3.25) 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. Our model that takes into account the evolution of the wave asymmetry $\mathfrak {S}_{0}$ over the shoal (cyan curves) due to the skewness of the surface elevation distribution (2.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. (2020) while keeping a vertical fixed asymmetry $\mathfrak {S}_{0}=1.2$ (blue curves in figure 5) 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 (3.21) leads to a maximal amplification between 75 % and 100 % for $\varGamma \approx 1.08\text{--} 1.10$. When we include the typical vertical asymmetry $\mathfrak {S}_{0}^{2} \sim 1.5$, the pre-shoal probability $\mathbb {P}_{\alpha }$ is transformed into $\mathbb {P}_{\alpha }^{2/3} \sim 10\mathbb {P}_{\alpha }$ within Regime II, seemingly becoming an alternative to Gram–Charlier models (Mori & Yasuda 2002). Finally, adjusting the homogeneous pre-shoal probability to the observed values instead of considering an initial Rayleigh distribution (dotted curve) improves the agreement (runs 6–12, figures 5e–j), demonstrating the applicability of our model to arbitrary probability distributions. This agreement over the whole range of experimental conditions reported by Raustøl (2014) is remarkable, as it requires no specific parameter tweaking.

Figure 5. Evolution of the probability of rogue waves $(\alpha \ge 2)$ over a shoal for Runs 1–2, 4–9, 11–12 from figure 5.8 of Raustøl (2014). Blue: model of (3.23) for a pre-shoal ($x < 0$) Rayleigh distribution with $\mathfrak {S}_{0}=1.2$. Cyan: same model, considering the evolution of skewness over the shoal. Dotted: model of (3.23) for a pre-shoal probability matched to the experimental data. Dashed: same as blue, but for symmetrical seas, $\mathfrak {S}_{0}=1$.

4.2. Discussion

In contrast, Tayfun (1980) (see the discussion following (36)–(38) of that paper) claims a Gaussian probability distribution for the second-order wave heights, as also discussed in Tayfun (1990) and Tayfun & Fedele (2007). In fact, these models fall under the broad category of quasi-determinism theories (Boccotti 2000), in which Longuet-Higgins (1980) and Naess (1985) 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 homogeneous conditions would not be able to describe the Raustøl (2014) and Trulsen et al. (2020) experiments. Also according to these models, any departure from a Rayleigh distribution of wave heights is due to third-order nonlinearities (Tayfun & Fedele 2007; Alkhalidi & Tayfun 2013). Nevertheless, it is important to remark that most of these theories were devised for deep waters, hence they fall within Regime III of our model. Moreover, following Marthinsen (1992) and Dingemans (1997), one can tentatively elaborate alternative finite-depth second-order mathematical structures as in § 3.5. Comparing with our model, we observe a numerical equivalence in deep water, but not in transitional and shallow waters where the alternative second-order models display a sharp departure from the observations of Goda (1983). Unfortunately, the finite-depth model of Tayfun & Alkhalidi (2020) is not a distribution of wave heights, hence is not applicable to our discussion.

Our $\varGamma$ model of the probability evolution over a shoal is based on the second-order correction, which is valid for ${Ur} \leqslant 8{\rm \pi} ^{2}/3$ (Dean & Dalrymple 1984). This limit of validity is well beyond the Ursell number ${Ur} \leqslant 0.22$ of the experiments considered in this work (Raustøl 2014; Trulsen et al. 2020). We can therefore expect that the $\varGamma$ model will still apply to more than 10 times larger waves and/or shallower waters.

The transient drop of $\mathbb {P}_{\alpha, \varGamma }/\mathbb {P}_{\alpha }$ to almost zero in the rising region of the shoal ($0 \leqslant x \le 1.6$ m) in runs 1, 2, 4 and 5, as well as its slow decay on the trailing side ($3.2 \le x \le 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 & Benoit (2021) described the exceedance probability of Run 3 of Trulsen et al. (2020) throughout the shoaling and de-shoaling stages by means of numerical simulations. However, an analytical 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 provides an explicit expression with no free fitting parameter. It is obtained directly from second-order correction and has been shown to be valid for all runs in Raustøl (2014). Although we were not able to analyse Run 3 (Jorde 2018) 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.

The modulational instability cannot account for the rise of the rogue wave probability due to shoaling in the considered experiments because they do not feature narrow-banded Stokes waves or lie in the optimal range $k_{p}h > 1.36$ (Zakharov & Ostrovsky 2009) in the shallower side, except for Runs 11 and 12, which showed no significant amplification. On the other hand, all runs start in the range $k_{p}h > 1.36$ before the shoal and showed no large deviation from the Gaussian sea, except for Run 12.

Furthermore, our model – and in particular the three regimes discussed in (3.15) and figure 1(b) – also allows us to understand the contradictory behaviours highlighted by Trulsen et al. (2020) and Zhang & Benoit (2021). 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 (Glukhovskii 1966), i.e. rogue wave likelihood lower than in transitional waters. In addition, Barbariol et al. (2015) show an increase of maximum crest up to some cut-off in dimensionless depth, upon which it starts to decrease sharply 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 rapid depth change for a dimensionless depth $0.5 < k_{p}h < 1.5$ but stops growing when it leaves this range. It is therefore less pronounced past the threshold $k_{p}h \leqslant 0.3$ until it starts to follow shallow water distributions such as in Glukhovskii (1966). Recently, after the submission of this work, experiments on the shoaling of irregular unidirectional waves have confirmed our prediction that the amplification of wave height statistics in Regime II will vanish in Regime I (Xu et al. 2021).

The evolution of the rogue wave probability highlighted by our model could be related to the process proposed by Li et al. (2021c) and confirmed experimentally in Li et al. (2021a), in which the generation of additional wave packets that interact with the original pre-shoal wave packet propagating over a step leads to a local peak some distance into the shallower region. Notably, Li et al. (2021c) have used a more general treatment for the surface elevation, containing both super- and sub-harmonics, as well as free and bound waves. Moreover, instead of bound super-harmonics in our model, one could investigate the effect of the nonlinear evolution of interacting free modes. For instance, similar experimental results were interpreted using a statistical matching of Korteweg–De Vries equilibrium states at the depth transition point (Majda et al. 2019; Moore et al. 2020). In retrospect to the ideas in Onorato & Suret (2016), these works obtained a connection between variance and skeweness with the dynamics. Therefore, subsequent investigations carrying out a comparison between the two types (free and bound) of mode evolution will be of great relevance.

5. Conclusion

This work presented a connection between statistical distributions and the fluid mechanics of the second-order perturbation in non-homogeneous conditions, providing successfully a physical explanation for the rogue wave probability increase over a depth change (Raustøl 2014; Trulsen et al. 2020). We have shown that our model reproduces very well the experiments of Raustøl (2014) and Trulsen et al. (2020) regarding the probability distribution as a function of the distance from the wavemaker. Moreover, we showed that the significant steepness and dimensionless depth affect the validity, and assessed numerically the extent of the deviation from the assumption of homogeneity. Furthermore, instead of introducing new physics (Haver & Andersen 2000), our model has demonstrated that an effective theory arises by challenging the homogeneity assumption. Introducing the $H_{1/3}/\sigma$ ratio (Goda 1983), we have established that the deep water regime produces no 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 decreases this large amplification to a level smaller than the initial stage in deep water.

Our model has been restricted to reformulating normalised moments up to the variance only. The evolution of either skewness or kurtosis as a function of the distance from the wavemaker (Trulsen et al. 2020) will 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 or kurtosis, as they are ‘symptoms’ of the dynamics and not the cause (Stansell 2004; Christou & Ewans 2014; Cattrell et al. 2018). Since our model relies on relations between steepness, slope and bandwidth that might be affected by a rapid depth change, the empirical findings in Mendes et al. (2021) regarding vertical asymmetry have to be extended to the current setting for an exact formulation. Furthermore, the generalisation of this work to multidirectional spectra would be needed, since Ducrozet & Gouin (2017) suggest that such a 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.