1. Introduction
Quantifying the impact of rogue waves on the stability of offshore and coastal structures has become a recent direction of research. An increased frequency of rogue waves in a time series may lead to an excess in wave loads (Tang et al. Reference Tang, Barratt, Bingham, van den Bremer and Adcock2022; He et al. Reference He, Kanehira, Mori, Gamaleldin, Babanin, Chauhan and Chabchoub2023; Li et al. Reference Li, Tang, Li, Draycott, van den Bremer and Adcock2023; Ma & Swan Reference Ma and Swan2023; Xin, Li & Li Reference Xin, Li and Li2023). The statistics of water waves define the design envelope for ocean vessels, and the possibility of avoiding or mitigating these extreme waves is important for marine and coastal safety. The risk of rogue waves in deep water is nowadays rather well known (see e.g. Pelinovsky & Kharif (Reference Pelinovsky and Kharif2008), Chabchoub, Hoffmann & Akhmediev (Reference Chabchoub, Hoffmann and Akhmediev2011), Chabchoub et al. (Reference Chabchoub, Hoffmann, Onorato and Akhmediev2012), Fedele et al. (Reference Fedele, Brennan, Ponce de León, Dudley and Dias2016) and references therein), but rogue waves are also observed in intermediate (Trulsen, Zeng & Gramstad Reference Trulsen, Zeng and Gramstad2012; Zeng & Trulsen Reference Zeng and Trulsen2012; Gramstad et al. Reference Gramstad, Zeng, Trulsen and Pedersen2013; Viotti & Dias Reference Viotti and Dias2014) and shallow waters (Sergeeva, Pelinovsky & Talipova Reference Sergeeva, Pelinovsky and Talipova2011; Didenkulova & Pelinovsky Reference Didenkulova and Pelinovsky2016; Mendes & Scotti Reference Mendes and Scotti2021; Karmpadakis, Swan & Christou Reference Karmpadakis, Swan and Christou2022; Didenkulova, Didenkulova & Medvedev Reference Didenkulova, Didenkulova and Medvedev2023) subject to bathymetry gradients. Among the mechanisms associated with rogue wave formation and amplification (Dysthe, Krogstad & Muller Reference Dysthe, Krogstad and Muller2008; Onorato et al. Reference Onorato, Residori, Bortolozzo, Montina and Arecchi2013; Adcock & Taylor Reference Adcock and Taylor2014), shoaling leads to the highest recorded excess in kurtosis (Li & Chabchoub Reference Li and Chabchoub2023).
When propagating nearshore, the hydrodynamics of wind-generated surface waves is affected by the relative water depth until wave breaking becomes dominant (Green Reference Green1838; Burnside Reference Burnside1915; Holthuijsen Reference Holthuijsen2007), and modifies the wave dispersion, wavelength, group velocity, wave height, among others. Several models describe the deformation and transformation of waves encountering bathymetric changes (Eagleson Reference Eagleson1956; Shuto Reference Shuto1974; Walker & Headlam Reference Walker and Headlam1982; Kweon & Goda Reference Kweon and Goda1996), but the rogue wave amplification by the non-equilibrium process of shoaling cannot be calculated directly from the shoaling coefficient and the behaviour of regular waves only. Connecting the knowledge of wave transformation with the implications thereof to anomalous wave statistics, theoretical, experimental and numerical studies investigated the possibility of the amplification of rogue waves subject to shoaling in the recent years (Kashima, Hirayama & Mori Reference Kashima, Hirayama and Mori2014; Ma, Ma & Dong Reference Ma, Ma and Dong2015; Ducrozet & Gouin Reference Ducrozet and Gouin2017; Majda, Moore & Qi Reference Majda, Moore and Qi2019; Zhang et al. Reference Zhang, Benoit, Kimmoun, Chabchoub and Hsu2019; Moore et al. Reference Moore, Bolles, Majda and Qi2020; Trulsen et al. Reference Trulsen, Raustøl, Jorde and Rye2020; Zheng et al. Reference Zheng, Lin, Li, Adcock, Li and van den Bremer2020; Doeleman Reference Doeleman2021; Fu et al. Reference Fu, Ma, Dong and Perlin2021; Gomel et al. Reference Gomel, Chabchoub, Brunetti, Trillo, Kasparian and Armaroli2021; Kimmoun et al. Reference Kimmoun, Hsu, Hoffmann and Chabchoub2021; Li et al. Reference Li, Draycott, Adcock and van den Bremer2021a,Reference Li, Zheng, Lin, Adcock and van den Bremerc; Lyu, Mori & Kashima Reference Lyu, Mori and Kashima2021; Xu et al. Reference Xu, Liu, Li and Jia2021; Zhang & Benoit Reference Zhang and Benoit2021; Lawrence, Trulsen & Gramstad Reference Lawrence, Trulsen and Gramstad2022; Mendes et al. Reference Mendes, Scotti, Brunetti and Kasparian2022; Lyu, Mori & Kashima Reference Lyu, Mori and Kashima2023; Zhang et al. Reference Zhang, Ma, Tan, Dong and Benoit2023). A key component of this amplification is the abruptness of the environmental transition, which often drives physical systems out of equilibrium (Lockwood Reference Lockwood2001; Sobolev Reference Sobolev2013; Steinbach, Gemming & Erbe Reference Steinbach, Gemming and Erbe2016; Passiatore et al. Reference Passiatore, Sciacovelli, Cinnella and Pascazio2022), thereby generating transient local phenomena until the new equilibrium state is established (Zhang, Ma & Benoit Reference Zhang, Ma and Benoit2024). Trulsen (Reference Trulsen2018) anticipates that variations in environmental or meteorological conditions over a rapid spatiotemporal scale will in general lead to out-of-equilibrium states and thus to anomalous wave statistics. This type of anomalous statistics is also observed in the generation of irregular fields in wave flumes or basins when the sea-states are forced by random input spectra. Indeed, Tang et al. (Reference Tang, Xu, Barratt, Bingham, Li, Taylor, van den Bremer and Adcock2020) found that the kurtosis behaves anomalously within short spatiotemporal scales for a randomly imposed initial wave spectrum.
While the shoaling process involves several physical parameters and the anomalous wave statistics thereof are strongly dependent on the wave steepness and relative water depth (Li et al. Reference Li, Draycott, Zheng, Lin, Adcock and van den Bremer2021b; Zhang & Benoit Reference Zhang and Benoit2021; Mendes et al. Reference Mendes, Scotti, Brunetti and Kasparian2022), the geometry of a linear shoal is simply characterized by three parameters: the slope; its spatial extension; the change in depth. These parameters, however, provide only two degrees of freedom. Here we consider the slope magnitude and the slope length (equivalently called ‘shoaling length’ hereafter), from which the depth change can be determined correspondingly. Although Zheng et al. (Reference Zheng, Lin, Li, Adcock, Li and van den Bremer2020) and Lawrence, Trulsen & Gramstad (Reference Lawrence, Trulsen and Gramstad2021) examined the effect of shoaling length on wave statistics and found that both skewness and kurtosis decreased as this length increased, such studies kept the depth difference constant, so that the effects of the bottom slope and shoaling length remain entangled. To date, an analytical expression of the effect of the shoaling length is missing. Mendes & Kasparian (Reference Mendes and Kasparian2022) provided an expression for the effect of the bottom slope, but they fixed the shoaling length equal to one wavelength. It is often argued that disturbances in wave statistics occur due to the ‘abrupt’ nature of the depth transition, which drives the system out of equilibrium. Especially for a submerged step, when wave packets evolving slowly in time and space meet an abrupt depth transition, a standing wave pattern between the superharmonic components of the carrier waves (free and bound) is modulated by their group structure (Li et al. Reference Li, Zheng, Lin, Adcock and van den Bremer2021c). However, the term ‘abruptness’ is generally used in an intuitive manner, without clearly defining whether it refers to a strong bathymetry change, a spatial extension shorter than the wavelength or a steep slope (Viotti & Dias Reference Viotti and Dias2014; Zheng et al. Reference Zheng, Lin, Li, Adcock, Li and van den Bremer2020; Lawrence et al. Reference Lawrence, Trulsen and Gramstad2021; Li et al. Reference Li, Draycott, Zheng, Lin, Adcock and van den Bremer2021b; Draycott et al. Reference Draycott, Li, Stansby, Adcock and van den Bremer2022). In an engineering context, the term ‘abrupt depth transition’ is also used for a variety of natural and man-made bathymetric features (Draycott et al. Reference Draycott, Li, Stansby, Adcock and van den Bremer2022) although their slopes are well below the threshold of 1/3 defined for abruptness by the mild slope equation (Booij Reference Booij1983). In the latter context, abruptness rather appears to refer to the depth difference. These ambiguities result in confusing and sometimes seemingly contradictory results.
There is widespread agreement that the wave steepness growth up to the second order (and sometimes to the third order) is the main driver of out-of-equilibrium phenomena and anomalous wave statistics in intermediate water depths. However, experimentations targeting the isolation of the effects of slope magnitude and shoaling length are still lacking. Here, we define a dimensionless shoaling length parameter 
$\ell$ as the ratio of the length of the shoal to the characteristic wavelength. We disentangle the effect on rogue wave probability of this measure of the wave shoaling ‘abruptness’ from the effect of slope magnitude. Relying on a non-homogeneous second-order theory and a set of numerical simulations using a fully nonlinear potential flow (FNPF) model, we vary the length of the shoal while keeping a fixed slope magnitude and the same sea-state conditions atop the shoal. We show that the effect of a shoal on the rogue wave probability amplification is not governed by the shoaling length, as soon as this length exceeds approximately half of the peak wavelength, thus leaving the slope magnitude as the dominant factor.
The remainder of this article is organized as follows. In § 2 an existing non-homogeneous second-order theory is further extended to provide a simple explicit approximation of the maximum value of kurtosis atop the shoal, while § 3 presents the FNPF model used for the present high-fidelity numerical simulations. The results from the two approaches are compared and discussed in § 4. Conclusions are summarized in § 5.
2. Non-homogeneous second-order wave theory
The present section investigates the effect of the shoaling length by disentangling it from the slope magnitude 
$|\boldsymbol {\nabla }h|$, with 
$h$ denoting the still-water depth (
$h > 0$) and 
$\boldsymbol {\nabla }$ the gradient operator in the horizontal plane. We define the shoaling length as 
$\ell \equiv L/\lambda$, with 
$L$ denoting the length of the shoal and 
$\lambda$ the dimensional wavelength, taken as the zero-crossing wavelength. We first recall the key outputs of Mendes et al. (Reference Mendes, Scotti, Brunetti and Kasparian2022) in § 2.1, then apply them to a compact expression for the kurtosis enhancement in § 2.2.
2.1. Background of non-homogeneous second-order wave theory
We rely on a non-homogeneous stochastic theory describing the energy redistribution among modes of irregular water waves travelling over a shoal of arbitrary slope magnitude (Mendes & Kasparian Reference Mendes and Kasparian2023), see figure 1 for a graphical description of the problem. Taking into account the disturbance due to the shoaling on the wave energetics, we may approximate the exceedance probability 
$\mathcal {R}$ (for wave height 
$H$ exceeding 
$\alpha$ times the significant wave height 
$H_s$) in a narrow-banded irregular wave train moving past a steep slope through the formula (Mendes et al. Reference Mendes, Scotti, Brunetti and Kasparian2022)
\begin{equation} \mathcal{R}(H>\alpha H_{s}) = {\rm e}^{{-}2\alpha^{2}/\varGamma}, \end{equation}
where the non-homogeneous parameter 
$\varGamma$ is defined as
\begin{equation} \varGamma \equiv \frac{ \mathbb{E}[\eta^{2}] }{ \mathscr{E} } , \end{equation}
with 
$\eta$ denoting the free surface elevation (FSE) and 
$\mathbb {E}[\eta ^{2}]$ denoting the ensemble average of 
$\eta ^2$. The ensemble average is computed through
\begin{equation} \mathbb{E}\left[ \eta^2 \right] = \int_{-\infty}^{+\infty} \eta^2 p(\eta) \, \textrm{d}\eta , \end{equation}
with 
$p(\eta )$ denoting the probability distribution function of 
$\eta$. In practice, 
$\mathbb {E}[\eta ^{2}]$ is approximated by the variance 
$\langle \eta ^2 \rangle$ of 
$\eta$ in a discrete time series. Here 
$\mathscr {E}$ is defined as the total mechanical energy of the waves averaged over one wavelength (see Dean & Dalrymple (Reference Dean and Dalrymple1991), for instance) and normalized by 
$\rho g$, where 
$\rho$ is the water density and 
$g$ the acceleration due to gravity. In the footsteps of Dean & Dalrymple (Reference Dean and Dalrymple1991), and considering a two-dimensional Cartesian coordinate system 
$(x,z)$ with its 
$z$-axis origin located at the still-water level, we calculate the energy 
$\mathscr {E}$ averaged over one zero-crossing wavelength (Dong et al. Reference Dong, Gao, Ma and Ma2020; Mendes et al. Reference Mendes, Scotti, Brunetti and Kasparian2022) as
\begin{equation} \mathscr{E}(x) = \frac{1}{2\lambda} \int_{x}^{x+\lambda} \eta^{2} \, \textrm{d}\kern0.05em x + \frac{1}{2g\lambda } \int_{x}^{x+\lambda} \int_{{-}h(x)}^{0} \left[ \left( \frac{\partial \varPhi}{\partial x} \right)^{2} + \left( \frac{\partial \varPhi}{\partial z} \right)^{2} \right] \textrm{d}z \, \textrm{d}\kern0.07em x , \end{equation}
where 
$\varPhi$ is the velocity potential. For linear regular waves, we know that 
$\mathbb {E}[\eta ^{2}]= \mathscr {E}=a^{2}/2$ (Airy Reference Airy1845), and thus 
$\varGamma =1$. Traditionally, the wave energy integration in (2.4) is implemented over space (Le Méhauté Reference Le Méhauté1976; Dean & Dalrymple Reference Dean and Dalrymple1991; Fredsøe & Deigaard Reference Fredsøe and Deigaard1992) for the convenience of studying spatial change of water depth, although time integration could also be used (Holthuijsen Reference Holthuijsen2007).
Sketch (not to scale) of the geometry of the shoaling problem, with incident waves coming from the left.
In the limit of a large number of wave components (Mendes et al. Reference Mendes, Scotti, Brunetti and Kasparian2022), for the purpose of the stochastic wave analysis one may use a simplified monochromatic velocity potential up to second order in steepness instead of the full second-order treatment as in Sharma & Dean (Reference Sharma and Dean1981). Consequently, without loss of generality, we start with a second-order Stokes wave solution. The velocity potential reads (Dingemans Reference Dingemans1997)
\begin{equation} \varPhi (x,z,t) = \frac{a\omega }{k} \left[ \frac{\cosh{\varphi} }{\sinh{ \mu }} \sin{\phi} + \left( \frac{3 ka}{8} \right) \frac{\cosh{(2\varphi)} }{\sinh^{4}{\mu}} \sin{(2\phi)} \right] , \end{equation}
where 
$\varphi \equiv k (z+h)$, 
$\phi \equiv kx - \omega t$, and 
$a$, 
$k$, 
$\omega$ denote the amplitude, wavenumber and angular frequency of the regular wave, respectively. Here 
$\mu \equiv kh$ denotes the relative water depth, and 
$ka$ is the wave steepness. The FSE solution of Stokes’ second-order theory reads (Dingemans Reference Dingemans1997)
\begin{equation} \eta (x,t) = a \left[ \cos{\phi} + \left( \frac{ ka }{ 4 } \right)\sqrt{\tilde{\chi}_{1}} \cos{(2\phi)} \right] , \end{equation}where the superharmonic term takes the form
\begin{equation} \tilde{\chi}_{1} = \frac{(3 - \sigma^2)^2} {\sigma^6} \quad \textrm{and} \quad \sigma \equiv \tanh\mu . \end{equation}
Furthermore, the horizontal and vertical velocity components 
$(u,w)$ are formulated as
$$\begin{gather} u \equiv \frac{\partial \varPhi}{\partial x} = a\omega \left[ \frac{ \cosh{\varphi} }{ \sinh{ \mu } } \cos{\phi} + \left( \frac{3 ka}{4} \right) \frac{\cosh{(2\varphi)} }{\sinh^{4}{ \mu }} \cos{(2\phi)} \right] , \end{gather}$$
$$\begin{gather}w \equiv \frac{\partial \varPhi}{\partial z} = a\omega \left[ \frac{ \sinh{\varphi} }{ \sinh{ \mu } } \sin{\phi} + \left( \frac{3ka }{4} \right) \frac{ \sinh{(2\varphi)} }{ \sinh^{4}{ \mu } } \sin{(2\phi)} \right] . \end{gather}$$
For an irregular wave field described by a second-order perturbation in steepness 
$\varepsilon = H_s / \lambda$ in a relative water depth 
$\mu = k_p h$, whose transition from a regular field is detailed in Mendes et al. (Reference Mendes, Scotti, Brunetti and Kasparian2022), the computation of (2.2) leads to
\begin{equation} \varGamma = \frac{ 1+ \dfrac{{\rm \pi}^{2} \varepsilon^{2} \mathfrak{S}^2 }{16} \, \tilde{\chi}_{1} }{ 1+ \dfrac{{\rm \pi}^{2} \varepsilon^{2} \mathfrak{S}^2 }{32} \left( \tilde{\chi}_{1} + \chi_{1} \right) +\check{\mathscr{E}}_{p2} } , \end{equation}
with a transformation 
$ka \rightarrow {\rm \pi}\varepsilon \mathfrak {S}$, where 
$\mathfrak {S}$ is the vertical asymmetry between the crests and troughs of the waves (see its definition in Mendes, Scotti & Stansell (Reference Mendes, Scotti and Stansell2021) and Mendes & Kasparian (Reference Mendes and Kasparian2023)). Here 
$\check {\mathscr {E}}_{p2}$ is the non-dimensional net change in potential energy due to the set-down over the shoal. The expression of 
$\tilde {\chi }_{1}$ is given in (2.7a,b), and the term 
$\chi _{1}$ reads
\begin{equation} \chi_{1} = \frac{9\,\textrm{cosh}(2\mu) }{\textrm{sinh}^{6} \mu} =\frac{9}{\sigma^6}(1-\sigma^2)(1-\sigma^4). \end{equation}
In the case of irregular waves, the relation between 
$\varepsilon$ and 
$k_pH_s$ arises from the difference between peak and zero-crossing wave periods, and therefore depends on the spectral peakedness (Ochi Reference Ochi1998). The spectral peakedness hardly changes a few wavelengths after the the end of the shoal (Zhang, Benoit & Ma Reference Zhang, Benoit and Ma2022), and consequently, we may use 
$\varepsilon \equiv (\sqrt {2}/{\rm \pi} ) k_{p}H_{s}$. Alternatively, one could try to integrate over the entire interval 
$[-L,0]$. This would give rise to corrections in the results of Mendes et al. (Reference Mendes, Scotti, Brunetti and Kasparian2022) containing 
$\varepsilon ^{2}$. Such corrections are expected to be quickly damped as 
$\ell > 0.5$ and much smaller than the main result of (2.10).
The effect of the slope magnitude is generally disregarded in investigations of wave fields evolving over a planar beach. Derivations for integral properties (momentum, energy, energy flux) often relied on the assumption that 
$|\boldsymbol {\nabla } h | \lesssim \varepsilon$ (Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1964; Turpin, Benmoussa & Mei Reference Turpin, Benmoussa and Mei1983; Iusim & Stiassnie Reference Iusim and Stiassnie1985; Porter Reference Porter2003; Mei, Stiassnie & Yue Reference Mei, Stiassnie and Yue2005). But for relatively steep slopes, as considered in the present work, Mendes & Kasparian (Reference Mendes and Kasparian2022) have shown that changes in potential energy are significant over a relatively steep slope and may affect wave statistics. The net change of the potential energy over a shoaling zone due to set-down 
$\check {\mathscr {E}}_{p2}$ is a function of the slope magnitude 
$|\boldsymbol {\nabla } h|$. We consider the slope function for deep water preshoal conditions approximated as (Mendes & Kasparian Reference Mendes and Kasparian2022)
\begin{equation} \check{\mathscr{E}}_{p2} \approx \frac{\varepsilon^{2}}{\mu^2} \, f(|\boldsymbol{\nabla} h|), \quad \textrm{with} \ f(|\boldsymbol{\nabla} h|) \equiv \frac{ 1 }{25 |\boldsymbol{\nabla} h|} - 5| \boldsymbol{\nabla} h| \left( 1 - |\boldsymbol{\nabla} h| \right) . \end{equation}
This approximation is valid for the range of slopes 
$1/100 \leqslant |\boldsymbol {\nabla } h| \leqslant 1/2$. Fixing the pair 
$(\mu, \varepsilon )$, at steep slopes (i.e. in the vicinity of 
$|\boldsymbol {\nabla } h| \sim 1/2$) the function 
$\check {\mathscr {E}}_{p2}$ reaches its maximum value.
2.2. Excess kurtosis of shoaling waves over arbitrary slopes
The skewness and kurtosis (denoted as 
$\lambda _3$ and 
$\lambda _4$, respectively) of the FSE are often used to indicate the degree of nonlinearity of the exceedance probability (Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1962). For a random variable 
$X$ with a zero-mean, they are defined as
\begin{equation} \lambda_3(X) = \frac{\mathbb{E}\left[ X^{3} \right]}{\left(\mathbb{E}\left[ X^2 \right]\right)^{3/2}} , \quad \lambda_4(X) = \frac{\mathbb{E}\left[ X^{4} \right]}{\left(\mathbb{E}\left[ X^2 \right]\right)^{2}}, \end{equation}
where 
$\mathbb {E}\ [\,]$ denotes the expected value (or ensemble average). For 
$X$ following a Gaussian distribution, 
$\lambda _3(X)=0$ and 
$\lambda _4(X)=3$. When 
$X$ represents 
$\eta$, skewness serves as an indicator of wave asymmetry in the vertical direction, and kurtosis is positively related to the occurrence frequency of extreme values. For mathematical convenience, we build our theory for the excess kurtosis 
$\hat {\lambda }_4 \equiv \lambda _4 - 3$ (thus 
$\hat {\lambda }_4 =0$ in a Gaussian sea).
To derive the excess kurtosis of FSE from the non-homogeneous second-order wave theory, we take advantage of an effective theory connecting the non-homogeneous parameter 
$\varGamma$ and the excess kurtosis 
$\hat {\lambda }_{4}$, namely (3.3) of Mendes & Kasparian (Reference Mendes and Kasparian2023):
\begin{equation} \hat{\lambda}_{4} \approx \frac{1}{9} \left[ \exp\left({ 8 \left( 1 - \frac{ 1 }{ \mathfrak{S}^{2}\varGamma } \right)} \right) - 1 \right] . \end{equation}
The preshoal value of the excess kurtosis 
$\hat {\lambda }_{4,0}$ (the subscript 0 is used henceforward to denote preshoal conditions) is also relevant to account for bias in initial conditions when varying 
$L$, such that
\begin{align} \Delta \hat{\lambda}_{4} &\equiv \hat{\lambda}_{4} - \hat{\lambda}_{4,0} = \frac{1}{9} \left[ \exp\left( { 8 \left( 1 - \frac{ 1 }{ \mathfrak{S}^{2}\varGamma } \right) }\right) - \exp\left({ 8 \left( 1 - \frac{ 1 }{ \mathfrak{S}_{{0}}^{2}\varGamma_0 } \right) }\right) \right] \nonumber\\ &\equiv(\lambda_4 - 3) - (\lambda_{4,0} - 3) = \Delta \lambda_4 . \end{align}
Having the exact values of the variables 
$(\varepsilon,\mu )$ prior and atop the shoal, we can compute the change of 
$\varGamma$ between the preshoal and postshoal regions and translate it into the change of kurtosis 
$\Delta \lambda _4$.
In order to build a more intuitive view on (2.15), we provide a manageable explicit expression approximating the excess kurtosis, by Taylor-expanding the non-homogeneous correction given in (2.10) in 
$\varepsilon ^2$. Note that the net potential energy is also of order 
$\varepsilon ^2$ in (2.12):
\begin{equation} \varGamma \approx 1 + \left( \frac{{\rm \pi} \varepsilon}{4} \right)^{2} \left( \frac{ \tilde{\chi}_{1} - \chi_{1} }{2}\right) - \check{\mathscr{E}}_{p2} . \end{equation}
As the water depth decreases towards shallow waters, the ratio between the trigonometric coefficients 
$\tilde {\chi }_{1}$ and 
$\chi _{1}$ can be fitted as a polynomial to leading order as
\begin{equation} \frac{\tilde{\chi}_{1}}{\chi_{1}} \approx 1 + 1.534 \mu^{4} \approx 1 + 1.5 \mu^{4} , \quad \text{for} \ \frac{1}{2} < \mu < \frac{3}{2} . \end{equation}
With such a simplification and the expression of 
$\check {\mathscr {E}}_{p2}$ given in (2.12), we can then rewrite the correction 
$\varGamma$ in (2.16) as follows:
\begin{equation} \varGamma \approx 1 + \left( \frac{{\rm \pi} \varepsilon}{4} \right)^{2} \frac{ \chi_{1}}{2} \times 1.5 \mu^{4}- \frac{\varepsilon^{2}}{\mu^{2}} \times f(|\boldsymbol{\nabla} h|) \approx 1 + \frac{\varepsilon^{2}}{4\mu^{2}} \left[ 2{\rm \pi}^{2} - 4f(|\boldsymbol{\nabla} h|) \right]. \end{equation}
As described in § 3.4 of Mendes et al. (Reference Mendes, Scotti, Brunetti and Kasparian2022), the effect of the vertical asymmetry of the wave profile can be parameterized for second-order waves at the region of the peak in amplification for 
$\mu \sim 1/2$:
\begin{equation} \mathfrak{S}^{2}\varGamma \sim \varGamma^{6} . \end{equation}
Hence, as soon as the steepness is sufficiently large (i.e. 
$\varepsilon \geqslant 1/50$), the argument of the exponential function in the last part of (2.14) can be expanded by using (2.18) and (2.19):
\begin{align} 1 - \frac{ 1 }{ \mathfrak{S}^{2}\varGamma } & \approx 1 - \left\{ 1 + \frac{\varepsilon^{2}}{4\mu^{2}} \left[ 2 {\rm \pi}^{2} - 4f(|\boldsymbol{\nabla} h|) \right] \right\}^{{-}6}\nonumber\\ &\approx 1 - \left\{ 1 - \frac{3\varepsilon^{2}}{2\mu^{2}} \left[ 2 {\rm \pi}^{2} - 4f(|\boldsymbol{\nabla} h|) \right] \right\}\nonumber\\ & \approx \frac{3\varepsilon^{2}}{2\mu^{2}} \left[ 2 {\rm \pi}^{2} - 4f(|\boldsymbol{\nabla} h|) \right] = \frac{3{\rm \pi}^2 \varepsilon^{2}}{\mu^{2}} -6 \check{\mathscr{E}}_{p2} . \end{align}
When the slope 
$|\boldsymbol {\nabla } h|$ is sufficiently steep, (2.20) saturates, so that the parameter 
$\varGamma$ becomes independent of the slope. We plug this expression into (2.14), finding
\begin{equation} \hat{\lambda}_{4} = \tfrac{1}{9} \left( \exp({ 24 {\rm \pi}^{2} \varepsilon^{2} /\mu^2 }) \exp({-48 \check{\mathscr{E}}_{p2}})- 1 \right) . \end{equation}
We can manipulate this algebraically by adding and subtracting 
$\exp ({ -48 \check {\mathscr {E}}_{p2} })$, which leads to
\begin{align} \hat{\lambda}_{4} &= \tfrac{1}{9} \left[{ \exp({ -48 \check{\mathscr{E}}_{p2} }) } \left( \exp({ 24 {\rm \pi}^{2} \varepsilon^{2} /\mu^2 })- 1 \right) + \left( \exp({ -48 \check{\mathscr{E}}_{p2}}) -1 \right) \right] ,\nonumber\\ & \equiv \exp({ -48 \check{\mathscr{E}}_{p2} }) \times \hat{\lambda}_{4 , b} + \tfrac{1}{9} \left( \exp({ -48 \check{\mathscr{E}}_{p2}}) -1 \right) . \end{align}
The second term is one order of magnitude smaller than the slope-independent kurtosis 
$\hat {\lambda }_{4 , b}$ since 
$|\check {\mathscr {E}}_{p2}| \sim \varepsilon ^2 \sim 10^{-3}$ from (2.12), thus 
$24 {\rm \pi}^{2} \varepsilon ^{2} /\mu ^2 \gg -48 \check {\mathscr {E}}_{p2}$. As a consequence,
\begin{equation} \hat{\lambda}_{4} \approx \frac{\exp({ -48 \check{\mathscr{E}}_{p2} })}{9} \left( \exp({ 24 {\rm \pi}^{2} \varepsilon^{2} /\mu^2 })- 1 \right) . \end{equation}Furthermore, the net potential energy form given by (2.12) can be approximated as follows (Mendes Reference Mendes2024):
\begin{equation} \check{\mathscr{E}}_{p2} \approx \frac{\varepsilon^{2}}{4\mu^2} \left( 7 - 20 \sqrt{ |\boldsymbol{\nabla} h| } \right) , \quad \frac{1}{100} < |\boldsymbol{\nabla} h| < \frac{1}{2} . \end{equation}
Figure 2 displays the comparison of (2.12) and (2.24a,b) in the range 
$|\boldsymbol {\nabla } h| \in [1/100,1/2]$ and 
$\varepsilon \in [1/100, 1/10]$ covering the wave conditions in the numerical part of the present study. It shows that the approximated expression (2.24a,b) represents well the original net potential energy formulation within the considered range of slope magnitudes. Using 
$\exp ({ -48 \check {\mathscr {E}}_{p2}}) \approx 1\unicode{x2013}48 \check {\mathscr {E}}_{p2}$, the background kurtosis can be approximated as 
$10{\rm \pi} ^{2}\varepsilon ^{2}/\mu ^2$ and thus the excess kurtosis up to second order in steepness reads (see Appendix B of Mendes (Reference Mendes2024))
\begin{equation} \hat{\lambda}_{4} \approx 20 {\rm \pi}^{2} \left( \frac{\varepsilon}{\mu} \right)^{2} \sqrt{|\boldsymbol{\nabla} h|} ,\quad \textrm{for} \ \varepsilon \lesssim 1/20 . \end{equation}Hence, the difference of excess kurtosis between preshoal and postshoal condition can be estimated, so the simplified form of (2.15) reads
\begin{equation} \Delta \lambda_4 \approx 20 {\rm \pi}^{2} \left[ \left( \frac{\varepsilon}{\mu} \right)^{2} - \left( \frac{\varepsilon_0}{\mu_0} \right)^{2} \right] \sqrt{|\boldsymbol{\nabla} h|} , \quad \textrm{for}\ \mu_{0} \gg 1 . \end{equation}
For a special but representative case of a Gaussian preshoal wave distribution, namely 
$\hat {\lambda }_{4,0}=0$, we have 
$\Delta \lambda _4= \hat {\lambda }_{4}$. Noteworthy, due to the Taylor expansions, (2.25) will not be able to capture all features of the spatial evolution of the excess kurtosis, which are better described by (2.14). As such, (2.25) is more suitable for estimating the maximum excess kurtosis atop the shoal relative to the preshoal condition. Furthermore, this approximation and the computation of 
$\varGamma$ in (2.10) are restricted to both linear and second-order waves, as we have assumed that 
$H_{s}(x) \ll h(x)$. In the next section, our theory and its inference will be validated by confronting it with fully nonlinear simulation results.
3. Numerical method: FNPF model
In order to investigate the effects related to the shoaling length and to validate the predictions of the non-homogeneous second-order theory, we performed fully nonlinear numerical simulations within the framework of two-dimensional FNPF theory. The FNPF theory assumes that the fluid is inviscid and incompressible, and the flow is irrotational. In addition, in this study, the free surface tension is ignored, and the seabed elevation is fixed in time. The water-wave problem can then be formulated in terms of the velocity potential 
$\varPhi (x,z,t)$ and the FSE 
$\eta (x,t)$
$$\begin{align} \frac{\partial^2 \varPhi}{\partial x^2}+\frac{\partial^2 \varPhi}{\partial z^2} &=0 {,\quad {\rm in}}\ -h(x) \leqslant z \leqslant \eta(x,t) , \end{align}$$
$$\begin{align}\frac{\partial \eta}{\partial t}+\frac{\partial \varPhi}{\partial x}\frac{\partial \eta}{\partial x}-\frac{\partial \varPhi}{\partial z} &=0 {,\quad {\rm on}}\ z = \eta(x,t) , \end{align}$$
$$\begin{align}\frac{\partial \varPhi}{\partial t}+\frac{1}{2}\left[\left(\frac{\partial \varPhi}{\partial x}\right)^2+\left(\frac{\partial \varPhi}{\partial z}\right)^2\right]+g\eta &=0 {,\quad {\rm on}\ } z = \eta(x,t) , \end{align}$$
$$\begin{align}\frac{\partial \varPhi}{\partial x}\frac{\partial h}{\partial x}+\frac{\partial \varPhi}{\partial z} &=0 {,\quad {\rm on}\ } z ={-}h(x) . \end{align}$$ Zakharov (Reference Zakharov1968) and Craig & Sulem (Reference Craig and Sulem1993) have shown the FNPF problem can be expressed using only the variables on the free surface, i.e. 
$\eta (x,t)$ and 
$\tilde {\varPhi }(x,t)\equiv \varPhi (x,z=\eta,t)$, which is also known as the Zakharov formulation,
$$\begin{gather} \frac{\partial \eta}{\partial t} ={-}\frac{\partial \tilde{\varPhi}}{\partial x}\frac{\partial \eta}{\partial x}+ \tilde{w}\left[1+\left(\frac{\partial \eta}{\partial x}\right)^2\right] {,\quad {\rm on}\ } z = \eta(x,t) , \end{gather}$$
$$\begin{gather}\frac{\partial \tilde{\varPhi}}{\partial t} ={-}g\eta-\frac{1}{2}\left(\frac{\partial \tilde{\varPhi}}{\partial x}\right)^2+\frac{1}{2}\tilde{w}^2\left[1+\left(\frac{\partial \eta}{\partial x}\right)^2\right] {,\quad {\rm on}\ } z = \eta(x,t) , \end{gather}$$
where 
$\tilde w(x,t) \equiv \partial \varPhi /\partial {z}(x,z=\eta (x,t),t)$ denotes the vertical velocity of the fluid particles on the free surface. To evaluate the temporal evolution of 
$\eta$ and 
$\tilde {\varPhi }$, one needs to solve 
$\tilde w$ from the given free surface variables 
$(\eta, \tilde {\varPhi })$. This is the core of the Zakharov formulation and is known as the Dirichlet–to–Neumann (DtN) problem. Various approaches have been proposed to solve the DtN problem (see e.g. Dommermuth (Reference Dommermuth2000), Madsen, Fuhrman & Wang (Reference Madsen, Fuhrman and Wang2006), Bingham, Madsen & Fuhrman (Reference Bingham, Madsen and Fuhrman2009), Papoutsellis, Charalampopoulos & Athanassoulis (Reference Papoutsellis, Charalampopoulos and Athanassoulis2018) and references therein).
Recently, the exact FNPF problem has been resolved by using a spectral approach in the vertical direction and a high-order finite difference method in the horizontal direction, with a code called Whispers3D (W3D). This model has been described and applied in various scenarios by Yates & Benoit (Reference Yates and Benoit2015), Raoult, Benoit & Yates (Reference Raoult, Benoit and Yates2016), Simon et al. (Reference Simon, Papoutsellis, Benoit and Yates2019), Zhang & Benoit (Reference Zhang and Benoit2021) and Zhang et al. (Reference Zhang, Benoit and Ma2022), showing high fidelity in modelling highly nonlinear non-overturning waves for arbitrary bottom profiles. In this model, a change of the vertical coordinate from 
$z$ to 
$s$ is adopted such that the physical domain 
$(x,z)$ with varying free surface and bottom boundaries 
$[-h(x), \eta (x,t)]$ can be transformed into a new domain 
$(x,s)$ with fixed upper and lower boundaries 
$[-1, 1]$:
\begin{equation} s(x,z,t)=\frac{2z+h(x)-\eta(x,t)}{h(x)+\eta(x,t)} . \end{equation}
Then, the model equations are reformulated in the 
$(x,s)$ domain and the vertical profile of the velocity potential is approximated using a basis of orthogonal Chebyshev polynomials of the first kind 
$T_n(s)$, truncated at a tunable maximum order 
$N_T$:
\begin{equation} \varPhi(x,s,t) \approx \sum_{n=0}^{N_T}a_n(x,t) T_n(s) . \end{equation}
The water-wave problem is solved once the 
$N_T+1$ unknown coefficients 
$a_n(x,t)$ are determined at each abscissa 
$x$ from the DtN problem. The detailed reformulation of the governing equations in the 
$(x,s)$ domain, the Chebyshev-tau method used to solve the DtN problem and the numerical algorithm adopted in this model have been reported in the above-cited references, and are not duplicated here.
The parameter 
$N_T$ plays a crucial role in balancing the accuracy and efficiency of W3D, and should be adapted for different scenarios to achieve optimal model performance. To illustrate the role played by 
$N_T$, we present the comparison of the phase speed 
$C(\mu \equiv kh,N_T)$ of sinusoidal waves in uniform water depth 
$h$ predicted by the linear wave version of W3D model with the analytical phase speed 
$C_{Airy}(\mu )$ of Airy linear wave theory. As shown by Benoit, Raoult & Yates (Reference Benoit, Raoult and Yates2017), the linear W3D model yields an analytical solution of 
$C(\mu,N_T)$,
\begin{equation} \frac{C}{\sqrt{gh}} = \sqrt{\frac{1 + \sum_{n=1}^{N_T-2}{A_n\mu^{2n}}}{1+\sum_{n=1}^{N_T-1}{B_n\mu^{2n}}}} , \end{equation}
where the 
$A_n$ and 
$B_n$ can be computed analytically (Benoit et al. Reference Benoit, Raoult and Yates2017). The reference phase speed from Airy theory reads
\begin{equation} \frac{C_{Airy}}{\sqrt{gh}}=\sqrt{\frac{\tanh\mu}{\mu}} . \end{equation}
Figure 3 shows that the relative error between 
$C$ and 
$C_{Airy}$ for linear waves reduces as 
$N_T$ increases, and it is of low level for deep-water waves (for instance, it is below 
$10^{-5}$ with 
$N_T \geqslant 7$ when 
$\mu ={\rm \pi}$). This relative error remains acceptable even for extremely deep water when 
$\mu$ reaches 
$100$: it does not exceed 
$2.5\,\%$ with 
$N_T \geqslant 10$ and decreases down to approximately 
$1\,\%$ with 
$N_T=15$. It could be further improved at the expense of computational resources (by increasing 
$N_T$). Figure 3 also provides guidance for the selection of 
$N_T$ in different cases. The present study on the shoaling length parameter effects benefits from the W3D model with adjustable accuracy/efficiency, as the relative water depth varies significantly among the different simulation cases. In the meantime, it is worth mentioning that the W3D model is capable of handling very steep slopes. Successful applications are reported in Zhang & Benoit (Reference Zhang and Benoit2021) for nonlinear irregular waves passing over a slope with gradient 
$|\boldsymbol {\nabla }h|\approx 0.26$, and in Benoit et al. (Reference Benoit, Raoult and Yates2017) for linear regular waves over a slope with mean gradient approximately 
$1.42$. This makes W3D a particularly relevant choice for the present study.
Relative error of the phase speed predicted by the linear wave version of the W3D model with different values of 
$N_T$ in comparison with the analytical solution of Airy linear wave theory. Vertical dash lines denote the generally adopted shallow- and deep-water limits, 
$\mu = kh={\rm \pi} /10$ and 
$\mu ={\rm \pi}$, respectively.
4. Numerical simulations and results analysis
4.1. Configurations of the simulation cases
The key to choosing the wave parameters of the incident wave train and the bathymetry set-up lies in two aspects. On one hand, the incident wave train parameters are manipulated to keep the steepness 
$\varepsilon$ and relative water depth 
$\mu$ after the shoal nearly unchanged among the various cases, such that the sea-states after the shoal differ only because of different degrees of non-equilibrium dynamics induced by different lengths of the shoal. On the other hand, the slope magnitude and shoaling length are disentangled by keeping 
$|\boldsymbol {\nabla }h|$ constant and varying only the slope length 
$L$. Furthermore, we choose a bottom profile with an upslope only but no downslope to avoid any confusion due to the disturbances of the reflected wave energy during the deshoaling process (Zhang & Benoit Reference Zhang and Benoit2021).
The anomalous wave statistics are very sensitive to the relative depth after the shoal 
$\mu _f = k_f h_f$ (Trulsen et al. Reference Trulsen, Zeng and Gramstad2012). Inspired by run 3 of the experiments of Trulsen et al. (Reference Trulsen, Raustøl, Jorde and Rye2020) and the simulations of Zhang et al. (Reference Zhang, Benoit and Ma2022), we use the water depth after the shoal 
$h_f=0.11$ m and incident peak period 
$T_p=1.1$ s in all cases, thus the relative water depth in the shallower region 
$\mu _f$ is also a constant. As a consequence, the preshoal water depth 
$h_0$ is fully determined by the value of the slope magnitude 
$|\boldsymbol {\nabla } h |$ and the shoaling length 
$L$. We acknowledge that the relative water depth difference 
$\mu _f/\mu _0$ also plays an effect on the anomalous wave statistics. Here we choose to fix 
$\mu _f$ and vary 
$\mu _0$, because the role played by 
$\mu _0$ is of secondary importance compared with that of 
$\mu _f$, and only becomes non-negligible when 
$\mu _0 \rightarrow \mu _f$. However, as 
$\mu _0 \rightarrow \mu _f$, wave shoaling is no longer the dominant physics and cannot excite anomalous wave statistics.
In theory, 
$\ell$ varies as 
$\lambda$ evolves over the slope and could be computed locally in space. However, we need a characteristic measure 
$\ell$ for the entire shoal when comparing the effect of varying 
$L$. Hereafter, we use its value at the end of the shoal, 
$\ell _{p} \equiv L/\lambda _{p,f}$, to better show the correspondence between the enhancement of kurtosis and the shoaling length parameter. As a result, the change in bathymetry now boils down to a longer shoal and a deeper water depth 
$h_0$ prior to the shoal, as illustrated in figure 4(a). Suppose the slope length 
$L$ is extended by a factor of 
$n$ to 
$L' = nL$, resulting in a new deeper region depth 
$h_0'$. Since 
$|\boldsymbol {\nabla } h|$ is kept the same in both cases, i.e. 
$(h_0'-h_f)/L' = (h_0-h_f)/L$, the preshoal depth is scaled correspondingly,
\begin{equation} \frac{h_0'}{h_0}=n \left[ 1 + \left(\frac{1}{n}-1\right)\frac{h_f}{h_0} \right] , \quad n \in \mathbb{R}^{{\ast}}_+ . \end{equation}
As demonstrated by Trulsen et al. (Reference Trulsen, Raustøl, Jorde and Rye2020), as long as 
$h_0$ is in deep water 
$(\mu _0 \gtrsim {\rm \pi})$, any increase in the water depth will have no significant impact on the anomalous wave statistics before and atop the shoal. This can be further understood from the context of modulational instability (Zakharov & Ostrovsky Reference Zakharov and Ostrovsky2009). By increasing the relative water depth 
$\mu _0 = k_{p0} h_0$, the initial value of excess kurtosis will increase by a few per cent due to the Benjamin–Feir index, see for instance the Benjamin–Feir index dependence on 
$\mu$ in Zhang et al. (Reference Zhang, Guedes Soares, Cherneva and Onorato2014). However, making the relative water depth 
$\mu _0 \gtrsim {\rm \pi}$ larger will not affect the wave statistics atop the shoal, and the difference between peak and preshoal statistics will be only slightly (by a few per cent) decreased. As the change of slope length will lead to a different shoaling coefficient (defined as 
$C_{shoal} = H_{s,f}/H_{s,0}$), the incident 
$H_{s,0}$ is tuned in each case, such that 
$\varepsilon _f$ atop the shoal is kept the same regardless of the length of the shoal.
Sketches of (a) the necessary change in the bottom profile for assessing the role played by the shoaling length effect and (b) the numerical wave flume (NWF) (not to scale).
We consider long-crested irregular wave trains described by a JONSWAP (joint North Sea wave project) spectrum,
\begin{equation} S(f)=\frac{\alpha_J{g}^2}{\left(2{\rm \pi}\right)^4}\frac{1}{f^5}\exp{\left[-\frac{5}{4}\left(\frac{f_p}{f}\right)^4\right]} \gamma^{\exp{\left[-{\left(f-f_p\right)^2}/{2\left(\sigma_Jf_p\right)^2}\right]}} , \end{equation}
where 
$\alpha _J$ denotes the adjustment factor for 
$H_s$, 
$\sigma _J$ the spectral asymmetry parameter (
$\sigma _J=0.07$ for 
$f \leqslant f_p$ and 
$\sigma _J=0.09$ for 
$f>f_p$), and 
$\gamma$ the peak enhancement factor. In this study, the spectral peak period 
$T_p=1/f_p=1.1$ s and peak enhancement factor 
$\gamma =3.3$ of the incident wave field are the same for all 11 cases. In each case, the incident wave train lasts for 5060 s, thus 
$4600T_p$. Table 1 summarizes the configurations of the incident wave fields and bathymetry information in 11 cases, as well as the key wave parameters in both deeper and shallower regions (by averaging the simulation results over the corresponding areas). These cases are considered representative for investigating the shoaling length effect, as they cover a relatively wide range of shoaling length parameter, 
$\ell _p \in [0.1, 11.9]$. Meanwhile, the relative water depth difference before and after the shoal also varies significantly, 
$\mu _0/\mu _f \in [1.1, 17.9]$.
Summary of the key parameters for the simulations. The incident sea-states are described by a JONSWAP spectrum of the same peak period 
$T_p=1.1$ s, and peak enhancement factor 
$\gamma = 3.3$. The slope length changes from case 1 to 11, yet the slope is kept constant 
$|\boldsymbol {\nabla }h| = 0.2625$ (approximately 
$1:3.81$). The incident 
$H_{s,0}$ is tuned in each case to keep 
$H_{s,f}$ more or less the same. The steepness measure can be converted to other common definitions through 
$\varepsilon = (\sqrt {2}/{\rm \pi} )k_pH_s$. Note that case 7 here shares the same wave and bottom configurations as case 4 in Zhang et al. (Reference Zhang, Benoit and Ma2022).
The incident wave train is constructed by linearly superimposing the harmonic components of the prescribed spectrum, and imposed at the wave-making boundary of the NWF. The low-energy components of the incident spectrum are ignored, keeping the non-trivial ones in the range 
$[f_{min}, f_{max}]=[0.5f_p,4f_p]$, as justified in the next subsection. The sketch of the NWF is provided in figure 4(b): the water depth 
$h_0$ is constant in the flat area prior to the shoal (6 m in length, roughly 
$3\lambda _{p,0}$), then it decreases to 
$h_f = 0.11$ m due to the presence of an up-slope (with a constant slope magnitude 
$|\boldsymbol {\nabla } h| =0.2625$ and length 
$L$), and remains constant in the flat area (27 m in length, roughly 
$25\lambda _{p,f}$) after the shoal. Such a long after-shoal flat area allows the out-of-equilibrium sea-state, induced by the depth change, to re-establish a new equilibrium state, based on the work in Zhang et al. (Reference Zhang, Benoit and Ma2022). In addition, the computation domain comprises a generation relaxation zone of 6 m in length (
${\approx }3\lambda _{p,0}$) on the left-hand end of the NWF and an absorption relaxation zone of 21 m in length (
${\approx }20\lambda _{p,f}$) on the right-hand end. As detailed in Appendix A, such a length of the damping zone is needed to keep negligible reflection of all wave components, including the low-frequency (LF) ones, in spite of the higher cost of computation effort.
4.2. Result convergence and statistical variability
Case 7 in the present study shares the same wave and bottom conditions as case 4 in Zhang et al. (Reference Zhang, Benoit and Ma2022), and other cases are also very similar to it (see table 1). Therefore, it is logical to adopt the same spatial grid interval and time step, 
${\rm d}\kern0.07em x = 0.01$ m and 
${\rm d}t=0.01$ s, which result in the Courant–Friedrichs–Lewy number 
$\text {CFL}_i = C_{p,i}\, \text {d}{t}/ \text {d}\kern0.07em {x} = (\lambda _{p,i}/ \text {d}\kern0.07em{x})/(T_p/ \text {d}{t})$, with 
$i=0$ or 
$f$, to be 
$\text {CFL}_0 =1.64$ prior to the shoal and 
$\text {CFL}_f=0.97$ atop the shoal.
The investigation of the shoaling length effect poses high requirements on the model accuracy. The increase of slope length results in larger 
$\mu _0$, and the demands on model accuracy increase sharply to describe the wave evolution in the deeper flat region. This is because, for irregular waves, the model should be reasonably accurate up to the cut-off frequency 
$f_{max}$ of the input spectrum, which is much more challenging than considering model accuracy up to the spectral peak frequency 
$f_p$. Table 1 lists the relative water depth 
$\mu (4f_p,h_0)$ for waves with 
$f=4f_p$ in depth 
$h_0$. In case 11 for instance, the relative water depth corresponding to 
$f_{max}$ in the NWF is actually higher than 180 prior to the shoal. We know from figure 3 that the accuracy of W3D is sensitive to the relative water depth, and that larger 
$N_T$ is needed for deeper water. We therefore performed a series of convergence tests for determining the cut-off frequency of the input spectrum and the choice of 
$N_T$ value, both would influence the performance of W3D.
We determined the extent of the spectral range 
$[0.5f_p,f_{max}]$ required to ensure accurate convergence of the simulations by running case 7 with 
$f_{max}$ taking 
$2f_p$, 
$3f_p$, 
$4f_p$ and 
$5f_p$ successively. Figure 5(a) shows the comparison of the spatial evolution of 
$\lambda _4$ with these four different cut-off frequencies. It is observed that the estimation of kurtosis is convergent for 
$f_{max} \geqslant 3f_p$. Figure 5(b) shows the relative error of 
$\eta (x)$ at 
$t=5060$ s with respect to the result with highest 
$f_{max}=5f_p$, and normalized by the incident significant wave height 
$H_{s,0}$. It is noticed that the convergence of FSE in the time domain is closely achieved for 
$f_{max}=4f_p$ with relative error within 1 %. Therefore, 
$f_{max}=4f_p$ is chosen for all subsequent simulations. Furthermore, the evolution of 
$\lambda _4$ in the case with 
$f_{max}=2f_p$ differs only slightly from the well-converged results with 
$f_{max} \geqslant 4f_p$. This seems to be contradictory to the conclusion of Tang et al. (Reference Tang, Barratt, Bingham, van den Bremer and Adcock2022), which states that removing the high-frequency spectral tail would result in substantial enhancement of 
$\lambda _4$. In fact, this is because the mechanisms are different in the study of Tang et al. (Reference Tang, Barratt, Bingham, van den Bremer and Adcock2022) and in the present one: in the former work, the spectral change manifests spontaneously under the flat bottom condition due to strong nonlinear wave–wave interaction. Whereas, in the present one, the spectral modulation is forced by the depth transition, and thus is less sensitive to the high-frequency tail of the incident spectrum but more to the rapid depth change.
(a) Spatial evolution of kurtosis 
$\lambda _4$ for case 7 considering four increasing cut-off frequencies 
$f_{max}$ of the incident JONSWAP wave spectrum; (b) relative error of 
$\eta (x)$ at final time 
$t = 5060$ s for the simulations with 
$f_{max} = 2 f_p, 3 f_p$ and 
$4 f_p$, with respect to the simulation with 
$f_{max} = 5 f_p$ used as reference.
All the cases listed in table 1 have been first tested for a short duration 
$50$ s (approximately 
$45 T_p$) to tune the 
$N_T$ values so that the balance between efficiency and precision can be achieved in each case. In figure 6(a), the tuned values of 
$N_T$ are displayed for all 11 cases. In figure 6(b), the calibration of 
$N_T$ in case 7 is provided as an example, it is observed that the time evolution of 
$\eta$ is identical for 
$N_T = 8$ and 
$10$, therefore in the final simulation of case 7 with long duration, 
$N_T=8$ is chosen.
(a) Calibrated values of 
$N_T$ for cases 1 to 11; (b) comparison of snapshots of duration 
$5 T_p$ of the time series of FSE 
$\eta$ at 6 m after the shoal in case 7 simulated with 
$N_T= 6$, 8, 10.
The skewness and kurtosis, as high-order moments, are subject to statistical variability. The duration of the record of 
$\eta$ plays an important role in obtaining (at least nearly) converged results of skewness and kurtosis. We checked that a time duration of 5060 s (equivalently, 4600
$T_p$) in the simulations is sufficient for obtaining statistically converged estimates of skewness and kurtosis. For such a purpose, the convergence of skewness and kurtosis is analysed with the time series computed at the location 
$x=0.75$ m where they assume their maximum values and are subject to the strongest variability. The time series of 
$\eta$ is separated into a series of individual waves through the zero-down-crossing method. Then, by evaluating the skewness and kurtosis of a section of the entire time series that contains the first 
$N$ waves, the skewness and kurtosis are functions of the number of waves in the time series. In figure 7, 
$\lambda _3$ and 
$\lambda _4$ computed from different numbers of waves in the time series are shown; it is noted that the estimates become nearly stable after including 4000 random waves, thus the selected simulation duration is long enough.
Convergence of (a) skewness and (b) kurtosis of FSE at 
$x=0.75$ m (equivalently 
$x/\lambda _{p,f} = 0.7$, where the skewness and kurtosis assume their maximum values) as a function of the number of waves in the simulated time series in cases 1 to 11.
4.3. Shoaling length effects on the spectral evolution along the NWF
The spatial evolution trends of two key wave parameters, the relative water depth 
$\mu$ and the wave steepness 
$\varepsilon$, are displayed in figure 8. They are evaluated locally, based on the spectral peak frequency and significant wave height extracted from spectral analysis of the time series of FSE. This figure further validates the methodology for determining the wave parameters of the input wave field, which ensures that 
$\mu _f$ and 
$\varepsilon _f$ after the shoal are basically the same whatever 
$\ell$. Before the shoal, 
$\mu$ varies significantly due to the change of depth 
$h_0$ among the cases, and 
$\varepsilon$ also differs in order to compensate for different shoaling effects on 
$H_{s,f}$ among the cases.
Spatial evolution along the NWF of non-dimensional wave parameters for cases 1 to 11: (a) relative water depth 
$\mu$ and (b) wave steepness 
$\varepsilon$. The vertical dash line at 
$x=0$ indicates the starting location of the shallower flat region.
The FSE time series is saved every 
$0.2$ m (0.2
$\lambda _{p,f}$ approximately) in the NWF, resulting in a relatively fine resolution of the spectral evolution in space, as is displayed in figure 9. It is noted that, in all cases, the spectral evolution after the shoal is almost the same: a beating pattern manifests for the second harmonics shortly after the shoal then gradually vanishes as waves propagate farther, and a broadened spectrum is eventually established. It indicates that the shoaling length effect is insignificant for wave spectral evolution as waves pass over a steep shoal. In cases 1 to 4, the spectral evolution is slightly different from the other cases: the beating pattern appears not only in the range after the slope but also in the flat area before it. The relative water depths before the shoal 
$\mu _0$ of these cases are low among all cases (see table 1) whilst the steepness values 
$\varepsilon _0$ are higher than other cases. Take case 1 as an example, since the Ursell number is proportional to 
${Ur}_0 \propto \varepsilon _0 / \mu _0^3$ and the relative water depth of case 1 has been decreased tenfold and the steepness increased by 50 % in comparison with case 11, case 1 has an Ursell number a thousand times larger than that in cases 9–11, resulting in higher relative importance of wave nonlinearity than dispersion.
Spatial evolution along the NWF of the wave spectrum for cases 1 to 11, displayed in (a–k), respectively. The colour scale depicts the value of log
$_{10}(S(f,x))$ with spectrum 
$S(f,x)$ in 
${\rm m}^2\ {\rm Hz}^{-1}$. The vertical white dash lines represent the extent of the plane slope.
The reflection is relevant in the simulations, it could be induced by the shoal and by the end of the damping zone. As the latter is considered as contamination to the wave field over the shoal, we have chosen a long damping zone to minimize its effect, see Appendix A for more detail. The reflection by the shoal is physical and the reflection rate in the simulations is around 8 % as in the simulations of Zhang & Benoit (Reference Zhang and Benoit2021). Furthermore, the interaction between the incident and reflected waves before the shoal affects the spectral peak frequency and results in the oscillatory behaviour of 
$\varepsilon$ prior to the shoal in figure 8(b).
4.4. Shoaling length effects on the spatial evolution of statistical parameters
Here, the spatial evolution of the asymmetry parameter, skewness and kurtosis is investigated, and the influences of the 
$\ell$ parameter on the evolution of these key statistical parameters are discussed. The asymmetry parameter measures the wave profile asymmetry in the horizontal direction. It is computed as the skewness of the Hilbert transform of the FSE, 
$\lambda _3[\mathcal {H}(\eta )]$, with 
$\mathcal {H}$ denoting the Hilbert transform operator. The asymmetry parameter assumes negative values when the wave front is steeper than its rear face, for instance, when a wave passes over a shoal. The results of the three statistical parameters are displayed in figure 10(a–c).
Spatial evolution of (a) asymmetry parameter 
$\lambda _3[\mathcal {H}(\eta )]$, (b) skewness 
$\lambda _3(\eta )$ and (c) kurtosis 
$\lambda _4(\eta )$, in cases 1 to 11 along the NWF. The vertical solid line at 
$x=0.7\lambda _{p,f}$ indicates the position where skewness and kurtosis achieve their maxima.
In figure 10(a), the evolution of the asymmetry parameter is displayed for all cases. The horizontal wave profile first leans forward due to shoaling (with 
$\lambda _3[\mathcal {H}(\eta )]$ approaching its negative extreme value) in the short area near the end of the slope, and then the wave profile develops reversely (with 
$\lambda _3[\mathcal {H}(\eta )]$ returning to 0 around 
$x=0.7\lambda _{p,f}=0.75$ m on its way to the positive maximum), and becomes backwards-leaning due to significant non-equilibrium wave evolution. Eventually, 
$\lambda _3[\mathcal {H}(\eta )]$ approaches a positive constant in the equilibrium state that is nearly established close to the end of the NWF. All cases 5 to 11 follow this trend, with their positive and negative global extreme values and the equilibrium values being almost the same and achieved at the same location. The picture is a bit different in cases 1 to 4, where the asymmetry parameter is higher than 0 initially due to the higher relative importance of nonlinearity as illustrated previously.
The spatial evolution of skewness 
$\lambda _3(\eta )$ in all cases is shown in figure 10(b). It has a finer resolution compared with figure 10(a), because W3D allows the computation of skewness and kurtosis at every grid point (with 
${\rm d}\kern0.07em x=0.01$ m) during the simulation. As for the asymmetry parameter, the evolution trends are very similar in all cases, with their maximum values atop the shoal converging at the same location (convergent maximum value 
$\lambda _3 \approx 0.6$). Furthermore, the equilibrium values of 
$\lambda _3$ in the far field atop the shoal are also at the same level, as the sea-states in all cases are of nearly the same levels of nonlinearity and dispersion. In cases 1 to 4, the skewness in the deeper flat region before the slope is higher than that in other cases, this is again related to the higher relative importance of wave nonlinearity.
In figure 10(c), the spatial evolution of kurtosis in all cases is shown. In all cases, 
$\lambda _4$ remains around 3 before the end of the shoal, and after the shoal, 
$\lambda _4$ is significantly enhanced to its local maximum at 
$x=0.7\lambda _{p,f}=0.75$ m, indicating a local increase of the rogue wave occurrence probability. Then, 
$\lambda _4$ declines to a level around 2.6, which indicates that the probability of rogue waves in the new equilibrium state is even lower than that in a Gaussian sea-state. We notice that 
$\lambda _4$ may not be converged yet at 
$x=25\lambda _{p,f}$. This is understandable because as a fourth-order moment, kurtosis would require a longer distance to be convergent than the third-order ones. As the trend of kurtosis is already clear in figure 10(c), we decided not to extend the length of NWF. Based on the evolution of these three statistical parameters shown in figure 10, we consider that the short-scale non-equilibrium wave evolution stage happens in the range 
$[0,5]\lambda _{p,f}$ and the long-scale in 
$[5,25]\lambda _{p,f}$ in the NWF. The spatial extent of these two scales is independent of 
$\ell _p$ (or of 
$\ell$, equivalently).
Now we focus on the local extremes of these statistical moments. From figure 10(a), it is noticed that the local extreme values of 
$\lambda _3[\mathcal {H}(\eta )]$ after the shoal in cases 1 and 2 are of lower modulus than other cases. From figures 10(b) and 10(c), we see that in cases 1 to 4, 
$\lambda _3$ and 
$\lambda _4$ have the lower maximum values shortly after the shoal. Furthermore, we observe that, despite an initial offset, the changes of the values before and after the shoal of the three variables (asymmetry, skewness, kurtosis) are quite small for the cases 1–4. Such behaviour in cases 1–4 is expected because the depth transition is not as strong as in cases 5–11, see the last column of table 1. Consequently, the depth changes in cases 1 to 4 are in a transitional regime between homogeneous evolution of the wave field (constant depth condition) and inhomogeneous evolution (rapidly depth varying condition). When the depth transition is strong enough, the wave propagation after the shoal is dominated by the non-equilibrium dynamics induced by the change of water depth. As a result, neither the preshoal values nor the maximum of both skewness and kurtosis feature any significant differences due to the change of 
$\ell _p$ in cases 5–11.
In figure 11, the change of excess kurtosis 
$\Delta \lambda _4$ between the over-shoal maximum value and the preshoal mean value is displayed as a function of 
$\ell _p$. The empirical kurtosis change in the simulations is obtained from the maximum value of 
$\lambda _4$ in the range 
$x\in [0,5]\lambda _{p,f}$ with its mean value in the range 
$x\in [-L-6 \text {~m},-L]$ subtracted. The theoretical predictions of 
$\Delta \lambda _4$ can be obtained with two levels of approximation:
(i) with
$\Delta \lambda _4$ in (2.15) and $\varGamma$ in (2.10), originally put forward in Mendes & Kasparian (Reference Mendes and Kasparian2023);(ii) with the simplified expression of
$\Delta \lambda _4$ in (2.26), as an approximation to option (i).
As a remark, the 
$\ell$ is computed from a zero-crossing wavelength from the perspective of (2.10) that would give rise to the kurtosis of models (i) and (ii), whereas the simulations have a peak wavelength counterpart. Following Figueras (Reference Figueras2010) for the relation between peak and mean periods, as well as Mendes & Scotti (Reference Mendes and Scotti2021) for the relation between mean and zero-crossing period, we may approximate 
$T_{p}^{2} \approx 2 T_{z}^{2}$ (with 
$T_z$ denoting the zero-crossing mean wave period) such that 
$\ell \sim 2 \ell _{p}$.
Comparison of the kurtosis enhancement 
$\Delta \lambda _4$ atop the shoal between simulated results and theoretical estimation with three options, displayed as functions of the shoaling length parameter 
$\ell _p$. The simulation results are marked as hollow circles, and the kurtosis predictions are computed in two ways (as indicated in the legend box).
For a strong depth transition (
$\mu _0 / \mu _f \gtrsim 2$ in cases 5–11), the two theoretical models provide similar results, both matching the simulations. Indeed, considering that for a total of 5000 waves (see figure 7), the expected 95 % confidence interval (
$\pm 2 \sigma$) has a width of 
$\pm 0.15$ (Joanes & Gill Reference Joanes and Gill1998; Cramér Reference Cramér1999; Wright & Herrington Reference Wright and Herrington2011), the two theoretical curves stay within this range near the simulated values. Thus, their agreement with the simulated excess kurtosis is excellent, which makes them an effective theory for the shoaling length variations, especially as option (ii) is of a much simpler formulation. In this range, 
$\Delta \lambda _4$ shows virtually no dependence on 
$\ell _p$ in the two theoretical models and the simulations. A relevant pondering arises regarding what would happen if we kept increasing the shoaling length beyond the range plotted in figure 11. Overall, the ever-increasing of 
$\ell _p$ will not change 
${\Delta }{\lambda }_4$ because the initial relative water depth 
$\mu _0$ would be so large (
$\mu _0>12$) that the waves cannot feel the ever-deeper water region, as discussed in § 4.3.
In the weak depth transition regime (
$\mu _0 / \mu _f \lesssim 2$, in cases 1–4), the two theoretical formulations are able to predict 
$\Delta \lambda _4$ in this regime as well. Although at first glance, one may interpret the growth of the excess kurtosis as being due to an increasing 
$\ell _p$, this growth is rather related to the strengthening of the depth transition because of the way we defined the change in water depth in (4.1). This is understandable, if one considers a case with 
$\mu _0 \rightarrow \mu _f$, no matter how long the shoal (with an infinitesimal gradient) is or how steep the shoal (with an infinitesimal length) is, waves would propagate without significant change of excess kurtosis, as in the flat bottom condition. In other words, in such a case, the dominant factor is neither the shoaling length nor the slope, but 
$\mu _0 / \mu _f$. Therefore, we conclude that the shoaling length parameter plays an insignificant role in the change of excess kurtosis in the weak depth transition regime.
5. Conclusion
In this work, we have demonstrated both theoretically and numerically that the enhancement of the excess kurtosis of FSE due to the shoal is not sensitive to the shoaling length. This has been shown by setting fixed values for wave steepness and relative water depth after the shoal and setting a constant slope magnitude. We provided an explicit definition of the abruptness of a shoal by introducing a shoaling length parameter 
$\ell _p$ as the ratio between the slope length and the characteristic wavelength. It allowed us to disentangle the effect of the shoaling length parameter from that of the slope magnitude by varying the preshoal depth while keeping the slope magnitude constant. By building on the recent theoretical work in Mendes et al. (Reference Mendes, Scotti, Brunetti and Kasparian2022), we derived estimates of the change of excess kurtosis from the second-order theoretical model with two levels of approximations. Both showed good agreement with the simulation results, indicating that the enhancement of kurtosis (i.e. rogue wave amplification) due to shoaling is independent of the shoaling length, as long as the wave evolution is dominated by the non-equilibrium dynamics due to the change of water depth. We therefore deduce that, for a fixed wave field condition after the shoal, the magnitude of the bottom slope is the main driver of the enhancement of the excess kurtosis.
Then, we showed numerically that not only the rogue wave amplification (characterized by the peak values of the statistical moments) but also the location of these maxima, the vertical asymmetry in surface elevation, wave steepness and spectral evolution are insensitive to the shoaling length. We observed that shorter shoaling lengths (typically 
$\ell _p \leqslant 0.5$) lead to a reduced increase of kurtosis over the shoal. This is because the depth prior to the shoal is insufficient to induce non-equilibrium dynamics. For the non-equilibrium wave evolution induced by depth transitions, three stages are characterized by different features: (1) the near-equilibrium evolution stage before the end of the slope; (2) the short-scale non-equilibrium evolution stage from the end of the up-slope up to a small distance atop the shoal, often called latency (Zheng et al. Reference Zheng, Lin, Li, Adcock, Li and van den Bremer2020); (3) the long-scale non-equilibrium evolution stage as waves propagate farther before reaching a new equilibrium state. Our FNPF simulation results have shown that the shoaling length parameter has no discernible effect on the spatial extent of short- and long-scale non-equilibrium wave evolution, nor on the locations where global extreme values or the equilibrium values of the statistical parameters are achieved.
From a more technical point of view, when designing the set-up of the FNPF simulations, we have shown that (i) a cut-off at higher frequencies of the wave spectrum has a significantly smaller impact on wave statistics than expected for a flat bottom in deep water, and (ii) insufficient attenuation of LF waves at the downstream boundary of the domain could have a notable influence on wave statistics and spectrum atop the shoal. Dedicated sensitivity studies and choices have been done to get rid of these two effects in the final set of simulations.
Finally, our study covers the range 
$0.1 \leqslant \ell _p \leqslant 12$, but equally applies beyond. Indeed, longer shoals correspond to even deeper preshoal conditions, where the waves do not interact with the bottom. The shoal therefore only starts to affect the propagation when the waves quit the deep-water region, thus limiting the effective length of the shoal. Conversely, very short 
$\ell _p$ implies negligible depth difference, hence, a vanishing shoal effect.
The outcomes of this work could also be useful for researchers investigating either experimentally or numerically the depth-induced non-equilibrium dynamics when designing the incident wave and bathymetry configurations. Also, the explicit and simple expression for the change of excess kurtosis of (2.25) can be used for fast estimation of freak-wave risk associated with a steep shoal in engineering practices. Further work would, however, be needed to assess the effect of 
$\ell$ on waves of larger steepness (beyond the second-order regime), especially reaching close to wave breaking.
Funding
In this work, J.Z. was supported by the National Natural Science Foundation of China (grant no. 52101301), and the China Postdoctoral Science Foundation (grant nos. 2023T160078, 2021M690523). S.M. and J.K. were supported by the Swiss National Science Foundation (grant nos. 200020-175697). We thank the three anonymous reviewers for their insightful suggestions that have significantly improved this manuscript.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Effects of LF waves in numerical simulations
Low-frequency subharmonics often occur as a result of nonlinear wave–wave interactions during nonlinear simulations. Particular attention should be paid to LF wave damping while running simulations with long duration, as LF wave energy may accumulate and affect the fundamental harmonics around 
$f_p$ in return. Conventionally, in W3D, a damping zone of three characteristic wavelengths at the end of the computational domain works well in absorbing waves. However, in the present work, as the duration is long and the nonlinear wave–wave interaction is quite active after the shoal, taking 
$L_{damp} = 4\ {\rm m} \approx 3 \lambda _{p,f}$ is not sufficient. Rather, the length of the damping zone should be set according to the LF wavelength, i.e. 
$L_{damp} = 21\ {\rm m} \approx 20 \lambda _{p,f} \approx 3 \lambda _{LF}$, with the LF waves having a characteristic frequency 
$f_{LF} \approx 0.05$ Hz.
Figure 12 displays the spatial evolution of the wave spectrum for case 7 with 
$L_{damp} \approx 3 \lambda _{p,f}$ (in figure 12a) and 
$L_{damp} \approx 20 \lambda _{p,f}$ (in figure 12b). The vertical dash lines indicate the range of the shoal, and the horizontal dash lines represent spectral peak frequency evaluated as 
$f_p \equiv ( \int _{0}^{f_{Nyq}} fS^4(f) \,\textrm {d}f ) / ( \int _{0}^{f_{\textrm {Nyq}}} S^4(f) \,\textrm {d}f )$ with 
$f_{Nyq}=50$ Hz according to the Nyquist sampling theorem. Clearly, with a short damping zone, the LF components receive a considerable amount of energy during the simulation and affect markedly the peak frequency atop and after the shoal in figure 12(a). With the longer damping zone, however, as shown in figure 12(b), the LF components are effectively suppressed. Consequently, the longer damping zone 
$L_{damp} \approx 20 \lambda _{p,f}$ was adopted in all simulations presented in this study.
Spatial evolution of the wave spectrum for case 7, simulated with two different extents of the damping zone after 
$x/\lambda _{p,f}=25$: (a) 
$L_{damp} = 4\ {\rm m} \approx 3 \lambda _{p,f}$ and (b) 
$L_{damp} = 21\ {\rm m} \approx 20 \lambda _{p,f}$. The colour scale depicts the value of 
$\log _{10}(S(f,x))$ with spectrum 
$S(f,x)$ in 
${\rm m}^2\ {\rm Hz}^{-1}$. The vertical white dash lines represent the extent of the plane slope. The horizontal dash line represents the local peak frequency of the wave spectrum.
Moreover, with these two choices of 
$L_{damp}$, the effect of LF components on the statistical parameters 
$\lambda _3$ and 
$\lambda _4$ can be discussed. In figure 13, the spatial evolutions of 
$\lambda _3$ and 
$\lambda _4$ are displayed. It is observed that strong LF components (i.e. with the reduced 
$L_{damp}$) result in lower levels of both parameters after the shoal. This means that, in nature, the LF components of coastal waves (possibly reflected from the shoreline or released during depth-induced wave breaking) could play a role in the mitigation of rogue-wave risk.
Spatial evolution along the NWF of (a) skewness and (b) kurtosis for case 7, simulated with two choices of relaxation zone lengths.
Open access
