Numerical modelling of wave overtopping discharges at rubble mound breakwaters using OpenFOAM®

Numerical modelling of wave interaction with rock-armoured rubble mound breakwaters has been performed to study wave overtopping. The influences of the slope angle, a berm in the seaward slope, a protruding crest wall, a recurved parapet, and the wave steepness have been studied using a validated CFD model (OpenFOAM). The numerical modelling confirms trends that have been observed in physical model tests while the validity of earlier developed guidelines has been examined outside the ranges of the physical model tests on which the guidelines are based. The numerical model results confirm that wave overtopping at rubble mound breakwaters depends on the wave steepness, that the influence of a berm is affected by the wave steepness, and that an earlier developed influence factor to account for the effects of a protruding crest wall can be applied to even larger crest walls than the tested crest walls on which the guidelines are based. The results indicate that the influence of the applied core material of the berm on the discharges is very limited. The numerical model also indicates that applying a recurved parapet on a crest wall of a rubble mound breakwater only has an effect for very small overtopping discharges. The numerical model results show that wave overtopping at rubble mound breakwaters strongly depends on the slope angle. Since this effect is so large that it cannot be neglected, while present guidelines for non-breaking waves do not include the effect of the slope angle, modified guidelines have been proposed. The observed effects of the slope on wave overtopping discharges at rubble mound structures still need to be verified based on physical model tests.


Introduction
The hydraulic performance of coastal structures such as rubble mound breakwaters determines whether the structures meet the functional requirements.Wave overtopping at rubble mound structures is one of the most important phenomena affecting the hydraulic performance.In addition to the design of coastal structures, also the adaptation of coastal structures has become more important due to sea level rise.For coastal structures that are in relatively shallow water, sea level rise can increase the wave loading on the structure because less wave dissipation occurs before the waves reach the structure.Adding a crest wall to an existing structure, increasing the height of a crest wall, adding a berm, or increasing the width or height of a berm, can be effective measures to account for effects of sea level rise.For this purpose, the individual effects of a crest walls and a berm need to be predicted, but also the combination of both (see for instance Hogeveen, 2021).
Estimates of wave overtopping are generally based on physical modelling in wave flumes and wave basins.Numerical modelling of wave overtopping provides additional opportunities to examine wave overtopping for a wide variety of structure geometries.To provide design guidelines for rubble mound structures with a crest wall and for structures with a berm in the seaward slope, Van Gent et al. (2022) provides design guidelines based on physical model tests.Although their tests have been performed for a wide range of structure geometries, numerical modelling provides opportunities to examine wave overtopping at structures with a crest wall and a berm to further extend guidelines for the design and (climate) adaptation of rubble mound structures.
Wave interaction with permeable coastal structures has been modelled numerically using the non-linear shallow water equations (NSWE) by Kobayashi and Wurjanto (1989), Wurjanto and Kobayashi (1993) as well as by Van Gent (1992, 1994and 1995a).The influence of the shape of wave energy spectra and the influence of berms on wave run-up and wave overtopping have been studied using non-linear shallow water equations (see Van Gent, 1999, 2000).However, other models are required for processes and structure geometries for which non-hydrostatic pressures are important, such as for wave overtopping at vertical walls at the crest of rubble mound breakwaters.
The Reynolds-Averaged Navier-Stokes equations (RANS) that are solved using the Volume-Of-Fluid (VOF) method were first applied for permeable coastal structures in Van Gent et al. (1994) and Van Gent (1995a).For an overview of RANS-VOF models for applications with permeable coastal structures reference is made to Losada et al. (2016).One of the more recent RANS-VOF models is the OpenFOAM model applied by for instance Chen et al. (2021Chen et al. ( , 2022) ) to examine the influences of roughness, berms and oblique wave attack on wave overtopping at coastal structures.Jacobsen et al. (2018) and Irías Mata (2021) examined forces on crest walls at rubble mound breakwaters using OpenFOAM.Castellino et al. (2021) examined forces and wave overtopping at vertical walls with a recurved parapet using OpenFOAM.
In the present study, OpenFOAM is validated and applied to examine the influences of several aspects such as the wave steepness, structure slope, berm, crest wall and recurved parapet, on wave overtopping at rubble mound breakwaters.
The structure of the paper is as follows.In Section 2 the most essential aspects of the physical model tests used for validation of the numerical model are described as well as the set-up of the numerical model.In Section 3 the validation of the numerical model is described.In Section 4 the numerical model is used to examine the influence of various parameters on wave overtopping discharges.In Section 5 the numerical model results are discussed in relation to existing empirical expressions.Section 6 provides conclusions and recommendations.

Description of physical model tests for validation
Physical model tests in a wave flume by Van Gent et al. (2022) were used for validation of the numerical model.The rubble mound breakwater configurations used for validation were a) without a crest wall and without a berm, b) without a berm and with a crest wall, and c) with a berm and with a crest wall.Fig. 1 illustrates some of the configurations that were used for validation.All cross-sections had a 1:2 slope, a permeable core (D n50 = 16 mm) and an armour layer (D n50 = 38 mm) with a thickness of 2D n50 .The width and level of the berm were varied.The berm consisted of the armour material.The crest walls were positioned on top of the core material.No recurved parapet (bullnose) was applied on the crest wall.
The incident wave height and wave steepness (s m-1,0 = 2π H m0 /gT m-1,0 2 ) were varied.In all tests a JONSWAP wave spectrum was used.All tests consisted of 1 000 waves.The mean overtopping discharges were measured.
Based on the physical model tests, the following expression was proposed to describe wave overtopping at rubble mound breakwaters: where the mean wave overtopping discharge q (m 3 /s/m) is made nondimensional by using Q = q/(gH m0 3 ) 0.5 , g is the acceleration due to gravity (m/s 2 ), R c is the freeboard (including the height of a crest wall, if present) relative to the still water level (m), H m0 is the spectral significant wave height of the incident waves at the toe of the structure (m), γ p is the influence factor for a recurved parapet on a crest wall (− ), and for the influence of roughness (γ f ), the influence of a crest wall (γ v ), the influence of a berm (γ b ), and the influence of oblique waves (γ β ), the following expressions were proposed: where D n50 is the diameter of the stones in the armour layer, R c -A c is the protruding part of the crest wall, B is the width of the berm, B L is the vertical distance between the level of the berm and the level of the armour layer at the crest, and β is the angle of wave incidence (β = 0 • for perpendicular wave attack).
The ratio of the protruding part of the crest wall and the freeboard (R c -A c )/R c was smaller than 0.35 in the tests on which the expression for the influence factor for the crest wall was based.This leads to a maximum influence factor for the crest wall of γ v = 1.16.The influence factor for the berm was within the range 0.5 ≤ γ b ≤ 1 based on a berm widths shorter than B/H m0 ≤ 5.24 and berm levels in the range 0 ≤ B L /H m0 ≤ 2.09.All data was obtained with slope angles of 1:2.Outside the tested ranges the validity of the expressions (Eqs.(1)-( 5)) remains unknown.Note that in the expression for wave overtopping discharges (Eq.( 1)) no influence of the slope angle is incorporated since all tests were performed with 1:2 slopes.Earlier expressions to described wave overtopping at rubble mound breakwaters also do not account for a potential influence of the slope angle (e.g.TAW, 2002, or EurOtop, 2018).

Description of numerical model set-up
The numerical model applied here consists of the coupling of the CFD model named OpenFOAM (Weller et al., 1998)  and outlet boundaries to generate and absorb waves.The offshore wave conditions until the vicinities of the breakwater were simulated in OCW3D.The wave breaking and the wave-structure interaction was computed by OpenFOAM.

Hydrodynamic model
The hydrodynamic model is based on the volume averaged twophase version of the Navier-Stokes equations.To account for permeable structures reference is made to Van Gent (1995a,b) and Jensen et al. (2014).In the Navier-Stokes equations, the velocity in the permeable parts is defined as the filter velocity (u f ).The filter velocity is the spatial averaged velocity over a certain area including the area occupied by the stones.Thus, the momentum and continuity equation are written as: where C m is the added mass coefficient that accounts for the transient interaction between grains and water, u f is the filter velocity vector in Cartesian coordinates, n is the porosity of the permeable structure, ∇ = y, z] is the Cartesian coordinate vector, μ is the dynamic molecular viscosity and F p accounts for the resistance force due to the presence of the porous media.The tracking of the free surface between the air and water interface is performed by the Volume-of-Fluid approach (VOF), see Hirt and Nichols (1981).This method captures the free surface by adding an advection equation for the indicator function of the VOF field.The indicator function F is used to estimate the fluids densities and viscosities: where the subscripts w and a refer to water and air, respectively.The fluids densities and viscosities are weighted average based on the distribution of water and air in each cell.To keep a sharp interface between water and air and to avoid wiggles around the free surface and overestimated velocities near the wave crests, the isoAdvector algorithm developed by Roenby et al. (2016) and suggested by Larsen et al. (2019) is used.The OpenFOAM version applied here is the OpenFOAM model including the waves2foam toolbox developed by Jacobsen et al. (2012) and used by Jacobsen et al. (2018), including the algorithm proposed by Roenby et al. (2016).

Mesh
The grid in OpenFOAM was generated by applying the utilities blockMesh and snappyHexMesh.Square cells were used as recommended by Jacobsen et al. (2012) since they give better solutions when modelling non-linear waves.A base mesh of 0.04 m in x and z direction was created using the blockMesh utility.Then, the mesh was only refined in the necessary regions where a good estimation of the physical processes is vital with the snappyHexMesh utility.They were prescribed near the water surface and around the crest wall.The refinement around the water surface was included to have nearly 8 to 11 cells per wave height, similar to the mesh defined by Jacobsen et al. (2018) and Chen et al. (2021), who defined 13 to 19 cells per wave height and 9 to 14 cells per wave height, respectively.In this way, the wave propagation and wave breaking processes are modelled as accurate as possible since they have a direct impact on the wave overtopping.Around the crest wall, an extra level of refinement was applied in order to model the flow around the crest wall in more detail, leading to cell sizes of 0.01 m in x and z direction, following the same approach taken by Jacobsen et al. (2018) and Molines et al. (2019), who applied two levels of refinement around the crest wall.Also, simulations with a recurved parapet were performed, and in those computations an even higher level of refinement was included to reproduce properly the curved part of the parapet.The impermeable crest wall was removed from the numerical domain by applying the snappyHexMesh utility.A summary and a visualization of the final mesh is shown in Fig. 3 and Table 1, which results from a mesh sensitivity study (see more details in Section 3.3).

Boundary conditions
Relaxation zones are applied at the inlet and outlet boundary to generate and absorb waves.Slip velocity conditions are included for the flume bottom and for the other impermeable parts.The upper boundary, i.e. the atmosphere, is treated as an open boundary where water and air can leave the domain but only air may enter the domain again.

Wave generation
In the numerical flume, the waves were generated in OCW3D by inputting the wave paddle velocity signal measured during the physical tests.Then, through the first relaxation zone (also known as coupling zone) displayed in Fig. 2, the wave signal from OCW3D was provided to the OpenFOAM model.

Flow resistance and turbulence
Turbulence is an important process in coastal waters.Wave breaking at breakwaters induces different levels of turbulence outside and inside the porous media.Van Gent (1995a,b)   M. Irías Mata and M.R.A. van Gent version of the Navier-Stokes equations where a constant eddy-viscosity is applied outside the permeable parts and the dissipation of energy in the permeable parts is accounted for by Forchheimer-type of resistance (F p in Eq. ( 6)).
The Forchheimer equation is applied to solve for the resistance force caused by the porous structure: The Forchheimer equation includes a linear and a non-linear term.When the first term dominates, the flow behaves as laminar (i.e.Darcy flow) while if the second term dominates, the porous media flow is turbulent.The parameters a and b are resistance coefficients for which the parameterisation by Van Gent (1995a,b) is used (with coefficients α F = 1 000 and β F = 1.1): Here, α F and β F are closure coefficients that depend on the grading and shape of the grains, ν is the kinematic viscosity and KC is the Keulegen-Carpenter number that accounts for effects of the oscillatory porous flow.It has been observed by several authors that it is possible to predict the bulk hydrodynamics accurately (such as wave height, wave transformation and dissipation through porous layers) without the inclusion of a detailed turbulence model (Van Gent, 1995a, 1995b;Jensen et al., 2014, Jacobsen et al., 2018;Molines et al., 2019, Irías Mata, 2021).The same approach is followed in this study.This means that there is no detailed turbulence model inside or outside the permeable structure because the same system of equations holds for the entire domain (inside and outside the porous media).Inside the permeable structure the Forchheimer equation applies, which already includes the effect of dissipation of wave energy due to turbulence because the resistance coefficients α F and β F are based on experimental results.
Nonetheless, a sensitivity analysis of the results due to the modelling of porous media flow is conducted and explained in Section 3.3.

Postprocessing of wave overtopping discharges
Wave overtopping discharges were obtained by placing an over-topping sheet ξ over a group of cell faces f on top of the crest wall, as illustrated in Fig. 4. Note that the base-plate in the model is narrow compared to actual base-plates of crest walls, but since the discharge at the top of the crest wall is the key-parameter here, the narrow base plate and the relatively large grid cells in the region at the back of the crest wall are considered not to affect the computed discharges.The model captures the flux of fluids going through each of the cell faces f.Each face fluid flux is then multiplied by the fraction of water to obtain the overtopping discharge per cell face Φ F,f .The total discharge is the sum of the water fluxes Φ F,f through all the cell faces f that are part of the overtopping sheet ξ.The instantaneous water flux is estimated by: where S f is the non-unit normal vector to the face (Jacobsen, 2017).

Validation of numerical model
In this section, the numerical model is validated by comparing the incident wave conditions and the overtopping discharges for a breakwater with a berm and a crest wall.17 cases were simulated in Open-FOAM, where four different wave conditions were included: 1) H m0 = 0.22 m and s m-1,0 = 0.013, 2) H m0 = 0.22 m and s m-1,0 = 0.026, 3) H m0 = 0.17 m and s m-1,0 = 0.013, and 4) H m0 = 0.17 m and s m-1,0 = 0.026.Hereafter, the first two conditions are referred to as "high waves" and the last two conditions as "low waves".For each wave condition variations of the berm width were included from structures   without a berm to structures with berm widths of B = 0.25 m, 0.5 m and 0.75 m.Each simulation consisted of 1 000 waves, as during the physical model tests.In all these validation cases the freeboard was R c = 0.3 m and the height of the protruding part of the crest wall was R c -A c = 0.05 m.One condition was added without a protruding crest wall (R c = A c = 0.25 m).First, the surface elevation time series and wave energy spectrum are compared against the experimental results, showing the capabilities of OpenFOAM to reproduce the wave conditions generated during the physical modelling campaigns.Then, the time series of wave overtopping discharges and the non-dimensional overtopping discharges are compared against the experimental results to assess the performance of OpenFOAM to predict both the events and magnitude of the overtopping discharges at breakwaters with crest walls and varying berm widths.

Wave conditions
The incident water surface elevation was computed in the physical and numerical flume at wave gauge WG1 (see Fig. 2) located 31.32 m from the wave maker.The measurements from the five gauges displayed in both physical and numerical flumes were used to decompose the incoming and reflected waves by applying the method of Zelt and Skjelbreia (1992).
Fig. 5 shows a comparison between the measured and computed time series of the incident water surface elevations for simulations for high waves with a relatively high wave steepness and for low waves with a relatively low wave steepness.The Pearson correlation coefficients are of 0.85 and 0.91, respectively.For the 17 simulations, the correlations are always higher than 0.83.It is concluded that the surface elevations are well captured by the numerical model.
Fig. 6 shows comparisons between measured and computed wave energy spectra.Both comparisons indicate that the model provides accurate reproductions of the wave energy spectra.A close look to the spectrum parameters for all 17 simulations reveals an average error of 0.7% and 8.6% for the H m0 and the T m-1,0 , respectively.However, it should be noted that the computed wave periods are systematically higher than those measured, leading to a wave steepness that is for these 17 simulations on average 15.6% lower than the measured wave steepness.Since a lower wave steepness is expected to lead to higher overtopping discharges, computed discharges are affected due to this underestimation of the wave steepness.

Wave overtopping
To verify the capabilities of OpenFOAM in predicting overtopping discharges at breakwaters, the measured and computed cumulative overtopping discharges were compared for all 17 simulations.In this section, two simulations with different overtopping volumes are illustrated in Fig. 7.In both cases, the OpenFOAM model overestimates the wave overtopping discharge (by a factor two and five, respectively).A more detailed inspection of the time series reveals that for some overtopping events larger overtopping volumes are computed and that for some conditions more overtopping events are captured in the numerical flume than in the physical flume.Nonetheless, the overtopping events present in the physical tests are mostly captured by the numerical model.
The comparison between the non-dimensional overtopping discharge for the experimental and numerical simulations is shown in Fig. 8. Overall, there is an overestimation of the wave overtopping discharges in OpenFOAM for these validation cases.The bias (average ratio of predicted over measured non-dimensional discharges) for these cases is a factor 4.3, i.e. the numerical model on average overestimates the measured discharges by a factor 4.3 in the 16 cases for which overtopping occurred; denoted by the dashed line in Fig. 8 with a RMSE value of 0.237 where the RMSE value is defined as: where n tests is the number of tests on which the RMSE is based.Note that differences between measurements and prediction methods also occur for other types of prediction methods like empirical expressions and machine learning techniques.For the mentioned 16 cases the method proposed in TAW ( 2002) provides a bias of a factor 0.670 (i.e. an average underestimate of a factor 1.5), with a RMSE = 0.747; EurOtop (2018) provides a bias of a factor 0.111 (i.e. an average underestimate of a factor 8.9), with a RMSE = 1.2021) also found that using OpenFOAM led to an overestimation of wave overtopping discharges.Fig. 9 displays the non-dimensional experimental and numerical overtopping discharges as a function of the berm width B for various wave conditions.The trends observed in the physical model tests include that a reduced amount of overtopping was found for wider berms (all lines in Fig. 9) and that an increased amount of overtopping was found for conditions with a low wave steepness (green lines above red lines and blue lines above black lines in Fig. 9).As indicated in Fig. 9 these trends were also observed in the numerical model.Also, the gradients of the computed lines (open symbols in Fig. 9) match reasonably well with the gradients of the lines from the experiments (filled symbols in Fig. 9).
Based on the results of the validation cases it can be concluded that the wave conditions that reach the structures are reasonably well described by the numerical model, although the computed wave steepness is generally lower than the measured wave steepness.
The overtopping discharges show that for the validation cases discharges are on average overestimated by the numerical model.From the validation it appears that trends that have been observed in the validation cases match with trends observed in the physical model tests.Therefore, it is concluded that it is likely that the numerical model cannot be used to accurately quantify the discharges (although more accurate than some existing methods), but the model can be used to analyse trends of various parameters.Note that reasonably accurate computations of mean overtopping discharges do not necessarily mean that individual overtopping events are all modelled with similar accuracy; reasonably accurate mean discharges can also be the result of a combination of overestimated and underestimated overtopping events.

Sensitivity analysis
A sensitivity analysis was conducted to gain more insight into the influence of certain parameters on the outcome of the numerical model, both in the water surface elevations and on the overtopping discharges.First, the mesh convergence was examined in detail to ensure the water surface elevations are independent of the grid size.Second, porous media parameters were varied (α F and β F ). Finally, a more detailed analysis on turbulence modelling is undertaken.

Mesh
The influence of the grid size on the wave propagation was analysed for a case with a crest wall and a berm with a width of B = 0.50 m.Three different mesh sizes at the water level were compared, corresponding to Δx = Δy = 0.04 m, 0.02 m, 0.01 m (named coarse, medium and fine  grid).Differences of around 5% and 8% for the significant wave height and wave steepness were encountered between the fine and medium grid.The medium grid size is selected in this research since it provides enough accuracy (error of 1% and 17% for the significant wave height and wave steepness between the numerical and experimental results) without increasing much more the computational time.The effect of the grid on the overtopping discharges is discussed in Section 4.4.

Porous media parameters
The flow inside the porous media (i.e.berm, armour and core) is simulated based on the Forchheimer equation, as mentioned in Section 2.2.5.The equation relies on the empirical terms (with coefficients α F and β F ) which are usually obtained from laboratory measurements, analytical expressions or doing a best fit to the experimental data when using a numerical model.Due to the presence of non-physical numerical dissipation in numerical models, effects of non-physical numerical dissipation can be reduced by using lower physical dissipation.Therefore, it is reasonable to assume that a calibrated numerical model will lead to optimal porous media coefficients that are somewhat lower than those obtained from physical experiments.Therefore, the sensitivity of the numerical model results has been analysed for lower porous media coefficients than the standard values (α F = 1 000 and β F = 1.1).
Due to the fact that no detailed turbulence model has been added to the numerical model, the calibration of the resistance coefficients (α F and β F ) would account for effects of turbulence inside the porous layers.
For the present sensitivity analysis the coefficient α F , which is the laminar contribution, is varied between 200 and 1 000 while the coefficient β F , which is the turbulence contribution, is varied between 0.5 and 1.1.Simulations with berm widths of B = 0 m, B = 0.50 m and B = 0.75 m were performed to study the influence of the porous media coefficients.
Overall, the results (see Table 2) show no influence of the porous media coefficients on the incident wave height (differences of less than 2.0% and 1.0% for the significant wave height and wave steepness between the simulations with different combinations of porous media parameters).Nonetheless, the porous media coefficients affect the overtopping discharges (see Table 2).By lowering both coefficients more water can flow through the permeable structure, reducing the overtopping discharges and reducing, by more than half, the difference between the experimental and numerical non-dimensional overtopping discharges.The best combination to fit the experimental data is obtained using α F = 200 and β F = 0.5.Note that the influence of the stone size D n50 and porosity n are already included by using the applied parametrization in Eq. ( 12).
The sensitivity study shows that there is an influence of the porous media coefficients in the wave-structure interaction and more specifically on the overtopping discharges.Although the study reveal that α F = 200 and β F = 0.5 return simulated overtopping discharges that are closer to the experiments, the standard values for the empirical coefficients (α F = 1 000 and β F = 1.1) were used for the validation and application of the model.As mentioned, coefficients being lower than those obtained from physical model tests, is attributed to the presence of numerical dissipation in OpenFOAM, which results in the need of lower coefficients to simulate the wave damping inside the porous media.

Turbulence modelling
To assess the validity of excluding a detailed turbulence model in the numerical model, simulations with a simple constant eddy viscosity model and with more complex turbulence models such as the stabilised versions of the k − ω and k − ω SST, developed by Larsen and Fuhrman (2018), were included in the model.Chen et al. (2021) used the stabilised k − ω model with good results for prediction of wave overtopping on (impermeable) dikes.Since the conventional two equation turbulence models (Menter, 1994;Wilcox, 2006) result in large over-estimations of the turbulence levels, Larsen and Fuhrman (2018) added a stress limiter λ 2 to overcome this problem (recommended value of 0.05).
The results showed that the simple constant eddy viscosity model does not influence the wave propagation (differences of less than 3.0% and 5.0% for the significant wave height and wave steepness between simulations with different eddy viscosity values).However, it has an impact on the flow through porous layers.By adding an eddy viscosity in the order of ∼ 10 − 3 − 10 − 4 , an extra resistance is included for the water inside the porous media.This is evidenced by higher water levels inside the porous structure along the entire simulation, which in turn leads to larger overtopping volumes.
The simulations with the stabilised turbulence models show again no significant effect on the wave propagation (differences of less than 0.5% and 2.0% for the significant wave height and wave steepness between simulations with different detailed turbulence closure models) but a severe effect on the wave overtopping and on the flow through the permeable structure.It is concluded that the numerical model does not  solve properly the combination of porous media and these detailed turbulence closure models.Thus, the original settings where no turbulence model is added, and where the turbulence inside the porous media is accounted for in the resistance coefficients of the Forchheimer equation, seems the most viable approach.In this case, no physical damping is added to OpenFOAM but the presence of the numerical dissipation already provides a stable and sufficiently accurate simulation of the flow.

Influence of parameters on wave overtopping
In this section the validated numerical model is applied to analyse the influence of various parameters on the wave overtopping discharges.The variations in wave conditions and cross-sections exceed the ranges of the applied physical model tests that were used to validate the numerical model.The focus is on analysing effects of a wider range of the wave steepness, various slope angles (not varied in the physical model tests), a wider range of heights of the crest walls, a wider range of berm widths and berm levels, variations in the permeability of the berm (not varied in the physical model tests) and the addition of a crest wall with a recurved parapet.Different combinations of the above parameters were included in the simulations to give a broader overview of not only the individual effect but also the combined effects of different breakwater layouts on the overtopping discharges.The slope was varied from cot α = 1.25 to 4 to cover a large range of rubble mound breakwaters.Although no berm was included in these simulations, other influence factors were added: a crest wall with different protruding parts of the crest wall, R c -A c = 0.0 m (no protruding wall), 0.05 m and 0.10 m above armour level and different incident wave conditions.A total of 35 cases were simulated in the numerical flume.Although the rear-side is expected not to affect the overtopping discharges (defined at the crest), here the structures were modelled with the same rear side slopes.

Slope angle
Fig. 11 shows the non-dimensional overtopping discharges as function of the slope angle (cot α) and surf-similarity parameter (i.e.Iribarren number) ξ m-1,0 .The left panels show the results of the computations without a protruding crest wall (and without a berm) but for varying wave steepness.The right panels show the results for one wave steepness but for varying heights of the protruding part of the crest wall (and without a berm).What stands out in this figure is the dependency between the slope angle and the overtopping discharge.The gentler the slope, the lower the overtopping discharges.This trend was captured regardless of the wave steepness and the presence of a crest wall.
The results indicate that not only the slope angle strongly affects the discharges, also the wave steepness strongly affects the discharges.These strong dependencies are found for geometries with (right panels in Fig. 11) and without (left panels in Fig. 11) a protruding crest wall.
The lower left panel of Fig. 11 (i.e.computations without a protruding wall and without a berm) indicates that the surf-similarity parameter can account for some effects of the slope and wave steepness but since not all data follow the same trend, there is an additional influence of the wave steepness and/or the slope angle, even for structures without a crest wall and without a berm (lower-left panel in Fig. 11).
The combined influence of slope angle and a crest wall obviously leads to lower overtopping discharges for higher crest walls.There is a steady decline in overtopping for gentler slopes, no matter the height of the crest wall; see right panels of Fig. 11.The influence of the slope angle  M. Irías Mata and M.R.A. van Gent will be discussed further in Section 5.

Berm
The influence of the berm on the overtopping discharges at breakwaters is analysed by varying the level, width and permeability of the berm.Fig. 12 exemplifies the four different configurations used to study the influence of each berm parameter on the overtopping discharges.First, the berm level was varied from 0.6 m to 1.0 m above the seabed (1.0 m corresponds to a berm at the crest of the armour) using for each of the simulations, a constant berm width of B = 0.5 m and a water depth of h = 0.75 m.Secondly, the berm width was varied with values of B = 0 m (no berm), 0.25 m, 0.50 m, 0.75 m and 1.0 m.A constant berm level (i.e.0.75 m above the seabed, which corresponds to a berm at the still water level) and a water depth of h = 0.75 m was applied.The wave steepness was varied from s m-1,0 = 0.011 to 0.039.Finally, the permeability of the berm was studied by including three different berm material configurations.The basic configuration consists of a homogenous berm as applied in the physical model tests with a stone diameter of D n50 = 0.038 m.Two additional configurations of the berm are computed.Both configurations include two different for the berm.For the outer section the same armour material as the armour layer (D n50 = 0.038 m) was used and for the inner section, the same material from the core (D n50 = 0.016 m).They differ in the fact that one represents the way  a new rubble mound with a berm is constructed while the second relates to an addition of a berm to an existing breakwater.In all the previous simulations, the material of the berm was kept as the armour material (as in the physical model tests).Now, by including these two materials into the berm, a berm with lower permeability is simulated in the numerical model.
A total of 90 simulations were included to analyse the influence of the different berm parameters on the overtopping discharges.
The effect of the berm level on the overtopping discharges is illustrated in Fig. 13.The left panel shows results for geometries without a protruding crest wall and the right panel shows results for geometries with a protruding crest wall.Both figures show that when the berm is located below still water level (i.e.B L /A c >1), there is hardly any influence on the amount of overtopping.On the other hand, above the still water level, the results reveal a sharp dependency of overtopping discharges on the berm level.The higher the berm level, the larger the reduction in overtopping discharges.In contrast to the influence of berms at dikes, the lowest overtopping discharges for breakwaters occur when the berm level is located at the same level as the crest or the same level as the armour in front of the crest wall.This effect was already noticed in the physical model tests.The numerical model indicates that the reduction due to the level of the berm depends strongly on the wave steepness; the higher the wave steepness the larger the reduction is for higher berms (i.e.low values of B L /A c ).This was also observed in the physical model tests.
Fig. 14 shows that the berm width has an influence on the overtopping discharges.This was found for structures without a crest wall (left panel of Fig. 14) and for structures with a crest wall (right panel of Fig. 14).The larger the berm width, the larger the reduction in the overtopping discharges.From these results, it also stands out that the decline in overtopping discharges is stronger for steeper waves, as also observed in the physical model tests.
The last berm parameter accounted for in the numerical model simulations is the berm permeability.The time series of wave overtopping discharges of a homogeneous berm with armour material and a non-homogeneous berm with armour material for the outer side and smaller core material for the inner side are compared in Fig. 15.On the left panel, the non-homogeneous berm represents the way a new rubble mound (with a berm) is constructed while in the right panel it represents the addition of a berm into an existing breakwater.Differences of less than 10% are encountered between the homogeneous and nonhomogeneous berm simulations.These results reveal that within the varied ranges, the influence of the berm permeability (i.e. with or without core material in the berm) on the wave overtopping discharges at breakwaters is negligible.

Crest wall
The influence of a crest wall on top of a breakwater on the amount of overtopping was tested by changing the height of the protruding part of the crest wall from non-protruding (R c -A c = 0) to 0.05 m, 0.10 m and 0.15 m above armour level.In all cases, the armour level was kept at the same level (A c = 0.25 m).
Fig. 16 shows the non-dimensional overtopping discharge as function of the protruding part of the crest wall for a condition with the high waves and a wave steepness of s m-1,0 = 0.019, with and without a berm with a width of B = 0.5 m at the still water level.What can be clearly seen in this figure is the rapid decrease in wave overtopping due to the presence of a crest wall.From no protruding wall to a protruding part that reaches half of the crest level, (R c -A c )/R c = 0.5, the overtopping discharge reduces about two orders of magnitude.Adding a berm to the breakwater cross-section does not impact the general trend in overtopping discharges due to the presence of a wall.

Recurved parapet
A recurved parapet can reduce the wave overtopping discharge, especially for relatively low overtopping discharges (for the influence of recurved parapets see for instance Oh et al., 2018 andMolines et al., 2019).Thus, to study the influence of the recurved parapet on a crest wall on a rubble mound breakwater, an additional configuration with a recurved parapet was tested in the numerical model.In Fig. 17, details of the parapet are presented.To obtain a proper representation of the recurved parapet, the mesh around the crest wall was further refined (see bottom right panel in Fig. 3).
During the numerical simulations, a grid dependency affecting the computed overtopping discharges was noticed.Due to this, to compare the effect of a crest wall with and without parapet on the overtopping volumes, the same level of refinement around the wall was applied also for the simulations without a recurved parapet (simulations with cell sizes around the crest wall of 2.5 mm and 1.25 mm instead of cell sizes of 10 mm).This level of refinement only applies for the results shown in this sub-section.
Fig. 18 shows the influence of the recurved parapet on the overtopping discharges.For large overtopping discharges, there is no reduction due to the parapet.It appears that once the space in front of the crest wall is filled with water, the rest of the wave easily overtops the crest wall.In contrast, for low discharges, the parapet does slightly reduce the amount of water passing over the structure.For the simulated cases, reductions of 16% and 32% were obtained by adding the recurved parapet.Since for the highest freeboard (R c /H m0 = 2.5) the computation with a parapet resulted in zero overtopping, while the computation without a parapet resulted in some overtopping, the reduction due to the parapet becomes infinitely high.To illustrate the difference between the computation with and without a parapet in Fig. 18, the computation with a parapet is set at a value lower than 10 − 7 for R c /H m0 = 2.5, although the actual value was zero.
Further investigation must be conducted to verify the dependency between the recurved parapet and the amount of overtopping.Nonetheless, these simulations provide insight into the effect of a recurved parapet for different wave heights, which is that the recurved parapet on a crest wall of a rubble mound breakwater is only effective for rather low overtopping discharges (similar to numerical simulations by Molines et al., 2019).

Discussion
The influence on wave overtopping discharges caused by different wave conditions and by different cross-section characteristics of the breakwater was studied in the numerical flume for wider ranges than the physical model tests on which Eqs. ( 1)-( 5) are based.First the parameters that are present in the empirical expressions are discussed, i.e. the crest wall and the berm.Thereafter, the parameter that is not present in the empirical expressions but appears to be important, i.e. the slope angle, is discussed.

Crest wall
The empirical expressions, in particular Eq. ( 1) in combination with the influence factor for crest walls as show in Eq. ( 3) were based on physical model tests within the range 0 ≤ (R c -A c )/R c ≤ 0.35 (thus a crest wall that covers about one third of the total crest elevation) while in the numerical model computations, crest wall heights up to (R c -A c )/R c = 0.5 were computed (thus a crest wall that covers half of the total crest height).
Fig. 19 shows the computational results as also shown in Fig. 16, but  now in combination with the empirical expressions (Eqs.( 1)-( 4)).The figure shows that for structures with and without a berm, the influence of the crest wall is adequately taken into account by the empirical expression (mainly Eq. ( 3)), not only qualitatively (the trends are the same) but also quantitatively (the values are very similar).This indicates that the empirical expressions can accurately be applied up to crest wall heights of up to (R c -A c )/R c = 0.5.

Berm
The empirical expressions, in particular Eq. ( 1) in combination with the influence factor for berms walls as show in Eq. ( 4) were based on physical model tests within a berm width between 0 and 0.75 m.In the numerical model computations also a berm width of B = 1.0 m is included (leading to B/H m0 < 5.85).Fig. 20 shows the computed and calculated non-dimensional discharges as function of the nondimensional berm width B/L m-1,0 .On average the computations overestimate the discharges, but the computations confirm that the wave steepness affects the influence of the berm width.However, for wide berms in combination with the steepest waves (i.e.light blue curves and symbols), the numerical model results indicate that Eq. ( 4) overestimates the reduction due to very wide berms, more than the overestimates for other berm widths; deviations between numerical model results and the empirical expression increase significantly for wide berms in combination with steep waves.This indicates that the range of validity of the empirical expression cannot be extended to very wide berms for applications of steep waves; this extension could lead to overestimates of the reductive effects of berms.
In the model tests the berm level was varied between the level of the armour at the crest (level of armour in front of a crest wall, if present), and a level of B L /A c ≤ 1.33.In the numerical model computations, also submerged berms with an even lower level were included, 0 ≤ B L /A c ≤ 1.6.The numerical model results (Fig. 13) show that for such low berm levels the effect of the berm is close to negligible.The empirical expressions also show a very limited effect on the discharges for low submerged berms (Fig. 21).Furthermore, the numerical computations confirm that there is a rather strong dependency of the wave steepness on the influence of the berm level, where the effects of the berm level are relatively strong for conditions with a high wave steepness.The numerical model computations show a strong decrease in discharges for high berms (close to the   level of the armour at the crest B L = 0).However, such strong dependency of the level of the berm is less pronounced in the numerical model results for conditions with a lower wave steepness, less than the empirical expressions indicate (such a sharp decrease close to B L = 0 is not present in the numerical model results).
Thus, the numerical model results confirm important effects of the wave steepness on the effects of the berm width and the berm level, as also shown in the physical model tests and empirical expressions.Whereas the numerical model computations indicate that the validity of the empirical expressions can be extended to lower berm levels than tested, the numerical model results indicate that the validity of the empirical expressions cannot be extended to very wide berms in combination with steep waves.

Slope angle
The numerical model results indicate that overtopping discharges increase for steeper slopes.For overtopping discharges in the relevant range (i.e.q/(gH m0 3 ) 0.5 > 10 − 6 ) the results indicate that the dependency is close to linear on a log-linear scale as shown in Fig. 11.A relatively simple way to incorporate the dependency of the slope in Eq. ( 1) would be as follows: For the slope angle cot α = 2 for which Eq. ( 1). was derived Eq. ( 1) and Eq. ( 15) result in the same discharges.However, physical reasoning on how an expression to estimate wave overtopping discharges could be obtained, can lead to a more complex expression.The following physical reasoning is applied here: • Wave overtopping discharges depend on the differences between the crest level and the level the wave run-up would have reached if the slope would be extended such that no overtopping would occur (i.e. the fictitious wave run-up level), see also Battjes (1974), Van Gent (2002) and Etemad-Shahidi et al. (2022).• The discharge decreases exponentially for higher crest levels, or for above mentioned difference between the crest level and the fictitious wave run-up level.• Except for the fictitious wave run-up level affecting the discharges, the overtopping volumes in a wave tongue increase for longer waves since the volumes in the crest of a longer wave are larger than the volumes in a shorter wave (see left panel of Fig. 22).This volume in a wave crest is linear-proportional to the wave length.Since the discharge is not expressed in volumes but in volumes per second, the effect of a larger volume in the crest of a longer wave can be expressed by the wave steepness to the power 0.5.• The (fictitious) wave tongue that would reach a certain level above the crest level is longer for a gentle slope than for a steep slope, see the right panel of Fig. 22.This makes it reasonable to assume that the volume in the (fictitious) wave tongue above the crest level is linearly proportional to the slope angle, thus resulting in a linear dependency of the discharge on the slope angle.• Influence factors can be accounted for as influence factors on the (fictitious) wave run-up levels, so should then be present in the exponential part of the expression.Above physical reasoning would result in the following wave overtopping expression, if the (fictitious) wave run-up level exceeded by 2% of the incident waves is used as a characteristic measure for the (fictitious) wave run-up: Applying Eq. ( 16) would require the development and calibration of influence factors for effects like roughness, oblique waves, berms, and crest walls to be derived based on such a concept of fictitious wave runup levels and the use of an accurate expression for wave run-up.Especially for structures with a permeable armour layer, measurements of wave run-up levels are complicated and therefore derived expression based on physical model tests are likely to go together with assumptions (like at which level above the amour the surface elevation in the wave tongue is defined) and inaccuracies in the derived expressions based on physical model tests.Alternatively, a conceptual approach to estimate (fictitious) wave run-up levels can be applied while coefficients in such wave run-up expression are calibrated based on measured overtopping discharges.Based on similarity with impermeable slopes it is reasonable to express the (fictitious) wave run-up level on permeable structures as a function of the surf-similarity parameter (i.e.Iribarren number) and linear proportional to the wave height.This leads to: where γ denotes the combined influence factors and f (γ, ξ m-1,0 ) denotes an expression for the non-dimensional (fictitious) wave run-up level.Instead of applying Eq. ( 17) and develop and calibrate influence factors for such an expression, influence factors and (a simplification of) the dependency on the surf-similarity parameter can be accounted in a somewhat different way as follows: Note that due to the presence of the surf-similarity parameter in the exponential part, the optimal powers of the wave steepness s m-1,0 and slopes angle cot α in Eq. ( 15) cannot be the same as in Eq. ( 18) (i.e.powers − 0.5 and 1 in Eq. 18 and − 1 and − 1 in Eq. ( 15) respectively).
The numerical model computations as discussed before show bias, and somewhat different dependencies on the berm width and berm level, compared to the results of the physical model tests.Therefore, calibrating coefficients c 1 , c 2 and c 3 based on the numerical model results can result in inaccuracies of the predicted wave overtopping discharges using the calibrated coefficients.Nevertheless, the coefficients can be calibrated based on the numerical model results to summarize the numerical model results, rather than to predict the actual overtopping discharges.The numerical model results can be summarized as follows: The left panel of Fig. 23 shows the comparison between the discharges computed with the numerical model and those calculated using Eq. ( 19) using Eqs.( 2)-(4) to account for the various influence factors.For larger discharges the comparison is rather good while for smaller discharges the spreading is large.
Using the tests by Van Gent et al. (2022) to calibrate the coefficients in Eq. ( 18), results in the same values for the coefficients c 1 , c 2 and c 3 .However, the wave steepness in Eq. ( 18) and Eq. ( 19) is present to a different power than in Eq. (1) (− 0.5 instead of − 1) and the wave steepness is also present via the surf-similarity parameter, while the influence factor for the berm also depends on the wave steepness.Therefore, the influence factor for the berm, that was derived based on Eq. ( 1), needs to be re-calibrated for applications in combination with Eq. ( 19).It appears that the original coefficients in Eq. ( 4) are not the optimal values if the expression is applied in combination with Eq. ( 19).The re-calibrated expression to account for the berm for application in   19)).Right panel: Measured discharges versus the empirical expression (Eq.( 19)).combination with Eq. ( 19) reads: The right panel of Fig. 23 shows the physical model test results by Van Gent et al. (2022) versus Eq. ( 19) in combination with Eqs.(2), ( 3) and ( 20) to account for roughness, a crest wall, and a berm, respectively.The RMSE between the measured and calculated non-dimensional discharges becomes 0.3333, which is somewhat larger (but still a rather low value) than the RMSE as reported by Van Gent et al. (2022;Fig. 11) using Eq.(1) (RMSE = 0.2038).Note that the advantage of Eq. ( 19) is that, unlike Eq. ( 1), the influence of the slope angle is accounted for in a way that is physically sound, despite the fact that for Eq.(15) (i.e.equal to Eq. ( 1) for cot α = 2) the RMSE value is somewhat lower than for Eq. ( 19), if compared to data from the applied physical model tests (RMSE = 0.2038 versus 0.3333).

Conclusions and recommendations
Numerical model computations have been performed for various geometries of rock-armoured rubble mound breakwaters after validation based on physical model tests.After the validation of the numerical model, the sensitivity of various parameters has been studied for a large variety of geometries, including geometries that have not been tested in the physical model.Although the numerical model results show deviations from measured discharges, the model provided valuable insights into the dependency of wave overtopping discharges on various structural parameters.The numerical modelling provided the following insights: • Wave steepness: The numerical model results confirm that wave overtopping at rubble mound breakwaters depends on the wave steepness, thus also for wave loading that can be characterised as non-breaking waves.• Crest wall: The numerical model results indicate that the earlier developed influence factor for crest walls is also valid for larger crest walls with protruding parts that cover half of the total crest height.• Recurved parapet: The numerical model results indicate that the effect of a recurved parapet on a crest wall of a rubble mound breakwater is small to negligible, except for very low overtopping discharges for which a small but noticeable effect was observed.• Berm: The numerical model results confirm that wave overtopping at rubble mound structures with a berm depends on the berm width, the level of the berm, and the wave steepness.The size of the core material inside the berm appears not to affect the overtopping discharge within the range of evaluated geometries.The numerical model results indicate that the earlier developed influence factor overestimates the effect of the width of the berm for very wide berms (wider than tested) in combination with steep waves.• Slope angle: The numerical model results show that wave overtopping at rubble mound breakwaters depends on the slope angle, thus also for wave loading that can be characterised as non-breaking waves.This effect is so large that it cannot be neglected (i.e. one or two orders of magnitude difference between structures with a slope of 1:1.5 and 1:4).The effect of the slope angle has been incorporated in an existing guideline for estimates of wave overtopping discharges at rubble mound structures (Eq.( 15)) and in a guideline developed based on physical reasoning (Eq.( 19)).The effects of the slope need to be verified based on physical model tests in which the slope angle is varied.
The numerical model results indicated that the model provides valuable insights and reproduces important trends in the dependency of wave overtopping discharges on various parameters.However, the actual overtopping values can differ from the measured values.
Therefore, it is recommended to further improve the accuracy of the numerical model as well as to enable studying 3D effects on wave overtopping numerically with a fast and accurate model.
Both the numerical model computations and the corresponding physical model tests have been derived for conditions without severe wave breaking on the foreshore.It is recommended to study the influence of wave breaking on the wave overtopping at rubble mound breakwaters.Furthermore, it is recommended to study the influence of crest walls, berms and slope angles on volumes per overtopping wave, on percentages of overtopping waves, and on flow velocities and the flow depth during overtopping events.
Furthermore, it is advised to study the influence of berms, crest walls, a shallow foreshore and offshore low-crested structures in more detail, since these can be effective measures for the adaptation of existing rubble mound structure due to sea level rise and more severe wave loading at the structure.
Fig. 1.Cross-sections of structure configurations used for validation of the model.
The relaxation zones work by weighing the computed and target solutions by means of a relaxation function α R (χ R ).The potential flow solver OCW3D provides the target solution Φ target and OpenFOAM solves for the computed solution Φ computed .At the end, the filter velocity u f and the indicator function F are updated at each time step.
and Jensen et al. (2014) used a

Fig. 2 .
Fig. 2. Layout of the numerical wave flume, with sections in which OceanWave3D and OpenFOAM® are used (measures in meters).

Fig. 3 .
Fig. 3. Basic mesh used in the OpenFOAM® numerical model (upper panel) and details of mesh around crest wall and recurved parapet.
397 and the Neural Network byVan Gent et al. (2007) provides a bias of a factor 0.42 (i.e. an average underestimate of a factor 2.4), with a RMSE = 0.518.Thus, the numerical model provides estimates that show bias within the range of other existing prediction methods but a smaller spreading than existing prediction methods (RMSE value for the numerical model is smaller).The empirical method byVan Gent et al. (2022) that is derived based on physical model tests including the mentioned 16 cases, gives a bias of a factor 1.02 and a RMSE = 0.177.Note thatChen et al. (

Fig. 9 .
Fig. 9. Computed and measured wave overtopping discharges as function of the berm width B. Filled symbols represent experimental results while nonfilled symbols represent numerical results.

Fig. 10
Fig. 10 illustrates the configuration of the breakwater used to study the influence of the seaward slope on the overtopping discharges.

Fig. 10 .
Fig. 10.Structure configurations to analyse the influence of the slope angle (1.25 ≤ cot α ≤ 4) in combination with various heights of the crest wall.

Fig. 11 .
Fig. 11.Non-dimensional overtopping discharges as function of the slope angle (cot α) and ξ m-1,0 .Left panels for varying wave steepness at structures without a wall.Right panels for varying wall heights for one wave steepness (s-mid).

Fig. 12 .
Fig. 12. Structure configurations to analyse the influence of the berm level, berm width and permeability of the berm.Upper left panel: Berm level.Upper right panel: Berm width.Middle panel and lower panel: Permeability variations.

Fig. 13 .
Fig. 13.Non-dimensional overtopping discharges as function of the level of the berm B L /A c ; left panel for structures without a protruding crest wall and right panel for structures a protruding crest wall.

Fig. 14 .
Fig. 14.Non-dimensional overtopping discharges as function of the berm width B/L m-1,0 ; left panel for structures without a protruding crest wall and right panel for structures with protruding crest wall.

Fig. 15 .
Fig. 15.Timeseries of wave overtopping discharges for a homogeneous berm and a berm with core material underneath (lower permeability).Left panel: new rubble mound is constructed (with a berm).Right panel: addition of a berm to an existing breakwater.

Fig. 16 .Fig. 17 .
Fig. 16.Non-dimensional overtopping discharges as function of the protruding part of the crest wall (R c -A c )/R c .

Fig. 18 .
Fig. 18.Non dimensional overtopping discharges with and without recurved parapet on the crest wall for different wave heights R c /H m0 .

Fig. 19 .
Fig. 19.Measured and calculated non-dimensional overtopping discharges as function of the protruding part of the crest wall (R c -A c )/R c for configurations with and without a berm.

Fig. 20 .
Fig. 20.Measured and calculated non-dimensional overtopping discharges as function of the non-dimensional berm width (B/L m-1,0 ), without a crest wall (left panel) and with a crest wall (right panel).

Fig. 21 .
Fig. 21.Measured and calculated non-dimensional overtopping discharges as function of the non-dimensional berm level (B L /A c ), without a crest wall (left panel) and with a crest wall (right panel).

Fig. 22 .
Fig. 22. Illustration of larger volumes in longer waves (left panel) and larger volumes in longer wave tongues at more gentle slopes (right panel).

Table 1
Mesh resolution for different regions.

Table 2
Wave characteristics and wave overtopping bias for different combinations of porous media parameters.