3.1. Oblique incident waves

Here, we focus on the wave propagation direction $\theta _{0}$ . In Lyu et al. (Reference Lyu, Mori and Kashima2023), it was assumed that $\theta _{0}$ is fixed at zero, which means the principal wave direction aligns with the direction of bottom topography variation, i.e. predominantly along the $x$ axis. Meanwhile, wave components also diverge in directions oblique to the principal direction. The direction of each component can be quantified by the ratio $k_{\zeta }/k_{0}$ , where $k_{\zeta }$ is the lateral (i.e. $y$ axis, and $k_{\zeta }$ and $k_{y}$ are equivalent in the following) component of the carrier wavenumber $k_{0}$ . Accordingly, the propagation direction of each sub-wave can be expressed as $\theta _{\zeta }$ = arctan $({k_{\zeta }}/{k_{0}})$ . In figure 1(a), we give the absolute value of $\ddot{A}$ at $x=x_{0}$ on $\omega _{\tau },k_{\zeta }$ , where $\tan \theta _{0}=0$ . In the absence of oblique components, the initial wave energy (approximately represented as $| \ddot{A}|$ ) is distributed as a two-dimensional Gaussian centred at the carrier frequency $\omega _{0}$ and $k_{\zeta }=0$ , symmetric in $\omega _{\tau },k_{\zeta }$ .

Figure 1.

Initial condition and wave refraction: (a) Fourier spectral at $x_{0}$ with $\tan \theta _{0}=0$ ; (b) Fourier spectral at $x_{0}$ , $\omega _{0}$ with different $\tan \theta _{0}$ ; (c) evolution of $\tan \theta$ on a slope for different $\tan \theta _{0}$ ; (d) validation by Snell number on a slope with different $\tan \theta _{0}$ . (Red line, bottom bathymetry $-kh$ in the principal direction.)

Figure 1. Long description

Panel A: A 3D surface plot shows the Fourier spectral at a specific condition. The x-axis is labeled k_y/k_x, the y-axis is labeled omega/omega_0, and the z-axis is labeled with a magnitude symbol. The plot features a peak in the center, indicating a concentration of spectral energy. Panel B: A line graph displays the Fourier spectral at another condition with different tan theta_0 values. The x-axis is labeled k_y/k_x, and the y-axis is labeled with a magnitude symbol. Multiple lines represent different tan theta_0 values, showing varying spectral intensities. Panel C: A line graph illustrates the evolution of tan theta on a slope for different tan theta_0 values. The x-axis is labeled x/L_0, the left y-axis is labeled tan theta, and the right y-axis is labeled -kh. Different lines represent different tan theta_0 values, showing how tan theta evolves over the slope. Panel D: A line graph validates the Snell number on a slope with different tan theta_0 values. The x-axis is labeled x/L_0, the left y-axis is labeled C_s(x)/C_s(x_0), and the right y-axis is labeled -kh. Different lines represent different tan theta_0 values, showing the relationship between the Snell number and the slope.

When $\theta _{0}$ is non-zero, the directional distribution of wave energy is centred along an axis rotated by $\theta _{0}$ from the $x$ axis. In figure 1(b), we present cross-sections of the two-dimensional Fourier spectrum at $\omega _{\tau }=\omega _{0},x=x_{0}$ , corresponding to different values of $\theta _{0}$ . As shown, changes in $\theta _{0}$ lead to a shift of the initial spectrum on the $k_{\zeta }$ direction, while keeping the spectral shape the same. As the simulation progresses, the wave train propagates primarily along the direction defined by $\tan \theta = k_{y}/k_{x}$ , where total wavenumber $k=\sqrt{{k}_{x}^{2}+{k}_{y}^{2}}$ , $\theta =\theta _{0}$ at $x=x_{0}$ . It is important to note that the effective step ${\rm d}x'$ now should be the hypotenuse in the computational grid, i.e. ${\rm d}x'=\sqrt{{\rm d}x^{2}+{\rm d}y^{2}}$ , where ${\rm d}y/{\rm d}x=\tan \theta$ . When ${\rm d}y/{\rm d}x$ is sufficiently small – i.e. when the wave incidence angle is relatively small – ${\rm d}x'$ can be approximated by ${\rm d}x$ . Under this assumption, (2.4) remains applicable as the governing equation for wave envelope evolution.

Based on this framework, we can also simulate the refraction of wave trains as they enter shallower regions. As the water depth $h$ gradually decreases, the wavenumber component normal to the bottom topography increases in accordance with the dispersion relation. In this study, the bottom is simplified as a slope varying only along the $x$ direction, i.e. $h=h(x).$ Consequently, during wave evolution on the slope, $k_{x}$ increases while the transverse component $k_{y}$ remains constant. As a result, the principal wave direction defined by $\tan \theta ={\rm d}y/{\rm d}x$ progressively decreases, aligning more closely with the gradient of bottom topography. This behaviour is consistent with the physical process of wave refraction observed in nearshore regions. This process of composite wavevector evolution agrees with both the mathematical formulation and the physical principles underlying wave refraction in coastal areas. In figure 1(c), we present the evolution process of $\tan \theta$ from wave trains starting with different $\theta _{0}$ . The bottom topography setting (red line) is the same as Lyu et al. (Reference Lyu, Mori and Kashima2023), where the dimensionless depth $kh$ linearly decreases with slope $\gamma =0.05$ from $kh=5$ to $kh=1.2$ , then decreases with a decaying slope $\gamma '=0.05(h/1.55)^{40}$ from $kh=1.2$ to $kh=1.1$ . When the $kh$ becomes smaller than 2.5, the decrease of $\tan \theta$ becomes obvious and gradually reaches a steady state in the shallow-water region, suggesting that the wave principal direction tends to align more closely with the gradient of depth $h$ .

It should be noted that, due to the approximation ${\rm d}x'={\rm d}x$ , the total wavenumber $k$ may deviate from its theoretical value. Therefore, it is essential to constrain this difference by appropriately selecting the value of $\theta _{0}$ within an acceptable range. According to Snell’s law, the ratio between phase speed $c_{p}$ and $\sin \theta$ remains constant throughout the refraction process:

(3.4) \begin{align} \frac{c_{p1}}{\sin \theta _{1}}=\frac{c_{p2}}{\sin \theta _{2}}, c_{p}=\frac{\omega _{0}}{k}. \end{align}

Since $\omega _{0}$ keeps constant in the propagation, this relationship at each numerical step can be expressed as

(3.5) \begin{align} k\left(x\right)\sin \left[\arctan \left(\frac{k_{y}}{k_{x}}\right)\right]=C_{s}\left(x\right)=C_{s}\left(x_{0}\right), \end{align}

where $C_{s}$ should be a constant. Deviation caused by ${\rm d}x'={\rm d}x$ also results in a deviation between $C_{s}(x)$ and $C_{s}(x_{0})$ as the calculation proceeds, and we estimate them by examining the $C_{s}(x)/C_{s}(x_{0})$ in figure 1(d) with the bottom topography the same as already stated. As the refraction occurs, $C_{s}(x)/C_{s}(x_{0})$ becomes smaller with the decrease of $kh$ and the increase of initial angle $\theta _{0}$ makes this process more obvious. Numerical results indicate that when the $\tan \theta _{0}$ is smaller than 0.2, the error introduced by the approximation ${\rm d}x'={\rm d}x$ remains under 1 %; when $\tan \theta _{0}$ is smaller than 0.1, the error remains under 0.2 %, suggesting acceptable accuracy.

Based on the previous assumptions, we compute the spatial evolution of the wave envelope $\overline{A}$ along the wave propagation direction using (2.4), with ${\rm d}x$ as the spatial step size. At each node, the envelope is resolved in the $y{-}t$ dimensions. The final surface elevation at step $x_{i}$ is given by

(3.6) \begin{align} \eta \left(x_{i},y,t\right) & =\varepsilon {Re}\left[\frac{1}{2}\overline{A}\left(x_{i},y,t\right)E\right] \nonumber \\ & \quad +\varepsilon ^{2}{Re}\left[\frac{k_{i}\cosh k_{i}h}{8\sinh ^{3}k_{i}h}\big(2\cosh ^{2}k_{i}h+1\big)\overline{A}^{2}\left(x_{i},y,t\right)E^{2}\right]\!, \end{align}

where $E= \exp \{\textrm{i} [ \int ^{{x_{i}}}k_{{x_{i}}} (x^{*}){\rm d}x^{*}+k_{y}y-\omega _{0}t \}.$ In figure 2, we give the transient surface elevation on $x{-}y$ dimensions at $t=40T_{0}$ starting from different $\theta _{0}$ with the same randomly generated $\overline{A}$ , where $\sigma _{\theta }=0.3$ , BFI = 0.5 and the bottom is flat at $kh=6$ . More details about numerical parameters will be provided in § 4. We only give the surface elevation over a flat bottom, because $\eta (x,y)$ at $t=40T_{0}$ is discretely integrated by (3.6) from $y{-}t$ dimensions and the variation of $c_{p}$ on an uneven bottom will increase the difficulty in constructing a transient surface, which is a limitation of pseudo spectral method. In figure 2, Since $\eta$ are computed from the same wave envelope $\overline{A}$ , the wave distribution and propagation trends are largely similar, with the primary difference being the degree of inclination of crest lines. With the increase of $\theta _{0}$ , the crest lines become progressively more inclined, indicating a growing angle between the principal wave propagation direction and the $x$ axis.

Figure 2.

Surface elevation $\eta$ at $t$ = 40 $T_{0}$ from different $\tan \theta _{0}$ over a flat bottom, where $\sigma _{\theta }=0.3,kh=6$ .

To ensure higher precision and minimise directional distortion, we restrict the $\tan \theta _{0}$ within 0.1 throughout this study. In addition to deviation due to simplifications, wave breaking is not considered in the present model, nor are the effects of refraction on wave height. To mitigate the influence of these unmodeled processes, the minimum water depth is set to correspond to $kh=1$ . Following Lyu et al. (Reference Lyu, Mori and Kashima2023), we can perform Monte Carlo simulations based on the initial random phase $\psi$ to statistically analyse the influence of $\theta _{0}$ on the two-dimensional directional wavefield propagation over uneven bottom topography, with particular focus on nonlinear evolution and wave deformation processes.

3.2. Extension to crossing waves

We next consider a more specialised case. Due to the simplifications introduced in the numerical treatment of refraction and bathymetric effect, the small oblique incident angle $\theta _{0}$ affects only the reconstruction of the free surface elevation based on (3.3) and (3.6). It does not directly contribute to the spatiotemporal evolution of the wave envelope $\overline{A}$ . It enables us to investigate special cases in a single dNLS equation including multiple wave systems with different $\theta _{0}$ . This situation does exist in a realistic ocean: wave groups generated under similar initial conditions can exhibit slight directional deviations due to a combination of external factors such as wind forcing, diffraction or bottom-induced reflection. We consider the crossing of multiple wave trains that evolve from the same conditions under the same boundary conditions and differ in their propagation direction.

Under the assumption of a small difference and with random phase shifts applied to suppress phase-induced interference on modulational instability, we neglect the coupling effect between different wave trains and describe this process in a single dNLS equation. If we denote $k_{ym}$ as the transverse wavenumber component of the wave train $ m $ , then the total surface elevation $\eta$ can be expressed as

(3.7) \begin{align} k_{ym}=k_{0}\tan \theta _{0m}, \quad \overline{A}E=\sum \overline{A}_{m} E_{m} , \\[-28pt] \nonumber \end{align}

(3.8) \begin{align} \eta =\varepsilon {Re}\left[\frac{1}{2}\sum \overline{A}_{m}E_{m}\right]+\varepsilon ^{2}{Re}\left[\frac{k_{i}\cosh k_{i}h}{8\sinh ^{3}k_{i}h}\big(2\cosh ^{2}k_{i}h+1\big)\left(\sum _{m}\overline{A}_{m}E_{m}\right)^{2}\right]\!, \\[0pt] \nonumber \end{align}

where $E_{m}= \exp \{\textrm{i}[\int ^{{x_{i}}}k_{{x_{i}}}(x^{*}){\rm d}x^{*}+k_{ym}y-\omega _{0}t +{\psi }_{m}^{\prime}]\}$ , ${\psi }_{m}^{\prime}$ provides the phase difference for each sample, which remains constant for all frequencies and directions in a given sample, but is independently and randomly uniformly distributed over $[0,2\pi )$ across different Monte Carlo realisations.

Due to the complexity of the analysis of a multi-system, this study focuses on the simplest and most commonly observed case – crossing between two symmetric wave trains with opposite principal directions. Specifically, we consider two wave trains (denoted as Train 1 and Train 2) that differ only in their initial incident angles $\theta _{0}$ (i.e. $k_{ym}$ ), and $\overline{A}_{1}=\overline{A}_{2}$ . Their initial crossing angle is $\theta _{12}=| \theta _{01}| +| \theta _{02}| =2| \theta _{01}|$ and we define $t_{c}=\tan \theta _{01}=-\tan \theta _{02}$ as the tangent-based parametrisation for convenience.

In a similar way to figure 2, we give the transient surface elevation on $x{-}y$ dimensions at $t=40T_{0}$ starting from different $t_{c}$ with the same randomly generated $\overline{A}$ , where $\sigma _{\theta }=0.3$ , BFI = 0.5 and the bottom is flat at $kh=6$ in figure 3. Since we use two wave train systems with perfectly symmetric incident directions, their transverse propagation components effectively cancel out, resulting in the principal direction being fully aligned with the $x$ axis. However, the surface elevations corresponding to different $t_{c}$ values exhibit significant differences in transverse inhomogeneity. As $t_{c}$ increases, this inhomogeneity decreases, while a series of more densely spaced regions with vanishing surface elevation, aligned with the primary propagation direction, become increasingly apparent. This behaviour arises from the symmetric construction with identical wave amplitudes $\overline{A}$ , for which (3.8) has periodic zero solutions satisfying . At these locations, the surface elevation $\eta =0$ identically as a result of transverse destructive interference between the two wave trains in idealised model.

Figure 3.

Surface elevation $\eta$ of crossing waves at $t$ = 40 $T_{0}$ from different $t_{c}$ over a flat bottom, where $\sigma _{\theta }=0.3,kh=6$ .

Changes brought about by $t_{c}$ shows the interference behaviour of crossing waves: as the crossing angle increases, energy dispersion within the wavefield becomes more pronounced, thereby reducing transverse inhomogeneity. However, the role of transverse wave group structures in the overall wavefield evolution and the likelihood of extreme events remains unclear, especially on a sloping bottom with varying $k_{i}$ and refraction. Moreover, a single random realisation is not sufficient to draw statistically robust conclusions.

Therefore, we proceed with Monte Carlo simulations to systematically investigate the evolution of crossing wavefields and to discuss the influence of $t_{c}$ on their statistical properties. We set ${\psi }_{m}^{`}$ uniformly distributed at $[0,2\pi )$ for each member. Different to the random phase $\psi$ equation (3.3), which is randomly assigned to each frequency and directional component, ${\psi }_{m}^{`}$ in (3.8) remains constant across all components within each individual wave train. Since the ensemble averages satisfy , the contribution of the second-order term in (3.8) in Monte Carlo results is dominated by the self-interaction component ${(\overline{A}_{m}}{E_{m}})^{2}$ . Considering the constraint of the small oblique angle imposed by the fundamental assumptions of the model, $t_{c}$ in our simulations varies from 0.01 to 0.07.

4.1. Oblique incident waves

Based on the previous assumptions and methodology already described, numerical experiments are conducted over uneven bottom topography. The simulations begin at the $x_{0}$ location on the $y{-}t$ plane, with $N_{t}=1000$ time samples and $N_{y}=900$ spatial samples in the $y$ direction. In subsequent computations, all time-related variables are non-dimensionalised by wave period $T$ and spatial variables by initial wavelength $L_{0}$ as was done in Lyu et al. (Reference Lyu, Mori and Kashima2023). The sampling intervals are set to ${\rm d}t=0.1$ , ${\rm d}y=1$ . The computational step size in the $x$ direction is ${\rm d}x=2\times 10^{-5}$ , with output recorded at every ${\unicode[Arial]{x0394}} x={\rm d}y$ . The wave period $T$ is held constant at $T=2.5$ throughout the simulation from medium water to shallow water and the initial water depth is set to $h=5.7$ . The corresponding initial wavelength, computed from the dispersion relation, is $L_{0}=9.723$ with $k_{0}=0.646$ . Initial wave steepness $\varepsilon =k_{0}a=0.1$ . Initial spectral bandwidth is fixed at $\sigma _{s}=0.285$ , then $\textrm{BFI}=0.5$ . Directional dispersion varies at $\sigma _{\theta }=0.3,0.4,0.5$ . Initial incident angle $\theta _{0}$ varies from $\tan \theta _{0}=0$ to $\tan \theta _{0}=0.1$ . Bottom topography is set as figure 4, where wave trains propagate from $kh=3.7$ to $kh=1.03$ with a constant slope $\gamma _{{s}}=0.05$ , and we selected some representative sections of $kh$ . Contributions from the initial BFI, $\sigma _{\theta }$ and bottom topography have been discussed by Lyu et al. (Reference Lyu, Mori and Kashima2021, Reference Lyu, Mori and Kashima2023), and this study concentrates on how the initial incident angle $\theta _{0}$ affects the evolution of wave trains. The role of directional dispersion in wave refraction remains unclear; therefore, we are also interested in examining the results under varying $\sigma _{\theta }$ and $kh$ conditions.

Figure 4.

Bottom topography from deep to shallow water with a constant slope $\gamma _{{s}}=0.05$ .

The primary output of the simulation is the surface elevation $\eta (y,t)$ at each spatial location $x_{i}$ , computed based on (3.6). By applying a Fourier transform, we obtain the power spectrum $P$ , which characterises the distribution of wave energy in the frequency and wavenumber domains. In this study, we are particularly interested in the lateral spatial structure of the wavefield, which reflects the effect from $\sigma _{\theta }$ . Therefore, instead of the conventional frequency-based spectra, we focus on the power spectrum $P$ as a function of the transverse wavenumber $k_{y}$ , which reveals the contribution from different directional components. The time length is $({N_{t}\,{\rm d}t})/{T}=$ 40 times wave periods, which is sufficient to ensure temporal ergodicity. The transverse spatial length is also validated through correlation length analysis to ensure adequate sampling coverage. For further details, refer to Lyu et al. (Reference Lyu, Mori and Kashima2023).

In figure 5, we present the transverse wavenumber spectra $P(k_{y})$ at different water depth $kh$ , obtained from a single random sample by performing a Fourier transform over all temporal points and subsequently averaged. In figure 5(a), we present the time-averaged spectra $P(k_{y})$ corresponding to $\tan \theta _{0}$ values ranging from 0 to 0.1, with an interval of 0.02. Despite the uncertainties introduced by random initial phases, the influence of $\tan \theta _{0}$ on the wave spectra remains clearly observable. Since wave train propagation is primarily along the $x$ axis, the transverse components in the $y$ direction come from directional dispersion, which is assumed to follow a Gaussian distribution centred on the principal propagation direction as (3.3). It is reflected in the $P(k_{y})$ spectra, where the energy is concentrated near $k_{y}=0\sim 0.1k_{0}$ and decays with increasing $k_{y}$ . As $\tan \theta _{0}$ increases, indicating a greater obliquity of the incident wave direction, the transverse wave component should become more significant theoretically. This trend is evident in the figure: the peak magnitude of $P(k_{y})$ gradually decreases, while the wave energy shifts towards higher $k_{y}$ values. It implies a broader directional spread and a higher proportion of wave energy propagating in the lateral direction.

Figure 5.

Time-averaged power spectra $P(k_{y})$ of surface elevation $\eta (y)$ with different $\tan \theta _{0},kh,\sigma _{\theta }$ over a slope, derived from a specific sample: (a) $kh=3,\sigma _{\theta }=0.3$ ; (b) $kh=3$ ; (c) $kh=2$ ; (d) $kh=1.36$ .

Figure 5. Long description

Panel A: A line graph shows the time-averaged power spectra of surface elevation with different values of tan theta 0 and a fixed value of sigma theta. The x-axis represents k y divided by k 0, and the y-axis represents P. The graph includes multiple lines representing different values of tan theta 0, ranging from 0 to 0.1. Panel B: A line graph shows the time-averaged power spectra of surface elevation with different values of tan theta 0 and sigma theta. The x-axis represents k y divided by k 0, and the y-axis represents P. The graph includes multiple lines representing different combinations of tan theta 0 and sigma theta. Panel C: A line graph shows the time-averaged power spectra of surface elevation with different values of tan theta 0 and sigma theta. The x-axis represents k y divided by k 0, and the y-axis represents P. The graph includes multiple lines representing different combinations of tan theta 0 and sigma theta. Panel D: A line graph shows the time-averaged power spectra of surface elevation with different values of tan theta 0 and sigma theta. The x-axis represents k y divided by k 0, and the y-axis represents P. The graph includes multiple lines representing different combinations of tan theta 0 and sigma theta.

Furthermore, figures 5(b)–5(d) present the time-averaged $P(k_{y})$ in the same form under various water depths $kh=3,2,1.36$ combined with different $\sigma _{\theta }=0.3,0.4,0.5$ . Based on previous dNLS simulations, the depths in this range are sufficiently representative to capture contribution from topography change and $kh=1.36$ is a key threshold for the Benjamin–Feir instability. By systematically comparing different sets of results, we can quantitatively assess the individual effects of the three variables: water depth $kh$ , directional dispersion $\sigma _{\theta }$ and initial wave direction $\theta _{0}$ . First, we observe that as $kh$ decreases, the distribution of $P(k_{y})$ remains nearly unchanged, indicating that the transverse transfer of wave energy is very small, which is one of the underlying reasons for wave refraction. A slight overall reduction in spectral amplitude is also noted, which can be attributed to that our computation of $P(k_{y})$ does not account for the increase in wave energy density associated with the decrease in $c_{\kern-1pt g}$ during shoaling. In reality, assuming no dissipation or breaking, the total wave energy should be conserved and the energy density increases as wave trains propagate from deep to shallow water. Next, we note that the effect of $\sigma _{\theta }$ is significant: as $\sigma _{\theta }$ increases, directional dispersion becomes stronger, leading to a reduction in the spectral peak of $P(k_{y})$ at lower $k_{y}$ (approximately $k_{y}\lt 0.25k_{0}$ ), while the energy at higher $k_{y}$ values increases. This indicates an enhanced contribution of transverse wave components and a broader directional spreading of wave energy. Moreover, an increase in $\sigma _{\theta }$ also reduces the differences induced by varying $\tan \theta _{0}$ , as the enhanced directional dispersion overwhelms the minor effects introduced by small-angle oblique incidence.

Figure 6.

Statistical parameters of wave surface morphology with different $\tan \theta _{0},kh$ over a slope derived from Monte Carlo simulation: (a) averaged maximum wave height at $\sigma _{\theta }=0.3$ ; (b) averaged run length defined by mean wave height.

To further investigate how differences in wave spectrum are reflected in the spatial structure of wave trains, we focused on two key parameters. The first is the statistical distribution of the maximum wave height, which helps evaluate the extreme events. Since our primary interest lies in the influence of the incident wave angle, we evaluated, for each sample, the transverse wave height sequences at each $(x,t)$ node and extracted the maximum wave height $H_{max }$ by the root-mean-square (r.m.s.) of surface elevation, $\eta _{\textrm{rms}}$ , from each sample. The second parameter is the wave run length, which provides a measure of the persistence or clustering of large waves (Medina & Hudspeth Reference Medina and Hudspeth1990). Specifically, we analysed the transverse wave height sequences and computed the lengths of consecutive wave heights that exceed a given threshold, following the run length methodology. This length is expressed as $R_{{m}}$ , the mean value of numbers of consecutive waves exceeding the mean wave height.

To reduce the randomness introduced by the initial phase in (3.3) and more clearly isolate the effects of different parameters, we conducted large-scale Monte Carlo simulations. For each set of initial conditions, we generated 600 wave evolution samples, in which only the initial phase distribution varied. In figure 6(a), we present the Monte Carlo-averaged values of $H_{max }/\eta _{\textrm{rms}}$ at different water depth $kh$ , as a function of $\tan \theta _{0}$ . The results indicate a clear reduction in maximum wave height as the water becomes shallower, with minimal sensitivity to the initial incidence angle of the waves. The observed dependence of wave height on water depth is consistent with previous studies (Lyu et al. Reference Lyu, Mori and Kashima2021, Reference Lyu, Mori and Kashima2023, etc.). The absence of a clear $\theta _{0}$ dependence of $H_{max }$ may be attributed to the relatively small magnitude of the effect, which is further limited by the basic assumptions of the dNLS model that restrict strongly oblique wave incidence. The transverse component primarily contributes through the first-order linear term in (3.6), while the second-order term is sometimes negligible and may be overshadowed by directional dispersion effects. In figure 6(b), we present a more detailed depiction of how $R_{{m}}$ varies with $\tan \theta _{0}$ under different conditions. The results reveal a weak increasing trend in $R_{{m}}$ with increasing $\tan \theta _{0}$ , but this trend becomes less significant as $\sigma _{\theta }$ increases. This suggests that a larger incident angle may slightly generate more persistent large-wave sequences in the transverse direction, indicating the increase of wave groupiness, and risk of cumulative loading on coastal and offshore structures. However, under strong directional dispersion, this effect becomes negligible, as the wave energy is more diffusely distributed across different directions, leading to reduced coherence in wave group formation and shorter large-wave sequences.

Figure 7.

Nonlinear interactions at $\sigma _{\theta }=0.3$ for different $\tan \theta _{0},kh$ over a slope derived from Monte Carlo simulation: (a) averaged kurtosis; (b) averaged skewness; (c) standard deviation of kurtosis; (d) standard deviation of skewness.

In previous studies on the four-wave interactions and the freak wave, the fourth-order moment and the third-order moment of the surface elevation, kurtosis $\mu _{4}$ and skewness $\mu _{3}$ of the surface elevation, are useful for analysing contributions from different sources of nonlinear effect to the occurrence probability of extreme wave heights (Janssen Reference Janssen2003; Mori & Janssen Reference Mori and Janssen2006; Lyu et al. Reference Lyu, Mori and Kashima2021, Reference Lyu, Mori and Kashima2023). To gain further insight into the role of $\theta _{0}$ in nonlinear wave evolution and refraction, we computed $\mu _{4}$ and $\mu _{3}$ of the surface elevation $\eta$ in time series at fixed spatial points, then took their averaged value in the transverse direction. In figures 7(a) and 7(b), the ensemble-averaged results from 600 Monte Carlo samples are plotted along the direction of bottom variation with different $\tan \theta _{0}$ at $\sigma _{\theta }=0.3$ . Over the varying bathymetry, the evolutions of $\mu _{4}$ and $\mu _{3}$ are in good agreement with the previously reported numerical results (Zeng & Trulsen Reference Zeng and Trulsen2012; Lyu et al. Reference Lyu, Mori and Kashima2021) and experimental observations (Kashima & Mori Reference Kashima and Mori2019; Kimmoun et al. Reference Kimmoun, Hsu, Hoffmann and Chabchoub2021). However, when attention is directed to the effect of oblique angle, both $\mu _{4}$ and $\mu _{3}$ remain essentially unchanged for different values of $\tan \theta _{0}$ , suggesting that the influence of the incident angle is mainly restricted to linear effects. Small variations in $\tan \theta _{0}$ do not contribute meaningfully to either three-wave or four-wave interactions. To capture statistical features that may be obscured by transverse averaging, we computed their standard deviation $s(\mu _{4})$ and $s(\mu _{3})$ across the transverse direction. Their ensemble-averaged values over the Monte Carlo simulations are presented in figures 7(c) and 7(d), plotted as functions of the $\tan \theta _{0}$ at $\sigma _{\theta }=0.3$ . Consistently, both $s(\mu _{4})$ and $s(\mu _{3})$ show clear dependence on the water depth $kh$ , but no significant trend with respect to $\tan \theta _{0}$ .

In summary, small initial values of $\tan \theta _{0}$ have a negligible effect on the nonlinear evolution of wave trains, and consequently on the contribution to the scale and the occurrence probability of freak waves. However, when the directional dispersion is relatively weak, an increase in $\tan \theta _{0}$ can influence the transverse distribution of wave energy, thereby enhancing wave groupiness in the transverse direction.

4.2. Crossing waves

For the simulation of crossing waves described by (3.8), we adopt the same initial parameters and bottom boundary conditions as in oblique incident waves. By definition, $t_{c}$ and $\tan \theta _{0}$ are conceptually similar, as both represent the transverse propagation component in the initial wavefield. To better understand their similarities and differences, we perform a quantitative comparison of their effects on wavefield evolution through similar simulation procedures. As in § 4.1, we investigate the evolution of the wavefield and the distribution of wave energy through the Monte Carlo simulation, focusing sequentially on the wave spectrum, surface elevation patterns and nonlinear interactions. However, unlike the previous result in a single system, we find that the influence of $t_{c}$ on the wave spectra is much smaller and cannot be clearly identified from individual random samples. Therefore, all results presented in this section are ensemble-averaged over 600 Monte Carlo simulations.

In figure 8, we present the transverse wavenumber spectra $P(k_{y})$ under various groups of water depth $kh$ , directional spreading $\sigma _{\theta }$ and crossing angle $t_{c}$ . It can be observed that the spectral distribution of the crossing wave systems is basically consistent with that of the single-wave system as shown in figure 5. Moreover, the effects of varying $kh$ and $\sigma _{\theta }$ on the energy distribution keep giving a similar process. As $t_{c}$ increases, the peak region of $P(k_{y})$ gradually decreases in magnitude, with a corresponding shift of energy towards larger $k_{y}$ values. It indicates enhanced transverse dispersion of wave energy. The effect is similar to that observed for increasing $\tan \theta _{0}$ in figure 5, although the magnitude of the impact is notably smaller. Both effects are modulated by the directional spreading parameter $\sigma _{\theta }$ : as $\sigma _{\theta }$ increases, directional dispersion becomes stronger, then the influence of $t_{c}$ and $\theta _{0}$ becomes increasingly less significant. Up to this point, we conclude that the degree of lateral interference in symmetric crossing waves shown in figure 3 has only a minor effect on the transverse energy distribution.

Figure 8.

Time-averaged power spectra $P(k_{y})$ of surface elevation $\eta (y)$ with different $t_{c},kh,\sigma _{\theta }$ over a slope, derived from Monte Carlo simulation: (a) $kh=3.7$ ; (b) $kh=3$ ; (c) $kh=2$ ; (d) $kh=1.36$ .

Figure 8. Long description

Panel A: A line graph depicts time-averaged power spectra of surface elevation with kh equal to 3.7. The x-axis is labeled k_y/k_0 ranging from 0 to 0.8, and the y-axis is labeled P ranging from 0 to 0.45. The graph includes multiple lines representing different values of t_c and sigma_theta. Panel B: A line graph depicts time-averaged power spectra of surface elevation with kh equal to 3. The x-axis is labeled k_y/k_0 ranging from 0 to 0.8, and the y-axis is labeled P ranging from 0 to 0.45. The graph includes multiple lines representing different values of t_c and sigma_theta. Panel C: A line graph depicts time-averaged power spectra of surface elevation with kh equal to 2. The x-axis is labeled k_y/k_0 ranging from 0 to 0.8, and the y-axis is labeled P ranging from 0 to 0.45. The graph includes multiple lines representing different values of t_c and sigma_theta. Panel D: A line graph depicts time-averaged power spectra of surface elevation with kh equal to 1.36. The x-axis is labeled k_y/k_0 ranging from 0 to 0.8, and the y-axis is labeled P ranging from 0 to 0.45. The graph includes multiple lines representing different values of t_c and sigma_theta.

As a more direct comparison, we examine the maximum wave height $H_{max }/\eta _{\textrm{rms}}$ and the transverse wave run length $R_{{m}}$ in figure 9. Compared with figure 6, $H_{max }/\eta _{\textrm{rms}}$ and $R_{{m}}$ in the multi-crossing system give higher values than the single system, which mainly comes from the linear superposition effect. This finding is consistent with previous laboratory and numerical results (Onorato et al. Reference Onorato, Osborne and Serio2006; Toffoli et al. Reference Toffoli, Bitner-Gregersen, Osborne, Serio, Monbaliu and Onorato2011; McAllister et al. Reference McAllister, Draycott, Adcock, Taylor and van den Bremer2019) at angles of approximately 40°–60°, which have shown that even relatively small crossing angles can enhance the maximum wave height. In figure 9(a), similar to the case of $\tan \theta _{0}$ in the single-wave system in figure 6(a), $t_{c}$ shows no significant effect on $H_{max }$ under any of the tested conditions. However, an unexpected result emerges in figure 9(b): as $t_{c}$ increases, $R_{{m}}$ decreases significantly, and the differences previously induced by water depth $kh$ and directional spreading $\sigma _{\theta }$ progressively diminish. When $t_{c}=0.07$ , the influence of both $kh$ and $\sigma _{\theta }$ on $R_{{m}}$ becomes negligible. This result highlights an important distinction between the multi-crossing system and the single wave system. Although both $t_{c}$ and $\theta _{0}$ provide the initial transverse propagation components of the wavefield, in crossing wave systems, the symmetric superposition of two wave trains with varying directional separation leads to transverse interference effects in a different way. As $t_{c}$ becomes larger, this interference suppresses coherent group formation on the $y$ axis, resulting in shorter transverse wave sequences and lower values of $R_{{m}}$ . It is consistent with figure 3 and provides more quantitative results: smaller $t_{c}$ leads to stronger wave groupiness and more persistent sequences of large waves, manifesting as increased transverse inhomogeneity in the surface elevation, which are similar to the larger crossing angle’s result in experiments by Waseda et al. (Reference Waseda, Kinoshita and Tamura2009) and Luxmoore et al. (Reference Luxmoore, Ilic and Mori2019), and the numerical result of Gramstad & Trulsen (Reference Gramstad and Trulsen2010). The diminishing influence of $kh$ and $\sigma _{\theta }$ further indicates that the spatial coherence of wave groups is largely dominated by the directional incoherence introduced by crossing. This is similar to the effect of increasing $\sigma _{\theta }$ , where enhanced dispersion reduces the apparent influence of other parameters.

Figure 9.

Statistical parameters of wave surface morphology with different $t_{c},kh,\sigma _{\theta }$ over a slope derived from Monte Carlo simulation: (a) averaged maximum wave height; (b) averaged run length defined by mean wave height.

With the establishment that $t_{c}$ makes a more pronounced influence on wave surface patterns compared with $\theta _{0}$ , we now aim to investigate in more detail how this influence is reflected in different forms of nonlinear interactions. In figure 10, we present kurtosis $\mu _{4}$ , skewness $\mu _{3}$ , and their standard deviations $s(\mu _{4})$ and $s(\mu _{3})$ with respect to $t_{c}$ and $x_{0}/L_{0}$ at different conditions from Monte Carlo simulation. As shown in figures 10(a) and 10(b), both $\mu _{4}$ and $\mu _{3}$ are primarily influenced by water depth $kh$ and directional spreading $\sigma _{\theta }$ , while the effect of $t_{c}$ remains negligible. Although $\mu _{3}$ exhibits some degree of fluctuation with increasing $t_{c}$ , the variation is not sufficiently systematic to attribute it to the influence of $t_{c}$ . However, we note that the values of $\mu _{4}$ are significantly higher than those presented in figure 7 and reported in previous studies (Lyu et al. Reference Lyu, Mori and Kashima2021, Reference Lyu, Mori and Kashima2023), indicating that the crossing wave system induces stronger four-wave interactions compared with the single-wave system. The values of $\mu _{3}$ are consistently negative, which may be attributed to the use of symmetrically propagating wave trains in our set-up.

Figure 10.

Nonlinear interactions with different $t_{c},kh,\sigma _{\theta }$ over a slope derived from Monte Carlo simulation: (a) averaged kurtosis; (b) averaged skewness; (c) standard deviation of kurtosis for $t_{c}$ ; (d) standard deviation of skewness for $t_{c}$ ; (e) standard deviation of kurtosis for $x/L_{0}$ ; (f) standard deviation of skewness for $x/L_{0}$ .

Figure 10. Long description

Panel A: A scatter plot showing averaged kurtosis with the x-axis labeled as t_c and the y-axis labeled as μ4. Different symbols and colors represent various kh and σg values. Panel B: A scatter plot showing averaged skewness with the x-axis labeled as t_c and the y-axis labeled as μ3. Different symbols and colors represent various kh and σg values. Panel C: A scatter plot showing the standard deviation of kurtosis for t_c with the x-axis labeled as t_c and the y-axis labeled as s(μ4). Different symbols and colors represent various kh and σg values. Panel D: A scatter plot showing the standard deviation of skewness for t_c with the x-axis labeled as t_c and the y-axis labeled as s(μ3). Different symbols and colors represent various kh and σg values. Panel E: A line graph showing the standard deviation of kurtosis for x/L0 with the x-axis labeled as x/L0 and the y-axis labeled as s(μ4). Different line styles and colors represent various t_c and σg values. Panel F: A line graph showing the standard deviation of skewness for x/L0 with the x-axis labeled as x/L0 and the y-axis labeled as s(μ3). Different line styles and colors represent various t_c and σg values.

Compared with the single-wave system, the increase in $\mu _{4}$ induced by wave crossing is consistent with previous findings (e.g. Toffoli et al. Reference Toffoli, Bitner-Gregersen, Osborne, Serio, Monbaliu and Onorato2011; Liu et al. Reference Liu, Waseda, Yao and Zhang2022), and our results further demonstrate that even small crossing angles are sufficient to cause a notable increase in kurtosis. However, when it comes to the effect of $t_{c}$ on $\mu _{4}$ , previous studies reveal that this effect is not universally consistent. Several works (Toffoli et al. Reference Toffoli, Bitner-Gregersen, Osborne, Serio, Monbaliu and Onorato2011; Sabatino & Serio Reference Sabatino and Serio2015; Brennan, Dudley & Dias Reference Brennan, Dudley and Dias2018; Gramstad et al. Reference Gramstad, Bitner-Gregersen, Trulsen and Borge2018) reported that $\mu _{4}$ first increased then decreased with the increase of crossing angle typically from 40° to 60°. In contrast, Luxmoore et al. (Reference Luxmoore, Ilic and Mori2019) did not find any significant relationship between $\mu _{4}$ and $t_{c}$ in their experiments, which may be attributed to the difference between long-crested and short-crested crossing seas (Liu et al. Reference Liu, Waseda, Yao and Zhang2022). As a supplement at small crossing angles, the ensemble-averaged $\mu _{4}$ shows limited sensitivity to $t_{c}$ under the Monte Carlo simulation conditions used in this study.

Upon further examination of the standard deviation $s(\mu _{4})$ and $(\mu _{3})$ , we observe that $t_{c}$ , in contrast to $\theta _{0}$ , gives a significant effect on their evolution. In figures 10(c) and 10(d), $s(\mu _{4})$ and $s(\mu _{3})$ pronouncedly decrease as the increase of $t_{c}$ . It indicates that in crossing wave systems, greater differences in propagation direction lead to weaker four- and three-wave interactions in the resulting wavefield. Although the ensemble-averaged values of $\mu _{4}$ and $\mu _{3}$ remain largely unchanged, the statistical homogeneity of the wavefield increases with growing $t_{c}$ . This is attributed to the enhanced transverse interference between the symmetrically propagating wave trains, which strengthens the overall dispersion and suppresses local energy focusing.

In figures 10(e) and 10(f), we present the evolution of $s(\mu _{4})$ and $s(\mu _{3})$ with respect to $x_{0}/L_{0}$ at different $t_{c}$ and $\sigma _{\theta }$ . The effects of $t_{c}$ and $\sigma _{\theta }$ have something in common, as they both reduce surface inhomogeneity and instability by enhancing directional dispersion. However, their differences are also significant: while $\sigma _{\theta }$ affects both the averaged value and standard deviation of $\mu _{4}$ , the impact of $t_{c}$ is primarily confined to the standard deviations. It suggests that the role of $t_{c}$ is more closely associated with fine-scale instability that may not significantly alter the ensemble mean value, but contribute to the structural variability of the wavefields. Moreover, from a quantitative perspective, the influence of $t_{c}$ on the standard deviations is more pronounced than that of $\sigma _{\theta }$ , especially in $s(\mu _{3})$ . In addition, figures 10(e) and 10(f) present the results for the single-wave system under normal incidence ( $\tan \theta _{0}=0$ ) with $\sigma _{\theta }=0.3$ , indicated by green lines. Here, $s(\mu _{4})$ in the crossing wave system is significantly higher than in the single system across all tested conditions, suggesting stronger nonlinear interactions and a higher probability of potential extreme events. The results for $s(\mu _{3})$ show a similar trend but are largely insensitive to variations in $\sigma _{\theta }$ . When $t_{c}=0.07$ , the $s(\mu _{3})$ value in the multi-crossing system decreases to a similar level to that in the single system.

In the investigation of freak wave generation mechanisms, the model is primarily based on the surface established from the linear superposition of the amplitude, which contributes to most of the observed increases in maximum wave height and kurtosis. However, the variations in the standard deviations of kurtosis and skewness also reflect the underlying contributions of second- and third-order nonlinear interactions during wave evolution. It should be noted that our analysis is based on two symmetric wave trains, which may limit the generality of some symmetry-related results, such as those involving skewness. To get a more comprehensive understanding of the evolution of nonlinear crossing wave systems, further simulations under broader conditions, along with the development of corresponding theoretical frameworks, are needed.

Overall, our results demonstrate that even small angular separations in propagation direction can have a non-negligible influence on the wavefield structure and nonlinear interactions in crossing wave systems. In the context of extreme event assessment, wave crossing with small directional differences may lead to more persistent large-wave sequences and increased non-Gaussianity of surface elevation, thereby enhancing the occurrence probability of freak waves and posing elevated risks for coastal and offshore hazard prevention work.