Wave loads on ocean infrastructure increase as a result of waves passing over abrupt depth transitions

Abrupt changes in water depth are known to lead to abnormal free-surface wave statistics. The present study considers whether this translates into abnormal loads on offshore infrastructure. A fully non-linear numerical model is used which is carefully validated against experiments. The wave kinematics from the numerical model are used as input to a simple wave loading model. We find enhanced overturning moments, an increase of approximately 20%, occur over a distance of a few wavelengths after an abrupt depth transition. We observe similar results for 1:1 and 1:3 slopes. This increase does not occur in linear simulations.


Introduction
Ocean waves are central to the study of ocean engineering and are a key design driver for many marine energy devices. Standard models of wave statistics assume that the underlying freely propagating components are uncorrelated. However, situations may arise where correlations are present between components, and this leads to a greater number of large waves than would otherwise be expected. Such waves can be considered rogue or freak waves (see reviews by Dysthe et al. (2008); Adcock and Taylor (2014)).
One mechanism which can cause an increased number of large waves is where waves propagate over a rapid change in water depth. This mechanism was identified a decade ago ) and has since motivated many researchers (Gramstad et al. 2013;Ducrozet and Gouin 2017;Bolles et al. 2019;Zhang et al. 2019;Zheng et al. 2020; B Tianning Tang tianning.tang@eng.ox.ac.uk 1 et al. (2021a) and Trulsen et al. (2020). The statistics of wave kinematics in long-crested irregular waves propagating over a shoal is studied by Lawrence et al. (2021b) using a Monte Carlo approach. The evolution of kurtosis and skewness of the horizontal velocity is discussed, and the statistical properties are compared to laboratory experiments. Trulsen et al. (2020) conduct the experiments of long-crested irregular water surface waves propagating over a shoal. The location of local maximum and local minimum of skewness and kurtosis are investigated in terms of the surface elevation and the horizontal fluid velocity. In addition, Bateman et al. (2003) investigate the calculation of the water particle kinematics generated by the propagation of surface gravity waves.
Our work seeks to examine the wave kinematics and loads induced on an idealised structure in the vicinity of an abrupt depth transition. We use fully non-linear potential-flow simulations which explicitly solve for the wave kinematics.

Experimental setup
In this paper, the experimental data are taken from Li et al. (2021a, b). In the experiment, only surface elevation is measured. The experimental setup is presented in Fig. 1. A total of 16 gauges were used to measure the waves and the locations of 16 gauges are shown in Table 1. Our study primarily focus on the results at gauges 1, 3, 8, and 11 and the measurements at these 4 gauges are selected for comparison with the numerical simulations. The width of the flume is 0.6 m and the length is 35 m. The top of the slope is 7.5 m away from the wave-maker. The water depth of the deeper water side, h d , is 0.55 m and the water depth of the shallower water side, h s , is 0.2 m. Thus, the height of the slope is fixed as 0.35 m. The slopes are indicated by the dashed lines in Fig. 1, showing a slope of both 1:1 and 1:3.

Numerical setup
OceanWave3D (OW3D) is used to model the waves (see Engsig-Karup et al. (2009)) and the governing equations of OW3D are shown in Appendix A.
Potential-flow equations for surface gravity waves are solved by OW3D in a three-dimensional Eulerian frame , which gives the kinematics information not only at the free-surface boundary but also at all the vertical nodes hence allowing full kinematic reconstruction. Thus, the wave kinematics are extracted directly during the time-marching process without any further assumptions. A classical fourthorder Runge-Kutta time-marching scheme is used. A σcoordinate transformation is applied to map the vertical solutions to a time-invariant grid, which can also be used to solve the fully non-linear potential-flow equations with variable bathymetry.
Wave breaking is implemented using a local smoothing function. This local smoothing function monitors the particle downward acceleration of the free surface ( w in the governing equations of OW3D, see Appendix A) and it removes a small amount of energy locally until the acceleration is below the threshold. In our simulation, the local smoothing function is triggered when the detected acceleration is greater than 0.4g (where g is the gravitational acceleration), which follows previous work (Tang et al. (2021)). Interestingly, in this study, we have experienced unusual breaking events especially around the abrupt depth transitions. These strong breaking waves caused difficulties to get convergence with this simplified breaking model, which helps further stabilise the numerical scheme. The simplified breaking model can only be applied once at each time step. Thus, we found a small number of events occurred where it is presumed, the wave was breaking caused the code to stop. This suggests that breaking at the top of the slope may have significant differences from waves on a flat bathymetry, which is consistent with the recent work of Draycott et al. (2022). When the code stops, a finer resolution in temporal domain is applied to allow the breaking model to trigger over more timesteps allowing the code to compute through the point at which it could not converge with the normal time step. After passing this point, we change to the original resolution until meeting the next crash point.
The model also has a 'linear' version which we use in this paper. This version solves the classic linear form of the water wave equations. In the linear model of OW3D, the surface elevation profile is set to zero in the σ -coordinate transformation and the non-linear terms in the free-surface boundary conditions were also neglected, which leads to generated waves follow the linear wave theory (see Bingham and Zhang (2007) for detailed implementation of this model).
The numerical domain is almost identical to that used in the experimental case, and the relative location of gauges is the same as in the experiment. In our numerical domain, we use gauge 3 as our datum (x = 0 m in the numerical domain), so locations upstream of this are negative. In our simulations, the peak frequency of the input waves is 0.8 Hz. The wavelength is 2.23 m and 1.6 m on deep and shallow side, respectively. The water depth on the deep side and shallow side is 0.55 m and 0.2 m, respectively, which is the same as experimental setup. we consider two different slopes: a slope of 1:1 and one of slope 1:3. The input waves and numerical domain are the same for both cases. The length of the numerical domain, L x , is 54 m (from x = − 21.88 m to x = 32.12 m a length of 24 wavelengths) and the number of nodes on the x-axis, N x , is 3500. The number of nodes on the z-axis, N z , is 20. The timestep, t, is 0.03125 s (40 time steps per wave period) and a total time duration of 3.5 h of wave propagation is simulated for both cases. Waves are generated using a double relaxation zone to help absorb reflected waves, which is located from − 21.88 to − 1.88 m (9 wavelengths). At the end of the domain, waves are absorbed using a pressure boundary condition. Thus, a damping zone is applied at the end of the domain and the length is 10 m (from x = 22.12 m to x =32.12 m and corresponding to 4.5 wavelengths). The top of the slope is located at x = 0 m in the numerical domain for both cases, which is 1.88 m away from the end of the relaxation zone (x = −1.88 m). The same numerical setup is adapted for both linear and non-linear model of OW3D simulations. The steepness of the slope is changed by changing the length of the slope. The discretisation used is based on the study of Barratt et al. (2020) and a grid convergence study on the same geometry used here is presented in Li et al. (2022).

Input waves
In this paper, we primarily consider random wave simulations. In the experiments, the waves are drawn from a JONSWAP spectrum with peak enhancement factor, γ = 3.3, and the JONSWAP spectrum is developed by Hasselmann et al. (1973) in which ω p is the peak frequency (ω p = 2π/T p ) and the parameter s =0.07 for ω < ω p and s = 0.09 for ω > ω p . The input measured significant wave height at the first gauge of the experiments was H s = 0.045 m giving a significant steepness of 1 2 H s k p = 0.064 (the peak wave number, k p , is 2.81m −1 on the deep side). The peak frequency of the input waves is 0.8 Hz giving a non-dimensional water depth on the deep side of k p h d = 1.55 and, assuming linear shoaling, a non-dimensional depth on the shallow side of k p h s = 0.79. In our simulations, the experimental data at gauge 1 are used as input for the OW3D simulations. We ignore the first 300 s of the experiment to allow the waves to reach a steady state.

Validation against experiments
Extensive verification and validation, in addition to a convergence study, of the model used in the present paper was undertaken in Li et al. (2022) and is not repeated here. This previous work focused on deterministic wave groups rather than random waves as are studied here. As such, we here focus on the aspects of the modelling that apply to random waves but not to deterministic groups. These aspects are principally the wave generation and wave breaking (as the groups in the previous work were designed not to break). We note that the wave generation in the numerical model is different from that in the experimental model. In addition, this potentialflow solver can never fully capture the non-potential-flow physics induced by wave breaking (e.g., Eeltink et al. (2022)) and vortex shedding due to the sharp edge at the top of the slope. The experiments and numerics will also differ in that some dissipation from the side walls will occur in the experiments and also differ in wave absorption with more complete absorption being expected in the numerical scheme. Another expected source of inaccuracy is that we use experimental data at gauge 1 as input to the numerical model. However, these data include waves reflected off the depth transition which are travelling in the opposite direction. In our validation, we present results for the 1:1 slope as similar results are observed for the 1:3 case, and intuitively, the latter case is expected to be marginally easier to simulate as non-linear physics is likely to be less severe. Despite the aforementioned limitations due to the potentialflow solver used, we generally find good agreement between experiments and numerics. Figure 2 shows a representative time series of the free surface η for the 1:1 slope case at several gauges. Given the differences between experiment and numerics, the agreement is generally good. Waves are well aligned in time and the magnitudes are adequately captured. The differences between the numerically generated waves and experimental results at Gauge 1 can be attributed to various causes. We assume that all the wave component measured at Gauge 1 during experiments are propagating downstream, whereas they are actually composed of both incident waves and reflected waves from the step and from the beach during the experiments. The numerical results measured at the end of relaxation zone also include reflections, especially from the slope, which also contribute to the deviation. We have applied a double relaxation zone to absorb these reflections, but there is also a compromise between being able to make the exact waves and reducing the influence of reflected waves. There also seems to be a systematic horizontal shift especially at gauge 11. One of the reasons for this can be the minor error between experiments and simulations in terms of the exact location of gauges or the slope. Bed and sidewall frictions could also contribute to horizontal shifts at far away gauges.
Next, we look at the significant wave height H s at our four gauges. The significant wave height is defined by the mean of the largest 1/3 of the waves and the individual wave is determined by zero up crossing. The significant wave height of numerical and experimental results are shown in Fig. 3. In terms of the experimental and numerical results at gauge 1, 3, 8, and 11, they are shown in Table 2. Our simulations are based on the experimental results at gauge 1. The 1% difference in H s at gauge 1 is a result of combined reasons mentioned above, and the discrepancy at all locations is small with both experiment and numerics showing similar trends.
We now turn to height statistics. Figure 4 presents the probability that a given wave height is exceeded at a particular location. We include both fully non-linear simulations as well as otherwise identical simulations using the linear model. The agreement between non-linear simulation and experiment is generally good, showing similar departures from the Rayleigh distribution and linear simulation after the step. The Rayleigh distribution without modification is based on a narrow-band assumption and it is the most commonly adopted distribution for the crest to trough wave heights as well as the starting point for many other models. The exceedance probability of wave height follows the Rayleigh distribution for narrowbanded linear waves. The Rayleigh distribution is given by  where m 0 is the zeroth-order moment of the spectrum and this equation describes the probability of a wave height is bigger than a given wave height H . Although within the confidence intervals, the numerical simulation does seem to slightly under-predict the number of large waves at gauge 8 (x = 0.90 m) and more clearly (although still within confidence intervals), does so at gauge 11 (x = 5.00 m) some way after the change in depth. This contrasts with a smaller over-prediction by the non-linear model at gauges 1 and 3. Given the multiple minor differences in setup between experiment and simulation, it is hard to identify the reason for these small differences. However, agreement is sufficient to have confidence in our simulation. There is also an unexpected discrepancy between the linear simulations and the Rayleigh distribution at gauge 11.
We also consider the average shape of an extreme event. Although, here, we are primarily using this for validation, this is, to our knowledge, the first time this has been examined for the wave over a steep slope problem. In our analysis, we are focused on the average shape of an extreme event at gauges 1, 3, 8, and 11. Based on our simulation results, there are approximately 12,000 waves at each gauge in both the 1:1 and 1:3 slope cases. We find the largest 40 waves measured at each gauge and find their average shape. The number of waves chosen is a trade-off between capturing the extremes and statistical variability. Under a theoretical linear model, the expected shape of an extreme event is given by the scaled autocorrelation function-the 'NewWave' (Lindgren 1970;Boccotti 1983;Tromans et al. 1991) in which E( f ) is the spectrum, and m 0 is the zeroth-order moment of the spectrum. f is the peak frequency, and t is time. η(t) indicates the NewWave profile. The linear model predicts a symmetric wave packet around an extreme event.
Deviation from the linear model is expected at all gauges as bound harmonics will significantly modify the shape. However, we include the linear-model prediction for reference. We analyse the resulting shape in Fig. 5. At the input gauge 1, there is a minor asymmetry in the numerical simulations. At gauge 3 (at the top of the slope), both numerics and simulation suggest some minor asymmetry with the trough preceding a large wave being slightly deeper than that following it. Gauge 8, roughly half of the wavelength after the top of the slope, shows some clear discrepancies from the standard 'NewWave shape'. The trough preceding the largest wave is relatively sharp, whereas, at this location, the trough after is far flatter. This is consistent with the shape observed for deterministic wavegroups and has been explained by the release of bound waves at the step (Li et al. 2021b). Both experiments and numerics show similar results, suggesting strongly that the simulations are capturing the physics of bound wave release. The results at the final gauge (gauge 11), after the point we are primarily interested in, show the waves returning to a more standard symmetrical shape, although some deviation from this can still be seen in both experiments and numerics. There are also some significant differences between the shapes predicted computationally and observed experimentally at gauge 11 (as results for wave height statistics also suggested). It is unclear why there is this difference. However, gauge 11 is beyond the main area we are interested in, so, whilst we note the discrepancy, we do not think this should invalidate the main conclusions of this paper.

Results
We consider the engineering aspects of the enhanced rogue wave activity at the top of slopes. There are many different wave/structure interactions one could consider all with their different characteristics. In the present paper, we consider a simplest model of waves interacting with a column as this is a useful general example. For simplicity, we follow the general approach taken recently by Klahn et al. (2021) of using a model based only on the fluid kinematics so as to avoid having to solve any complex wave/structure interaction problem. The choice of this model is not critical to the conclusions of this paper-the aim is to do a calculation which relies on the Fig. 6 Comparison of moment statistics for slope 1:1 and slope 1:3 at gauge 3 and 8. Two gauges are shown with 90% confidence intervals. Moment is the total moment around the sea-bed after the slope kinematics rather than the free surface considered in most studies.
We model the structure using the drag term of the Morison's equation (Morison et al. 1950). In this case, a monopile of a diameter (D) of 0.5 m (D/λ = 0.225) is selected, and the load per unit length on the monopile is a function of the undisturbed horizontal velocity, which is given by where F(x, t) is the force per unit length on a fixed vertical cylinder located at a specific spatial location x, ρ is the water density, C D is the drag coefficient (C D = 1.2), and u is the horizontal velocity at corresponding time and location. We choose to analyse the moment around the sea-bed as this is probably the most important feature for the design of a typical offshore structure. From the OW3D simulations, we calculate the horizontal velocity. We can calculate the total moment around the sea-bed by integration before calculating the peak moment in each loading cycle. We only consider this after the top of a slope as the different water depths before this point would confuse the analysis as the change in the seabed depth would dominate the change in moments. In Fig. 6, moment exceedance is presented for gauges 3 (top of slope) and 8 (slightly after top of slope) and both 1:1 and 1:3 slopes showing both linear and non-linear simulation results. The linear simulations predict very similar moments in all cases and these are always lower than that predicted by the nonlinear model. However, the departure from the linear is very much larger for gauge 8 than just at the top of the slope at gauge 3. This difference is greater for the steeper slope (1:1 case). Interestingly, at gauge 3, there is a bigger difference between linear and non-linear for the milder slope (1:3 case), likely due to the different phase shifts of second-order bound and free waves between the two cases. At the top of a slope, a steeper slope may lead to the second-order waves more out of phase and, therefore, smaller amplitude of second-order waves. As a consequence, it would lead to smaller differences between the non-linear and linear simulations for a steeper slope case, as noted. This finding is consistent with the observations in both Li et al. (2021b) and Zheng et al. (2020) due to the following. The former found both a larger skewness and kurtosis at gauge 3 for the case of 1:3 slope (i.e., their Figs. 2 and 3) compared to the case of 1:1 slope. Zheng et al. (2020) (i.e., their Table VI) suggests a slope gradient plays an important role in the location of the peak kurtosis and skewness relative to the top of a slope, which is the furthest and nearest for the steepest slope (a step) and the mildest slope, respectively.
To analyse extremes, we consider the overturning moment which corresponds to the 1-in-1000 wave event. Thus, in Fig. 7, the moment is the moment corresponding to a probability of 10 −3 in the moment exceedance figure as a function of spatial location.
We see very similar behaviour for both slopes except in the non-linear simulations very near the top of the slope where higher loads are predicted for the flatter (1:3) case. This may be the result of a slightly higher significant wave height for Fig. 7 Overturning moment about the sea-bed predicted by the numerical model. 90% confidence intervals are shown. The moment is the moment corresponding to the probability of 0.1% (one-in-1000 wave) in the moment exceedance figure for different spatial locations. The locations of gauge 3 (At the top of the slope) and gauge 8 (0.9 m after the top of the slope) are indicated by dotted black lines the flatter case due to less breaking in OW3D. However, the moment behaviour is consistent with the previous work on this setup where very similar physics was found for the two slopes. There appears to be enhanced moment around the top of the step, which increases approximately 20%, consistent with the increased number of large waves observed in both simulation and experiment. The 'peak' in this lasts for several wavelengths beyond the top of the slope before decreasing and plateauing. The non-linear simulations consistently predict larger moments than the linear ones. However, the linear simulations do not produce a significant peak and are essentially uniform across the region considered. Thus, we can have confidence that the peak in the moments observed is the result of non-linear physics.

Conclusions
This paper has related the theoretical problem of excess rogue waves occurring at the top of slopes to a practical problem in ocean engineering. In doing so, this enables us to draw the following conclusions regarding the wave kinematics of the problem.
We have run a fully non-linear potential-flow model which directly solves for the wave kinematics. This model has been found to be in good agreement with experimental results for surface elevation including for the expected shape of an extreme wave events which to our knowledge has not previously been studied. We then use these kinematics to drive a simple wave loading model. We find enhanced loads for several wavelengths after the abrupt depth transition. This study suggests that extra care is needed when designing ocean infrastructure at the top of slopes and that the abnormal free surface wave statistics previously observed at the top of slopes also leads to unusual moment or unusual kinematics occurred at the top of the slope.