Statistics of Simulated Storm Waves over Bathymetry

: This paper shows simulations of high waves over different bathymetries to collect statistical information, particularly kurtosis and crest exceedance, that quantiﬁes the occurrence of exceptionally extreme waves. This knowledge is especially pertinent for the design and operation of marine structures, safe ship trafﬁcking, and mooring strategies for ships near the coast. Taking advantage of the ﬂexibility to perform numerical simulations with HAWASSI software, with the aim of investigating the physical and statistical properties for these cases, this paper investigates the change in wave statistics related to changes in depth, breaking and differences between long- and short-crested waves. Three different types of bathymetry are considered: run-up to the coast with slope 1/20, waves over a shoal, and deep open-water waves. Simulations show good agreement in the examined cases compared with the available experimental data and simulations. Then predictive simulations for cases with a higher signiﬁcant wave height illustrate the changes that may occur during storm events.


Introduction
This paper contributes to the study and statistics of storm waves by analysing the results of long-term simulations of high waves. To investigate the usefulness of specific model computations, a comparison with experimental data should include a quantitative comparison of statistical properties, particularly kurtosis, as a measure of occurrence for extreme wave events.
Images obtained from X-band radars on a ship or from the coast can be processed to detect the bathymetry and collect extensive wave data and current over an area of 5 km 2 , following the methodology developed more than 30 years ago, mainly using spectrum analysis [1,2]. At present, more deterministic methods of analyzing the radar images [3] make it possible to predict waves from sailing ships and to determine the bathymetry near coastal areas. This provides a simple way of collecting detailed wave information over large areas of several square kilometers, and increasing our knowledge of realistic coastal waves. Nevertheless, laboratory experiments and numerical simulations based on advanced wave theories are still important methods of quantifying the environmental wave properties in a well-controlled, systematic way, including their statistics. This paper aims to obtain improved insight into the physics, with simulations of extreme wave cases for which there are no experimental/previous simulations. We investigate the statistics of extreme waves over bathymetry in terms of the skewness, asymmetry, kurtosis, and crest exceedance.
The methodology adopted here is to justify the results of predictive simulations for which no laboratory or numerical simulation data are available, as well as is reasonably possible. Before executing a predictive simulation, simulations with different parameter settings (grid size, value of breaking parameter, bottom roughness, etc.) produce results that indicate sensitivity to parameter variations. Provided that changing the parameters of the simulations only causes small deviations in the agreement between simulation results and the available data, the most optimal setting provides a good guide for the ultimate simulation, with results that are accepted and presented.
In this paper, three different cases of wave transformation are studied: wave run-up to the coast, waves over a shoal, and waves in deep water.
The first case, in Section 3, deals with one-dimensional (1D) wave run-up on a slope of 1/20. Starting with extreme waves with a steepness of 0.06 and 0.08, we then take more extreme cases with a steepness of 0.13. Data are available from experiments in the BGO-First offshore tank, located in la Seyne sur Mer, as reported in [4]. The good agreement of the simulation for the most extreme experiment provides information for the set-up used to design successive predictive simulations for more extreme waves; the wave transformations and their statistical properties are then investigated.
The second case, in Section 4, concerns waves over a shoal. From a series of experiments with non-breaking waves [5], we successfully compare a 1D simulation with the result of the experiment with the highest kurtosis. Then, we investigate predictive simulations of higher waves that will break over the shoal. To investigate the effect of directionality, we also simulate short-crested waves over the shoal and observe the effect on the statistics.
In Section 5, the third case uses simulation results for short-crested waves in an extreme open sea, called the Draupner sea in [6], for which exceedance data are available. When restricted to 1D waves, the simulations show an essential change in the wave statistics, mainly caused by much more severe breaking.
For the simulations, HAWASSI software was used, based on the Hamiltonian formulation for inviscid and irrotational fluid flow above bathymetry, and has a pseudo-spectral implementation and kinematic breaking condition. Some details are provided in the next section.
In the discussion section, we will comment on the obtained results.

HAWASSI-Modelling and Simulation
HAWASSI software was developed in recent years for the efficient modelling of waves and ship motions in laboratories and environmental situations. The basis of HAWASSI is the Hamiltonian formulation of the wave dynamics for inviscid, irrotational fluid [7,8] For 1D waves, with obvious extensions for short-crested waves, the free surface elevation is denoted by η(x, t), with x being the position and t the time. The value of the fluid potential Φ(x, z, t) at the free surface is denoted by φ(x, t) = Φ(x, z = η(x, t), t). The Hamiltonian is the sum of kinetic and potential energy.
H F (φ, η) = K(φ, η) + P(η), with K(φ, η) = 1 2 F η −D |∇Φ| 2 dzdx, and P(η) = gη 2 dx in which D = D(x) is the depth. The kinetic energy K is approximated to obtain an explicit, necessarily approximate, expression in (η, φ). In this paper, we use approximations for the second-and third-order accurate equations that are given for the corresponding approximation of the Hamiltonian by Here, ∂ t denotes time differentiation, and δ φ and δ η are the variational derivatives of the functional concerning φ and η, respectively. The kinetic energy functional can be written as where C is the phase velocity operator that depends on depth D and elevation η. For second-and third-order codes, the integrand is Taylor-approximated in η to the thirdand fourth-order, respectively [9], for the simulations of long-and short-crested waves in this paper. For numerical discretization, HAWASSI uses pseudo-spectral methods. For varying bottoms D = D(x), this leads to Fourier integral operators in the expression of the Hamiltonian and its derivatives. To reduce computational efforts, a partitioning of the Fourier integral operators over selected bottom depths is constructed, see [9][10][11][12]. The number of Fourier modes is chosen by considering the main wavelength.
For wave generation, a specified spectrum can be influxed from one or more given points in 1D, and from specified curves in 2D. This influxing uses the methods described in [13], including an adaptation zone for nonlinear effects to develop.
At the endpoints of the simulation interval in 1D, and at the boundary of a square simulation area in 2D, damping zones are applied to make waves vanishingly small, to prevent artificial physical connections from the sides caused by the spectral periodicity. For that reason, in cases of run-up, the run-up area above the water level is followed by a run-down to a zero level at the boundary.
The breaking of waves is modelled, as in [14]: for initiation, a kinematic condition is used, for which the value of U/C, the quotient of fluid velocity and local phase speed at the crest has to be specified as smaller than one; after initiation, the energy is reduced using an eddy-viscosity model. Most of the results in this paper were obtained from simulations with U/C = 0.8.
More than a hundred simulation results were successfully compared to the result of experiments in laboratory wave tanks. HAWASSI was extended to include ship motions [15,16] and is also used for bathymetry and wave reconstruction from X-band radar images; see, e.g., [17].

D Long-Crested Wave Run-Up on Coast
In this section, we consider the run-up of waves on a slope of 1/20 with 1D simulations. As a reference, to quantify the accuracy of the HAWASSI simulation, we started with the two most extreme cases of a laboratory experiment performed in the BGO-First offshore tank located in la Seyne sur Mer, as described in [4]. The available experimental data are spatially scaled, by a factor of 100, from the laboratory scale to describe the environmental situation. The two most extreme cases in the experiment have a Pierson-Moskowitz spectrum, with significant wave height Hs = 6 m and period Tp = 12 s (H6T12), and Hs = 7.5 m and Tp = 16 s (H7.5T16). The input signal with a wave steepness of 0.08 and 0.06, respectively, for the two cases, is selected from the measurement data at the first buoy (x = 640 m).
Taking some reported sagging of the bottom at several places into account, we changed the flat slope to a catenary-shaped bottom profile. This led to an improvement in the correlation between simulated and measured wave traces of from 5 to 10%, depending somewhat on the position of measurement and wave case.
Going beyond the wave heights of the experiment, two extreme cases with a steepness of 0.13 are simulated: irregular waves with Pierson-Moskowitz spectrum with Hs = 12 m, Tp = 12 s (H12T12) and Hs = 15 m, Tp = 16 s (H15T16). With this steepness value, some extreme waves (wave height > 2.5 Hs) still occur.
All simulations were conducted over the domain length of 3 km, of which 400 m length was used for the right and left damping zones. The grid size was taken as 3 m for H6T12, 4 m for H7.5T16, 5 m for H12T12, and 6 m for H15T16. The choice of grid size was based on the main wavelength at the influx point. We took approximately 16 points per wavelength, and the number of grid points was taken to be a power of two to support faster Fourier transform computation. For wave run-up, a smaller grid size is required, to accommodate shortwaves entering the coast.
The HAWASSI setting for the kinematic breaking condition was taken to be U/C = 0.8 for the Tp = 12 s cases; for the Tp = 16 s cases, the breaking condition was U/C = 0.6 and, in addition, a bottom friction 0.02 was applied. These values were chosen to obtain the optimal comparison results with the experimental data. For the predictive simulations, we used the insight we obtained based on our experience with the experiments.

Moderate Wave Run-Up: Comparison of Simulations with Experiment
For this case, and in the rest of this paper, we show the so-called "Wave-print" of a simulation; see Figure 1 for the two cases. The plot provides a compact overview of the simulation domain, and various characteristic results for wave quantities over the full domain for the full time trace.
MTT( ) = min ( , ) The time signal of 13,000 s, as measured in the experiment at a depth D = 76 m, contains approximately 1000 waves. After filtering to remove long reflected waves, the signal is used as the influx signal for the simulation. For case H6T12, except for some isolated cases of breaking starting at D = 17 m, all waves broke at D = 10 m, which is Hs/D = 0.61. A nearly full break for case H7.5T16 at D = 17 m is for Hs/D = 0.41.
The breaking criterion used is U/C = 0.8 for case H6T12. For H7.5T16, this value gives a slightly too high MTC near x = 1700 at depth D = 20 m; the result shown above is for s simulation with U/C = 0.6 over the whole area, but the reduction in the MTC only occurred in a small area, and close to the measurement values. For H7.5T16, we used bottom friction 0.02, which slightly reduced the run-up and offered a better comparison with the experimental spectrum for infra-gravity waves.
(a) (b) Figure 1. The wave-print for the two cases for which experimental data are available: in (a) for case H6T12, and in (b) for case H7.5T16, with bottom friction. At the left of the influx point, before the damping zone, the difference in the value of MTC and that of MTT show quite substantial reflected waves.
The graphs of the significant wave height (Hs), the maximal temporal crest (MTC) and minimal temporal trough (MTT), average temporal crest (ATC), average temporal trough (ATT), and the mean water level (MWL) were plotted, all calculated over the full time trace, and shown over the full simulation area. The area where 95% of all breaking took place is indicated at the bottom of the plot. For waves over a non-flat bottom, a schematic view of the bottom is given. If available, measurement data are also shown with circles, triangles, and crosses. Given the free surface elevation η(x, t), crest elevation, η c (x, t), and trough elevation η t (x, t) in a time interval t ∈ [t 1 , t n ], the MTC, MTT, ATC, and ATT are calculated as follows: The time signal of 13,000 s, as measured in the experiment at a depth D = 76 m, contains approximately 1000 waves. After filtering to remove long reflected waves, the signal is used as the influx signal for the simulation. For case H6T12, except for some isolated cases of breaking starting at D = 17 m, all waves broke at D = 10 m, which is Hs/D = 0.61. A nearly full break for case H7.5T16 at D = 17 m is for Hs/D = 0.41.
The breaking criterion used is U/C = 0.8 for case H6T12. For H7.5T16, this value gives a slightly too high MTC near x = 1700 at depth D = 20 m; the result shown above is for s simulation with U/C = 0.6 over the whole area, but the reduction in the MTC only occurred in a small area, and close to the measurement values. For H7.5T16, we used bottom friction 0.02, which slightly reduced the run-up and offered a better comparison with the experimental spectrum for infra-gravity waves.
To quantify the quality of the simulation, we use the data of measurements at several points along the slope to calculate the correlation between the simulated and measured time trace; in addition, we show for four statistical quantities the quotient of the simulation and the experimental data.
The correlation between the simulated, η simul (x, t) and measured time trace, η meas (x, t) at a certain position is computed as in which η meas and η simul are the mean of the simulated and measured time trace. The Skewness, a measure of crest-trough shape, and Asymmetry, a measure of left-right differences in a wave trace, η are computed ( [14]) as follows: where H is the Hilbert transform and . is the mean operator. Kurtosis, which is a measure of the occurrence of extreme waves, is defined by: The results presented in Table 1 indicate the quality of second-order HAWASSI simulations. Note that the nearly vanishing values for Skewness and Asymmetry in the experiment may lead to relatively large values for the quotient. Statistical quantities along the slope are presented in Figure 2a,b for the two cases, and also the Infra-Gravity (IG) wave spectra of H7.5T16 in Figure 2c at the positions of measurement. Table 1. Comparison between experiment and simulations at 6 positions: the correlation 'Corr' of the time traces and the quotients of statistical quantities of simulation and measurement for: 'Hs' significant wave height, 'Sk' skewness, 'As' Asymmetry, and 'Kurt' Kurtosis for the two wave cases H6T12 and H7.5T16.

Run-Up of Storm Waves
In this section, we provide results for two cases that are more extreme than the laboratory cases, namely, H12T12 (Hs = 12 m, Tp = 12 s) and H15T16 (Hs = 15 m, Tp = 16 s). In the simulation for H15T16, we used bottom friction 0.02, as in the laboratory case H7.5T16.
We begin to show the wave-prints in Figure 3a,b, respectively. Both cases have approximately the same position at the start of the breaking, but note that this is much earlier than for the moderate run-up cases of the experiment. The breaking leads to the reduction in the MTC after the breaking position, and results in somewhat lower heights in the shallower area. We can also observe from the deviation in Hs (see Table 2) that the evolution of the significant wave height of more extreme wave run-up differs from the moderate cases; Hs continues to decrease towards the shoreline.
As there are no measurement data for the two extreme cases, we show the numerical values of the simulations at six depths, towards the shoreline. For comparison with the two cases of the previous section, we provide the results of all four cases in Table 2. In addition, the steepness values for different depths are listed. The statistical quantities of waves, calculated for 1000 wave periods and for longer time traces of 4000 wave periods, were simulated. Both H6T12 and H7.5T16 showed that the longer simulation does not significantly change the statistics, so that 1000 waves can represent the four statistical quantities of Hs, Skewness, Asymmetry, and excess Kurtosis (Kurt-3). Figure 2c shows some discrepancies in the IG wave spectra of H7.5T16 compared to the measurement in the shallower area, but the overall comparison is satisfying; see

Run-Up of Storm Waves
In this section, we provide results for two cases that are more extreme than the laboratory cases, namely, H12T12 (Hs = 12 m, Tp = 12 s) and H15T16 (Hs = 15 m, Tp = 16 s). In the simulation for H15T16, we used bottom friction 0.02, as in the laboratory case H7.5T16.
We begin to show the wave-prints in Figure 3a,b, respectively. Both cases have approximately the same position at the start of the breaking, but note that this is much earlier than for the moderate run-up cases of the experiment. The breaking leads to the reduction in the MTC after the breaking position, and results in somewhat lower heights in the shallower area. We can also observe from the deviation in Hs (see Table 2) that the evolution of the significant wave height of more extreme wave run-up differs from the moderate cases; Hs continues to decrease towards the shoreline.   The four simulations with increasing significant wave height show a consistent increase in Skewness towards the shallower area. It can be observed that, for the same wave period, the enhancement of the significant wave height can lower the excess Kurtosis, because more extreme waves break earlier, and the occurrences of an extreme wave will be fewer. A graphic illustration of the statistical data over the whole area is given in Figure  4.  As there are no measurement data for the two extreme cases, we show the numerical values of the simulations at six depths, towards the shoreline. For comparison with the two cases of the previous section, we provide the results of all four cases in Table 2. In addition, the steepness values for different depths are listed.
The four simulations with increasing significant wave height show a consistent increase in Skewness towards the shallower area. It can be observed that, for the same wave period, the enhancement of the significant wave height can lower the excess Kurtosis, because more extreme waves break earlier, and the occurrences of an extreme wave will be fewer. A graphic illustration of the statistical data over the whole area is given in Figure 4.

D Long and Short-Crested Waves over a Shoal
The list of experiments of long-crested irregular waves, as reported in [5], deals with different configurations for non-breaking waves over an underwater shoal; see Figure 5. One of the cases in that experiment, identified as "Run 3", had the largest maximal kurtosis at the beginning of the shallowest part of the shoal. The shoal is symmetric above the flat bottom of 53 m deep; the shallowest part has a length of 160 m and is 11 m under the water level; the length of the slopes on each side to the bottom is also 160 m, with a slope approximately of 1:4. For this run, the waves have a Jonswap spectrum, with = 3.3 and peak period Tp = 11 s. The maximal wave height for which no breaking occurs was to be Hs = 2.5 m. HAWASSI simulations with the kinematic breaking criterion of U/C = 0.8 also did not have wave-breaking, and the statistics for the simulation almost exactly correspond with this experiment, as summarised in Section 4.1. By simulating higher waves that will break over the shoal, we investigate the effect on the statistics in Section 4.2.
For all cases, we used the breaking criterion U/C = 0.8. All simulations used the second-order nonlinear model, except for Hs = 2.5 m, Tp = 11 s. In case Hs = 2.5 m, Tp = 11 s, the third-order simulation provided a better fit to the peak of kurtosis, while, in other cases, there was no significant difference between second-and third-order.

D Long and Short-Crested Waves over a Shoal
The list of experiments of long-crested irregular waves, as reported in [5], deals with different configurations for non-breaking waves over an underwater shoal; see Figure 5. One of the cases in that experiment, identified as "Run 3", had the largest maximal kurtosis at the beginning of the shallowest part of the shoal. The shoal is symmetric above the flat bottom of 53 m deep; the shallowest part has a length of 160 m and is 11 m under the water level; the length of the slopes on each side to the bottom is also 160 m, with a slope approximately of 1:4.

D Long and Short-Crested Waves over a Shoal
The list of experiments of long-crested irregular waves, as reported in [5], deals with different configurations for non-breaking waves over an underwater shoal; see Figure 5. One of the cases in that experiment, identified as "Run 3", had the largest maximal kurtosis at the beginning of the shallowest part of the shoal. The shoal is symmetric above the flat bottom of 53 m deep; the shallowest part has a length of 160 m and is 11 m under the water level; the length of the slopes on each side to the bottom is also 160 m, with a slope approximately of 1:4. For this run, the waves have a Jonswap spectrum, with = 3.3 and peak period Tp = 11 s. The maximal wave height for which no breaking occurs was to be Hs = 2.5 m. HAWASSI simulations with the kinematic breaking criterion of U/C = 0.8 also did not have wave-breaking, and the statistics for the simulation almost exactly correspond with this experiment, as summarised in Section 4.1. By simulating higher waves that will break over the shoal, we investigate the effect on the statistics in Section 4.2.
For all cases, we used the breaking criterion U/C = 0.8. All simulations used the second-order nonlinear model, except for Hs = 2.5 m, Tp = 11 s. In case Hs = 2.5 m, Tp = 11 s, the third-order simulation provided a better fit to the peak of kurtosis, while, in other cases, there was no significant difference between second-and third-order. For this run, the waves have a Jonswap spectrum, with γ = 3.3 and peak period Tp = 11 s. The maximal wave height for which no breaking occurs was to be Hs = 2.5 m. HAWASSI simulations with the kinematic breaking criterion of U/C = 0.8 also did not have wave-breaking, and the statistics for the simulation almost exactly correspond with this experiment, as summarised in Section 4.1. By simulating higher waves that will break over the shoal, we investigate the effect on the statistics in Section 4.2.
For all cases, we used the breaking criterion U/C = 0.8. All simulations used the second-order nonlinear model, except for Hs = 2.5 m, Tp = 11 s. In case Hs = 2.5 m, Tp = 11 s, the third-order simulation provided a better fit to the peak of kurtosis, while, in other cases, there was no significant difference between second-and third-order.

Simulations Compared to Experiment
For comparison with the experimental observation, 1D wave simulations were performed on a scale with a factor of 100 in space. The total simulation domain of 2390 m includes damping zones of around 250 m at the boundaries, and is discretized by a grid size of 2.34 m (1024 grid points); for the total simulation time of 54,000 s, outputs are saved every 0.3 s. To investigate the effect of dimensionality, we also performed simulations for short-crested waves, with a spreading of approximately 20 degrees. The wave-print in Figure 6 provides an overview of the simulation results for the 1D long-crested waves and the short-crested waves. For the short-crested wave simulation, to reduce the computation time, the simulation is given over a smaller domain that still covers the shoal area. Since the waves have some directional spreading, there are waves that go out of the domain in a direction perpendicular to the wave main propagation; hence, the Hs slightly reduces along the propagation direction (see Figure 6b).

Simulations Compared to Experiment
For comparison with the experimental observation, 1D wave simulations were performed on a scale with a factor of 100 in space. The total simulation domain of 2390 m includes damping zones of around 250 m at the boundaries, and is discretized by a grid size of 2.34 m (1024 grid points); for the total simulation time of 54,000 s, outputs are saved every 0.3 s. To investigate the effect of dimensionality, we also performed simulations for short-crested waves, with a spreading of approximately 20 degrees. The wave-print in Figure 6 provides an overview of the simulation results for the 1D long-crested waves and the short-crested waves. For the short-crested wave simulation, to reduce the computation time, the simulation is given over a smaller domain that still covers the shoal area. Since the waves have some directional spreading, there are waves that go out of the domain in a direction perpendicular to the wave main propagation; hence, the Hs slightly reduces along the propagation direction (see Figure 6b). To simulate the 1D waves, we used the measurement data at the left-most probe position (−350 m) as the influx signal, and moved this to the position of wave-maker (−790 m). The length of the time trace of 54.000 s contains approximately 4900 waves; our simulations over two-times-longer time intervals did not show any differences in the statistics. As in the experiment, no breaking occurred in the simulation that used third-order nonlinearity, and kinematic breaking condition U/C = 0.8. A snapshot of the time trace is given in Figure 7. The short-crested waves are simulated over an area of 1200 m × 1200 m, including damping zones of 100 m around the boundaries, discretized with a grid size of 4.71 m in the x-axis (main wave propagation direction) and 9.45 m in the y-axis. The analysis data To simulate the 1D waves, we used the measurement data at the left-most probe position (−350 m) as the influx signal, and moved this to the position of wave-maker (−790 m). The length of the time trace of 54.000 s contains approximately 4900 waves; our simulations over two-times-longer time intervals did not show any differences in the statistics. As in the experiment, no breaking occurred in the simulation that used third-order nonlinearity, and kinematic breaking condition U/C = 0.8. A snapshot of the time trace is given in Figure 7.

Simulations Compared to Experiment
For comparison with the experimental observation, 1D wave simulations were performed on a scale with a factor of 100 in space. The total simulation domain of 2390 m includes damping zones of around 250 m at the boundaries, and is discretized by a grid size of 2.34 m (1024 grid points); for the total simulation time of 54,000 s, outputs are saved every 0.3 s. To investigate the effect of dimensionality, we also performed simulations for short-crested waves, with a spreading of approximately 20 degrees. The wave-print in Figure 6 provides an overview of the simulation results for the 1D long-crested waves and the short-crested waves. For the short-crested wave simulation, to reduce the computation time, the simulation is given over a smaller domain that still covers the shoal area. Since the waves have some directional spreading, there are waves that go out of the domain in a direction perpendicular to the wave main propagation; hence, the Hs slightly reduces along the propagation direction (see Figure 6b). To simulate the 1D waves, we used the measurement data at the left-most probe position (−350 m) as the influx signal, and moved this to the position of wave-maker (−790 m). The length of the time trace of 54.000 s contains approximately 4900 waves; our simulations over two-times-longer time intervals did not show any differences in the statistics. As in the experiment, no breaking occurred in the simulation that used third-order nonlinearity, and kinematic breaking condition U/C = 0.8. A snapshot of the time trace is given in Figure 7.  The short-crested waves are simulated over an area of 1200 m × 1200 m, including damping zones of 100 m around the boundaries, discretized with a grid size of 4.71 m in the x-axis (main wave propagation direction) and 9.45 m in the y-axis. The analysis data for the short-crested waves are taken as the average over an interval of 200 m, perpendicular to the main propagation direction, within a time interval of 14,400 s. We provide a comparison of the main statistical quantities of the two simulations with the measurement data in Figure 8.
for the short-crested waves are taken as the average over an interval of 200 m, perpendicular to the main propagation direction, within a time interval of 14,400 s. We provide a comparison of the main statistical quantities of the two simulations with the measurement data in Figure 8. Notable differences between the simulations for the long-and short-crested waves include the faster decay of Hs, similar to the experiment at the start of the down-run slope, and the smaller value of the kurtosis for short-crested waves above the flat part of the slope. We did further investigate the similar decay in the Hs, but the faster decay in the short-crested simulation must be caused by the directionality. This directionality effect is also shown in the lower kurtosis value.
To improve the understanding of the dynamics of long-crested waves running over the shoal, it is noted that the wavelength decreases from 179 m in the deep part to 107 m above the flat part of the slope. The theoretical phase speed of the waves decreases monotonically, from approximately 16 m/s before the start of the shoal to approximately 9.6 m/s at the flat part. The change in the group speed is much lower: from 9.5 m/s in the deep area to 8.6 m/s at the flat part of the slope, which is close to the value of the phase speed above the shallow depth.
Both the simulation and the experiment show an oscillating Hs (see Figure 9) with approximately three oscillations over the 160 m length of the flat part, of length 53 m, which is half of the wavelength, similar to a "period doubling" phenomenon. However, observe that the pattern starts slightly before the top of the shoal. Notable differences between the simulations for the long-and short-crested waves include the faster decay of Hs, similar to the experiment at the start of the down-run slope, and the smaller value of the kurtosis for short-crested waves above the flat part of the slope. We did further investigate the similar decay in the Hs, but the faster decay in the short-crested simulation must be caused by the directionality. This directionality effect is also shown in the lower kurtosis value.
To improve the understanding of the dynamics of long-crested waves running over the shoal, it is noted that the wavelength decreases from 179 m in the deep part to 107 m above the flat part of the slope. The theoretical phase speed of the waves decreases monotonically, from approximately 16 m/s before the start of the shoal to approximately 9.6 m/s at the flat part. The change in the group speed is much lower: from 9.5 m/s in the deep area to 8.6 m/s at the flat part of the slope, which is close to the value of the phase speed above the shallow depth.
Both the simulation and the experiment show an oscillating Hs (see Figure 9) with approximately three oscillations over the 160 m length of the flat part, of length 53 m, which is half of the wavelength, similar to a "period doubling" phenomenon. However, observe that the pattern starts slightly before the top of the shoal.
It is observed that the short-crested simulation improves the comparison with the measured Hs in the down-run area behind the flat part of the slope. The down-run area also shows a noticeable change in the crest exceedance plot compared to values before and at the top of the shoal, shown in Figures 10 and 11. It is observed that the short-crested simulation improves the comparison with the measured Hs in the down-run area behind the flat part of the slope. The down-run area also shows a noticeable change in the crest exceedance plot compared to values before and at the top of the shoal, shown in Figures 10 and 11. Figure 10 shows the crest-and trough-exceedance plots over the interval of interest. Observe that the Maximal Temporal Crest (MTC) and Minimal Temporal Trough (MTT) curves in the wave-print plots are the upper and lower boundaries of these exceedance plots, but lack information on statistical occurrence.
In Figures 10 and 11, the higher value of MTC above the flat part of the slope for the 1D long-crested simulation, compared to the experiment and short-crested simulation, is noticeable, as is the smaller value of troughs in the run-down area for both simulations compared to the experiment.   It is observed that the short-crested simulation improves the comparison with the measured Hs in the down-run area behind the flat part of the slope. The down-run area also shows a noticeable change in the crest exceedance plot compared to values before and at the top of the shoal, shown in Figures 10 and 11. Figure 10 shows the crest-and trough-exceedance plots over the interval of interest. Observe that the Maximal Temporal Crest (MTC) and Minimal Temporal Trough (MTT) curves in the wave-print plots are the upper and lower boundaries of these exceedance plots, but lack information on statistical occurrence.
In Figures 10 and 11, the higher value of MTC above the flat part of the slope for the 1D long-crested simulation, compared to the experiment and short-crested simulation, is noticeable, as is the smaller value of troughs in the run-down area for both simulations compared to the experiment.   It is observed that the short-crested simulation improves the comparison with the measured Hs in the down-run area behind the flat part of the slope. The down-run area also shows a noticeable change in the crest exceedance plot compared to values before and at the top of the shoal, shown in Figures 10 and 11. Figure 10 shows the crest-and trough-exceedance plots over the interval of interest. Observe that the Maximal Temporal Crest (MTC) and Minimal Temporal Trough (MTT) curves in the wave-print plots are the upper and lower boundaries of these exceedance plots, but lack information on statistical occurrence.
In Figures 10 and 11, the higher value of MTC above the flat part of the slope for the 1D long-crested simulation, compared to the experiment and short-crested simulation, is noticeable, as is the smaller value of troughs in the run-down area for both simulations compared to the experiment.   Figure 10 shows the crest-and trough-exceedance plots over the interval of interest. Observe that the Maximal Temporal Crest (MTC) and Minimal Temporal Trough (MTT) curves in the wave-print plots are the upper and lower boundaries of these exceedance plots, but lack information on statistical occurrence.
In Figures 10 and 11, the higher value of MTC above the flat part of the slope for the 1D long-crested simulation, compared to the experiment and short-crested simulation, is noticeable, as is the smaller value of troughs in the run-down area for both simulations compared to the experiment.

More Extreme Breaking Waves
We now consider three more extreme cases of 1D long-crested waves for which no data from measurements are available. The three cases are for Hs = 4 m, Tp = 11 s (same period as experiment, but two times higher Hs), Hs = 6 m, Tp = 12 s, and Hs = 7.5 m, Tp = 15 s, denoted by H4T11, H6T12, and H7.5T15, respectively, using the same numerical parameters as mentioned in the previous subsection. In all cases, the waves break at the shoal run-up before the top, as can be seen in the wave-prints; for the two most extreme cases, all waves break at nearly the same depth (see Figure 12a-c).

More Extreme Breaking Waves
We now consider three more extreme cases of 1D long-crested waves for which no data from measurements are available. The three cases are for Hs = 4 m, Tp = 11 s (same period as experiment, but two times higher Hs), Hs = 6 m, Tp = 12 s, and Hs = 7.5 m, Tp = 15 s, denoted by H4T11, H6T12, and H7.5T15, respectively, using the same numerical parameters as mentioned in the previous subsection. In all cases, the waves break at the shoal run-up before the top, as can be seen in the wave-prints; for the two most extreme cases, all waves break at nearly the same depth (see Figure 12a Waves will break in the run-up of the shoal. However, the breaking of 1D longcrested waves is substantially different than the breaking for short-crested waves. To investigate this, we also performed a short-crested simulation for wave case H6T12. The short-crested waves are simulated over an area of 1500 m × 1800 m, including damping zones of 100 m around the boundaries, discretized with a grid size of 5.88 m in the x-axis and 7.06 m in the y-axis. Similar to the previous subsection, the data analysis is carried out using the data average over an interval of 200 m perpendicular to the main propagation direction, within a time interval of 14,400 s. The wave-print of the short-crested wave simulation is shown in Figure 12d.
All 1D long-crested waves break at the start of the flat part and continue during rundown until the deep part. However, of all the short-crested waves passing an interval of 200 m perpendicular to the propagation direction, only approximately 20% of the shortcrested waves break, all at the start of the flat part; see Figure 13. Note also that the shortcrested breaking continued for much longer downstream, even after the end of the slope. Due to the directionality, the breaking mechanisms in long-and short-crested waves are slightly different. Waves will break in the run-up of the shoal. However, the breaking of 1D long-crested waves is substantially different than the breaking for short-crested waves. To investigate this, we also performed a short-crested simulation for wave case H6T12. The short-crested waves are simulated over an area of 1500 m × 1800 m, including damping zones of 100 m around the boundaries, discretized with a grid size of 5.88 m in the x-axis and 7.06 m in the y-axis. Similar to the previous subsection, the data analysis is carried out using the data average over an interval of 200 m perpendicular to the main propagation direction, within a time interval of 14,400 s. The wave-print of the short-crested wave simulation is shown in Figure 12d.
All 1D long-crested waves break at the start of the flat part and continue during run-down until the deep part. However, of all the short-crested waves passing an interval of 200 m perpendicular to the propagation direction, only approximately 20% of the shortcrested waves break, all at the start of the flat part; see Figure 13. Note also that the short-crested breaking continued for much longer downstream, even after the end of the slope. Due to the directionality, the breaking mechanisms in long-and short-crested waves are slightly different. The statistics plots of 1D long-crested waves in Figure 14 show similar patterns compared to case H2.5T11, except for Hs. The Hs of short-crested waves remains nearly constant above the flat part of the shoal and in the run-down area, while the Hs of 1D decreases monotonically approximately 20% from its highest value at the start of the flat to the lowest value at the end of the flat part. This decrease is much larger than for case H2.5T11, because there is much more severe breaking. The kurtosis in that area is also large for short-crested waves, even larger than for 1D long-crested waves in most of the area. Figure 15 compares the crest and trough distribution for 1D long-and short-crested waves over the shoal for H6T12.  The statistics plots of 1D long-crested waves in Figure 14 show similar patterns compared to case H2.5T11, except for Hs. The Hs of short-crested waves remains nearly constant above the flat part of the shoal and in the run-down area, while the Hs of 1D decreases monotonically approximately 20% from its highest value at the start of the flat to the lowest value at the end of the flat part. This decrease is much larger than for case H2.5T11, because there is much more severe breaking. The kurtosis in that area is also large for short-crested waves, even larger than for 1D long-crested waves in most of the area. Figure 15 compares the crest and trough distribution for 1D long-and short-crested waves over the shoal for H6T12. The statistics plots of 1D long-crested waves in Figure 14 show similar patterns compared to case H2.5T11, except for Hs. The Hs of short-crested waves remains nearly constant above the flat part of the shoal and in the run-down area, while the Hs of 1D decreases monotonically approximately 20% from its highest value at the start of the flat to the lowest value at the end of the flat part. This decrease is much larger than for case H2.5T11, because there is much more severe breaking. The kurtosis in that area is also large for short-crested waves, even larger than for 1D long-crested waves in most of the area. Figure 15 compares the crest and trough distribution for 1D long-and short-crested waves over the shoal for H6T12.  To summarize the results, Figure 16 shows the statistical quantities of all four cases of 1D long-crested waves. The most extreme wave heights are found just before or at the start of the flat part, with heights of 8, 12, and 15 m; in addition, the troughs are more extreme there, as seen in the wave-prints: from 4 and 7 to 10 m.  To summarize the results, Figure 16 shows the statistical quantities of all four cases of 1D long-crested waves. The most extreme wave heights are found just before or at the start of the flat part, with heights of 8, 12, and 15 m; in addition, the troughs are more extreme there, as seen in the wave-prints: from 4 and 7 to 10 m.
In Figure 16a, large changes in Hs over the horizontal flat part are seen, and the higher the Hs, the larger the decrease over that area. This is caused by the increased breaking, which reduces the wave height in the shallowest part and even above the downward slope. The other three statistical properties shown in Figure 16b

Short-Crested Open Sea Statistics Compared to 1D Long-Crested Waves
The appearance of rogue waves has attracted much attention in the last two decades [19,20]. For waves in open seas, a substantial stimulus was the availability of reliable measurements from the Draupner platform in 1995 in the North Sea, showing an 18 m high wave passing under, and partly damaging, the platform [21]. This motivated intensified studies on the nonlinear effects of wave propagation.
The influence of quasi-resonant four-wave interactions in long-crested nonlinear wave statistics was studied, and shown to compare successfully with wave tank data, showing major differences from Gaussian behaviour; see, e.g., [22], and local processes [19,23]. Extending this to more realistic short-crested waves, it was found that modulational instability does not play an important role. Instead, the main generation mechanism for rogue waves is the constructive interference of elementary waves enhanced by secondorder bound nonlinearities; see [24] for an extensive investigation of three rogue waves in the North Sea.
In this section, we compare the simulations of high short-crested waves with simulations of 1D long-crested waves. As a specific example, we use the short crested simulations of the so-called 'Draupner Sea' [6]. This open sea simulation, above a depth of 70 m, used the 2D spectrum obtained from data of the European Centre for Medium-Range Weather Forecast, rotated over an angle of 130 in the western direction to obtain the mean energy flux directed from North to South, as reported in [6]. The sea has wave parameters were In Figure 16a, large changes in Hs over the horizontal flat part are seen, and the higher the Hs, the larger the decrease over that area. This is caused by the increased breaking, which reduces the wave height in the shallowest part and even above the downward slope. The other three statistical properties shown in Figure 16b-d provide a summary of the statistical properties of the waves.

Short-Crested Open Sea Statistics Compared to 1D Long-Crested Waves
The appearance of rogue waves has attracted much attention in the last two decades [19,20]. For waves in open seas, a substantial stimulus was the availability of reliable measurements from the Draupner platform in 1995 in the North Sea, showing an 18 m high wave passing under, and partly damaging, the platform [21]. This motivated intensified studies on the nonlinear effects of wave propagation.
The influence of quasi-resonant four-wave interactions in long-crested nonlinear wave statistics was studied, and shown to compare successfully with wave tank data, showing major differences from Gaussian behaviour; see, e.g., [22], and local processes [19,23]. Extending this to more realistic short-crested waves, it was found that modulational instability does not play an important role. Instead, the main generation mechanism for rogue waves is the constructive interference of elementary waves enhanced by second-order bound nonlinearities; see [24] for an extensive investigation of three rogue waves in the North Sea.
In this section, we compare the simulations of high short-crested waves with simulations of 1D long-crested waves. As a specific example, we use the short crested simulations of the so-called 'Draupner Sea' [6]. This open sea simulation, above a depth of 70 m, used the 2D spectrum obtained from data of the European Centre for Medium-Range Weather Forecast, rotated over an angle of 130 in the western direction to obtain the mean energy flux directed from North to South, as reported in [6]. The sea has wave parameters were Hs = 12 m, Tp = 14 s, and had a very wide spectrum of 120 degrees, as shown in Figure 17. Hs = 12 m, Tp = 14 s, and had a very wide spectrum of 120 degrees, as shown in Figure 17.
The uni-directional spectrum in the main energy propagation direction to the south is also shown. The 1D long-crested waves at peak frequency have a wavelength of 270 m, and the phase and group velocities are 19.8 m/s and 12.4 m/s, respectively. This section investigates the differences between 1D long-and short-crestedness on the wave statistics of the crest exceedance. The short-crested waves were simulated over an area of 10.4 × 6 km 2 ; a smaller observation area was used for the statistics over an area of 6 × 2.5 km 2 , with at least 8000 waves passing each point of this domain.
To investigate the statistics for 1D long-crested waves, third-order simulations were performed over an interval of 19 km, excluding 500 m damping zones at the boundaries, with a grid size of 5 m. More than 10,000 wave periods (150,000 s) produced a total of 2.2 × 10 9 datapoints. The computation time, compared to the physical time, was 5% for the non-breaking wave simulation and 15% for the simulations with breaking waves, all using a desktop with Intel Core i9 @2.8 GHz processor.
For the short-crested sea simulation, no breaking took place, although the breaking condition was set as U/C = 0.8; hence for all short-crested waves, U/C did not exceed 0.8.
For 1D long-crested waves, three simulations were performed, with different settings for breaking: one simulation without breaking criterion, and two with breaking criterion U/C =0.9 and 0.8, respectively; in the last two cases, breaking took place, so long-crested waves with U/C > 0.9 occurred. Wave-prints are given in Figure 18; the plot of the wave profile at the time the highest wave occurred in the three 1D long-crested simulations is included. Note that the highest wave crests, in all cases, reached an elevation of more than 1.5 significant wave height, while the surrounding crest elevations were approximately the average crest value over the whole simulation time, indicating a freak wave occurrence.
In the 1D long-crested simulation, without breaking, an extreme event of four successive freak wave heights took place, caused by three successive waves. All four crest heights are greater than 1.6 Hs, namely, (referring to Figure 19 Figure 19. It is notable that the steepness of the highest long-crested wave is approximately 0.113, which is not exceptional, but the steepness of This section investigates the differences between 1D long-and short-crestedness on the wave statistics of the crest exceedance.
The short-crested waves were simulated over an area of 10.4 × 6 km 2 ; a smaller observation area was used for the statistics over an area of 6 × 2.5 km 2 , with at least 8000 waves passing each point of this domain.
To investigate the statistics for 1D long-crested waves, third-order simulations were performed over an interval of 19 km, excluding 500 m damping zones at the boundaries, with a grid size of 5 m. More than 10,000 wave periods (150,000 s) produced a total of 2.2 × 10 9 datapoints. The computation time, compared to the physical time, was 5% for the non-breaking wave simulation and 15% for the simulations with breaking waves, all using a desktop with Intel Core i9 @2.8 GHz processor.
For the short-crested sea simulation, no breaking took place, although the breaking condition was set as U/C = 0.8; hence for all short-crested waves, U/C did not exceed 0.8.
For 1D long-crested waves, three simulations were performed, with different settings for breaking: one simulation without breaking criterion, and two with breaking criterion U/C =0.9 and 0.8, respectively; in the last two cases, breaking took place, so long-crested waves with U/C > 0.9 occurred. Wave-prints are given in Figure 18; the plot of the wave profile at the time the highest wave occurred in the three 1D long-crested simulations is included. Note that the highest wave crests, in all cases, reached an elevation of more than 1.5 significant wave height, while the surrounding crest elevations were approximately the average crest value over the whole simulation time, indicating a freak wave occurrence.  In the 1D long-crested simulation, without breaking, an extreme event of four successive freak wave heights took place, caused by three successive waves. All four crest heights are greater than 1.6 Hs, namely, (referring to Figure 19 Figure 19. It is notable that the steepness of the highest long-crested wave is approximately 0.113, which is not exceptional, but the steepness of the crest height is exceptional, being approximately four times larger (0.45). Figure 19. The four successive freak wave heights observed in the simulation without breaking, shown with parts of their profile at different appearance times.
In Figure 20, the quantitative statistical properties are presented over the length of the simulation interval. The plots show that Hs slowly decreases over the length of the simulation interval: approximately 2, 3, and 6% for a simulation with no breaking and U/C= 0.9 and U/C = 0.8, respectively. The amplitude of the influx signal was designed in such a way that, for each breaking condition, the average Hs is close to 12 m over the length of the observation domain. The values of Skewness and Asymmetry are nearly independent of the value of the breaking criterion, but the Kurtosis decreases with a more stringent breaking condition. All the statistical values averaged from all buoys' data in 1D simulations are pretty close to the previously calculated 2D data ensemble statistics, averaged over eight buoys [6]; see Table 3. The skewness of the 1D case without breaking is closer to the 2D case, but In Figure 20, the quantitative statistical properties are presented over the length of the simulation interval. The plots show that Hs slowly decreases over the length of the simulation interval: approximately 2%, 3%, and 6% for a simulation with no breaking and U/C= 0.9 and U/C = 0.8, respectively. The amplitude of the influx signal was designed in such a way that, for each breaking condition, the average Hs is close to 12 m over the length of the observation domain. The values of Skewness and Asymmetry are nearly independent of the value of the breaking criterion, but the Kurtosis decreases with a more stringent breaking condition. In Figure 20, the quantitative statistical properties are presented over the length of the simulation interval. The plots show that Hs slowly decreases over the length of the simulation interval: approximately 2, 3, and 6% for a simulation with no breaking and U/C= 0.9 and U/C = 0.8, respectively. The amplitude of the influx signal was designed in such a way that, for each breaking condition, the average Hs is close to 12 m over the length of the observation domain. The values of Skewness and Asymmetry are nearly independent of the value of the breaking criterion, but the Kurtosis decreases with a more stringent breaking condition. All the statistical values averaged from all buoys' data in 1D simulations are pretty close to the previously calculated 2D data ensemble statistics, averaged over eight buoys [6]; see Table 3. The skewness of the 1D case without breaking is closer to the 2D case, but All the statistical values averaged from all buoys' data in 1D simulations are pretty close to the previously calculated 2D data ensemble statistics, averaged over eight buoys [6]; see Table 3. The skewness of the 1D case without breaking is closer to the 2D case, but the 1D case with breaking criterion U/C = 0.8 has the same significant wave height and provides the closest asymmetry (absolute value) and kurtosis with the 2D case. The number of breaking events (55,295) in the simulation with U/C = 0.9 is more than tripled (135,084) in the simulation with U/C = 0.8. The number of breaking events decreases somewhat with increasing distance from the influx point (see Figure 21), as can be expected with the slight decay of Hs downstream (see Figure 20). the 1D case with breaking criterion U/C = 0.8 has the same significant wave height and provides the closest asymmetry (absolute value) and kurtosis with the 2D case. The number of breaking events (55295) in the simulation with U/C = 0.9 is more than tripled (135,084) in the simulation with U/C = 0.8. The number of breaking events decreases somewhat with increasing distance from the influx point (see Figure 21), as can be expected with the slight decay of Hs downstream (see Figure 20). For the short-crested waves, it holds that U/C < 0.8 for all waves. The fact that the simulation of 1D long-crested waves could be carried out without breaking indicates that their crest and elevation exceedance will be much higher, which is quantified in Figures 5  and 6. For specified breaking conditions of 0.9 and 0.8, actual breaking took place, and the exceedance probability decreases. The graphs in Figures 5 and 6 suggest that the 1D longcrested waves with U/C = 0.8 have the crest exceedance most similar to the short-crested sea. Hence, using U/C = 0.8 in simulations of 1D long-crested waves will produce statistical results that are close to those for short-crested waves, at least for these wave conditions. In Table 4, we provide the Elevation Exceedance probability for the short-crested waves and for the 1D long-crested waves, as simulated with different U/C values in Figure 22. For the short-crested waves, it holds that U/C < 0.8 for all waves. The fact that the simulation of 1D long-crested waves could be carried out without breaking indicates that their crest and elevation exceedance will be much higher, which is quantified in Figures 5 and 6. For specified breaking conditions of 0.9 and 0.8, actual breaking took place, and the exceedance probability decreases. The graphs in Figures 5 and 6 suggest that the 1D long-crested waves with U/C = 0.8 have the crest exceedance most similar to the short-crested sea. Hence, using U/C = 0.8 in simulations of 1D long-crested waves will produce statistical results that are close to those for short-crested waves, at least for these wave conditions. In Table 4, we provide the Elevation Exceedance probability for the short-crested waves and for the 1D long-crested waves, as simulated with different U/C values in Figure 22. Table 4. Quantitative data of Elevation Exceedance for the four waves cases. Note that, comparing the short-crested waves with 1D long-crested waves, for η H s < 1.5 , the probability of the 1D longcrested waves with U/C = 0.8 is, on average, a factor two higher, but for η H s = 1.6 the probability of short-crested waves becomes much higher. (a) (b) Figure 22. Elevation exceedance plots for the short-crested Draupner sea (black), and the 1D long-crested waves without breaking (blue, highest), with breaking for U/C = 0.9 (green, middle) and with breaking U/C = 0.8 (red). (a) for the cumulative probability distribution, and (b) for the probability distribution. The lines in the elevation exceedance probability distribution plot are the crest exceedance probability of the Rayleigh distribution (solid green) and Forristall 3D (dashred). Table 4. Quantitative data of Elevation Exceedance for the four waves cases. Note that, comparing the short-crested waves with 1D long-crested waves, for < 1.5 , the probability of the 1D longcrested waves with U/C = 0.8 is, on average, a factor two higher, but for = 1.6 the probability of short-crested waves becomes much higher.

Discussion
In this paper, we have shown the possibility of studying high waves using numerical simulations with HAWASSI. For the interpretation of simulation results, it is most useful to present the main results in a compact and illuminating way; therefore, we designed the "Wave-prints" as an extension of plots in [25] to show the development of the most relevant quantities over the whole simulation area. In addition, we showed that the change in wave parameters, breaking phenomenon, and directionality can lead to an essential change in the statistics, especially Skewness, Kurtosis, and crest exceedance.
Comparing the simulations with the experiment for wave run-up to the shore and the wave over shoal showed that the simulation provided reliable results. This offers some Figure 22. Elevation exceedance plots for the short-crested Draupner sea (black), and the 1D long-crested waves without breaking (blue, highest), with breaking for U/C = 0.9 (green, middle) and with breaking U/C = 0.8 (red). (a) for the cumulative probability distribution, and (b) for the probability distribution. The lines in the elevation exceedance probability distribution plot are the crest exceedance probability of the Rayleigh distribution (solid green) and Forristall 3D (dash-red).

Discussion
In this paper, we have shown the possibility of studying high waves using numerical simulations with HAWASSI. For the interpretation of simulation results, it is most useful to present the main results in a compact and illuminating way; therefore, we designed the "Wave-prints" as an extension of plots in [25] to show the development of the most relevant quantities over the whole simulation area. In addition, we showed that the change in wave parameters, breaking phenomenon, and directionality can lead to an essential change in the statistics, especially Skewness, Kurtosis, and crest exceedance.
Comparing the simulations with the experiment for wave run-up to the shore and the wave over shoal showed that the simulation provided reliable results. This offers some confidence in the results for the two predictive simulations of more extreme storms. In the simulations, the initiation of breaking depends on the value of the kinematic quantity U/C. In all four cases of wave run-up to the shore, the deepest depth at which most waves broke was around Hs/D = 0.34; this is less than half the value of Hs/D = 0.8, which is usually defined as the start of the breaking zone. The earlier breaking waves occurred in the case of more extreme waves, in which an increase in significant wave height gave lower Kurtosis. In the case of H6T12, the breaking waves occurred later, and the Kurtosis increased towards the shoreline.
Even for the waves without breaking, high Kurtosis was measured in the area of wave transformation over the shoal. The experiment and simulation showed substantial changes in Kurtosis and Hs over the structure. The nonbreaking 1D long-crested waves of the experiment had a crest height of 2Hs (5 m) at the beginning of the flat part. Although the simulated value of Hs followed the value of the experiment during the run-up, the remarkable oscillations above the flat part were smaller in the simulation and might be the cause of the somewhat higher value of Hs in the successive run-down area, as shown in Figures 8a and 9. The results of more extreme and breaking waves over the shoal showed differences compared to short-crested waves. The main, large differences are a consequence of the occurrence of breaking. Smoothing the corners of the shoal had little effect.
For the Draupner sea, we showed that both the short-and the 1D long-crested simulations could run without breaking. For the short-crested sea, all waves had U/C < 0.8, while the 1D long-crested simulation contained substantial waves, with U/C > 0.8. Simulating without breaking leads to very high waves: the statistics of the 1D long-crested waves without breaking show an elevation exceedance far beyond known distributions, such as those described by Rayleigh or Forristall 3D distributions. However, using the mild kinematic breaking condition U/C = 0.9 or 0.8, the crest distribution of the 1D long-crested waves agreed quite well with the short-crested results, even though the short-crested waves have U/C < 0.8.
Concerning the methodology needed for predictive simulations to be able to provide sufficient confidence in the validity of the methods and methodology, in Appendix A, an example was given to illustrate the practical strategy that was used to deal with predictive simulations. For all the extreme cases shown in this paper, the robustness of the results was tested in a similar way, to achieve the best possible reliable results.
Finally, for practical considerations, it should be kept in mind that for most resultsnumerical and laboratory measurements-the highest waves will be higher than in nature because the effect of local wind gusts is ignored completely. In nature, the effect of local wind gusts may lead to additional breaking that reduces wave heights, even if the inclusion of gustiness may not be needed for storm surge calculations [26].

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author on reasonable request.

Acknowledgments:
The authors thank Molin for making available the measurement data of the run-up experiments in Section 3 and Trulsen for the measurement data of the shoal in Section 4. We are very grateful to Jang Whan Kim of Technip Energies to involve us in the run-up simulations of the Oceanide experiment.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
As an example of the challenges present in numerical simulations, we consider the case of run-up, H7.5T16. In Section 3.1, results are provided for the setting with breaking criterion U/C = 0.6, and with bottom friction 0.02. This was expected to be more 'realistic' ("long waves are more affected by the bottom in shallower water") than a simulation with standard settings U/C = 0.8, and no bottom friction. To make this reasoning more explicit, two different setting results are shown in Figure A1.
The comparison with experiments seemed to be reasonably 'good'. More precisely, referring to Figure A1: The simulation in a,b has too-high MTC, Hs, and Kurtosis between 20 and 30 m depth. The simulation shown in (c) with U/C = 0.5 seems to lead to better MTC values, but the other statistical quantities, particularly the Kurtosis at the two shallowest measurement points (0.4 simulated, 0.6. experiment), are too low. This led to the choice of U/C = 0.6 with some bottom friction.
We consider HAWASSI software as a 'modelling tool'; its relatively fast simulations support this methodology. The use of HAWASSI is, therefore, more practical for longterm simulations, as used here to study the statistics of extreme waves. Such simulations require more effort when carried out in fully three-dimensional CFD simulations with high computational costs (see [27]). There are many alternative Boussinesq models and codes that eliminate the vertical coordinate from the basic three-dimensional Euler equations; well-known examples that use a spectral implementation are HOSM codes, which are more efficient than CFD for waves in intermediate and deep water, but with limited abilities to perform wave run-up [28].
MTC values, but the other statistical quantities, particularly the Kurtosis at the two shallowest measurement points (0.4 simulated, 0.6. experiment), are too low. This led to the choice of U/C = 0.6 with some bottom friction.
We consider HAWASSI software as a 'modelling tool'; its relatively fast simulations support this methodology. The use of HAWASSI is, therefore, more practical for longterm simulations, as used here to study the statistics of extreme waves. Such simulations require more effort when carried out in fully three-dimensional CFD simulations with high computational costs (see [27]). There are many alternative Boussinesq models and codes that eliminate the vertical coordinate from the basic three-dimensional Euler equations; well-known examples that use a spectral implementation are HOSM codes, which are more efficient than CFD for waves in intermediate and deep water, but with limited abilities to perform wave run-up [28].
(a) (b) (c) (d) Figure A1. Wave print and statistical quantities over the area for wave case H7.5T16 for two simulation settings: (a,b) for simulation without bottom friction and U/C = 0.8, and (c,d) for U/C = 0.5, without bottom friction.