Numerical modelling of wave overtopping at dikes using OpenFOAM ®

Accurate calculation of wave overtopping is important for determining the required crest height and geometry of a dike. Berms and roughness elements are widely used to reduce the average overtopping discharge at dikes while the reductive effects of berm and roughness are still not fully understood. Several empirical formulae are available to predict the overtopping rate at coastal structures. However, the extrapolation of these empirical formulae is not always applicable for complex structures (e.g. a dike that has a berm and/or roughness elements on the waterside slopes) or wave conditions that are outside the applicability of the empirical predictors. A 2D numerical model based on OpenFOAM ® is set up in this study for predicting wave overtopping at dikes that have complex configurations with berms and roughness elements. The validation results show that this OpenFOAM ® model is capable of reproducing the incident waves accurately and predicting the wave overtopping discharge with good accuracy. Subsequently, the numerical model is applied to study the reductive influence of a berm and protruding blocks on the mean overtopping discharge at dikes. The roughness of protruding blocks is incorporated by explicitly modelling the protrusions using refined mesh. The model shows reasonable behaviour of the reduction of wave overtopping influenced by a berm and roughness. This indicates the capabilities of the numerical model in the design and safety assessment of dikes.


Introduction
An accurate prediction of wave overtopping discharge is an important requirement for the design and safety assessment of dikes.The average overtopping discharge at coastal structures should be below acceptable limits under specified wave conditions and water levels in the design of dikes to ensure the stability of coastal structures thereby protecting infrastructure and people in the hinterland (CIRIA et al., 2007).In order to satisfy these criteria, berms and roughness elements are often applied at the seaward side of dikes in areas where it is not feasible, undesirable or uneconomic to raise the crest level of dikes.However, the effects on the average overtopping discharge are still not fully understood.Several empirical formulae are available to estimate the reductive influence on wave overtopping at dikes.The TAW (2002) and EurOtop (2018) manuals provide empirical methods for calculating the influence factors of berms and roughness.Nevertheless, Chen et al. (2020a;2020b) made clear that both the TAW (2002) and EurOtop (2018) equations showed weak performance for estimating these two influence factors by comparing them to experimental data.Therefore, Chen et al. (2020b) developed empirical equations for berm and roughness influence factors, based on the analysis of experiments.However, the validity of empirical methods that are mainly based on physical model tests within a limited range of conditions and configurations is not always certain for other conditions and configurations.For example, the validity of the empirical equations developed by Chen et al. (2020b) is uncertain for dike configurations or wave conditions outside the tested ranges.Additionally, new detailed physical model tests can be costly and it can be time consuming to prepare a new physical model setup.
With the development of numerical methods and computational resources, numerical modelling has become an important complementary tool with experiments for modelling interaction of waves with coastal structures.There are different types of numerical models depending on the governing equations and solving technique.A computationally efficient approach to model wave interaction with coastal structures is based on the nonlinear shallow water equations (NSWE) (see for instance Kobayashi and Wurjanto, 1989;Van Gent, 1995, 2000;Ma et al., 2012), such as the SWASH model (e.g.Zijlema et al., 2011;Suzuki et al., 2017).Models based on NSWE can be very efficient, providing the possibility of simulating wave trains of 1000 waves rapidly (Losada et al., 2008).However, one restriction of the use of NSWE is associated with the difficulty of dealing with complex geometries of coastal structures (Suzuki et al., 2017).There exist many numerical models based on the Reynolds-Averaged Navier-Stokes (RANS) equations.In the numerical models based on RANS equations, there exist two main free surface tracking techniques, i.e. volume of fluid (VOF) and particle methods like Smooth Particle Hydrodynamic (SPH).SPH is capable of tracking large deformations of the free surface accurately and has been applied to simulate wave overtopping on impermeable coastal structures (Didier andNeves, 2009 andAkbari, 2017).Nevertheless, the SPH method is expensive in terms of computational time since it requires a large number of particles and small time steps (Ο(10 − 5 s)) to obtain enough accuracy (Klapp et al., 2016).Models based on the VOF method are capable of dealing with large free surface deformation and are less computationally expensive than the models based on SPH, such as the models by Van Gent (1995), Losada et al. (2008), Fang et al. (2010) and Higuera et al. (2013).
OpenFOAM® is an open-source computational fluid dynamics framework that allows users to run tasks in parallel on multiple processors.The open-source nature often leads to useful libraries and toolboxes that are freely shared in the public domain (Davidson et al., 2015).Additionally, OpenFOAM® has recently been developed with the capability of modelling interaction of waves and structures and it appears to become increasingly popular for coastal engineering (Jacobsen et al., 2015).Therefore, OpenFOAM® is adopted in this study to model wave overtopping at dikes that have complex configurations.
Several studies have been conducted on numerical modelling of wave overtopping at coastal structures making use of OpenFOAM®.Higuera et al. (2013) applied a three-dimensional Navier-Stokes (NS) solver called IHFOAM built on OpenFOAM® to investigate the interaction between waves and a high mound breakwater.A 3D model is computationally expensive and therefore only a total of 40 s is simulated for wave overtopping with irregular waves, which is insufficient for estimations of average overtopping discharge.Higuera et al. (2013) also mentioned that long simulations should be carried out to obtain reliable results of wave overtopping.Jensen et al. (2014) simulated wave overtopping over a smooth straight impermeable structure and a porous breakwater that has a straight waterside slope.In their 2D numerical model, waves were generated and absorbed using waves2Foam package (Jacobsen et al., 2012) which is a toolbox based on OpenFOAM® that applies the relaxation zone technique to generate and absorb free surface water waves.The numerical model in Jensen et al. (2014) was found to give a good agreement between the simulated and measured average overtopping discharges.Since not too much wave breaking was observed during the physical model tests that were used to validate the numerical model in Jensen et al. (2014), no detailed turbulence model was applied in the model.Patil (2019) also used OpenFOAM® with waves2Foam to model wave overtopping over a smooth straight slope and there was a reasonable agreement between the modelled and experimental results while the overtopping discharge was slightly underestimated by the numerical model.Previous research mainly focused on the simulation of the overtopping process at structures that have simple configurations, such as structures with straight and/or smooth seaward slopes.Validation of the capability of OpenFOAM® in modelling overtopping discharge with long time series of irregular waves at dikes that have complex configurations with a berm and roughness elements is to the authors' knowledge still not available.
The objectives of this study are to explore the capability of the OpenFOAM® model for accurately estimating wave overtopping discharges at impermeable dikes that have a berm and roughness elements at the seaward side as well as to explore the applicability of the numerical model to investigate the influence of several configurations of berms and roughness on average overtopping discharge at dikes.The OpenFOAM® version 1812 is used in this study.The existing wave-s2Foam package is applied to generate waves and the solver wave-IsoFoam is applied to solve the numerical model.Physical model tests from Chen et al. (2020b) are used to validate the OpenFOAM® model for different configurations including a smooth straight slope (SS), a smooth slope with a berm (SB), a straight slope covered by protruding blocks (PbS), a bermed slope covered by protruding blocks (PbB) and a bermed slope with only the upper slope covered by protruding blocks (PbB_up).Following the validation, the berm width and berm level relative to the still water level (SWL) are varied in the numerical model to study the reductive influence of a berm on the average overtopping discharge at dikes for conditions and configurations outside the tested ranges.Additionally, the coverage length of protruding blocks on the upper slope is varied to investigate how the coverage length of roughness elements affects the total roughness factor.
This paper is organised as follows: A description of the methodology is given in Section 2. In Section 3, the validation of the OpenFOAM® model is presented.Section 4 focuses on the application of the validated model to investigate the berm influence on the overtopping discharge.Model results of roughness influence are presented in Section 5. Applicability, limitations and efficiency of the numerical model are further discussed in Section 6. Section 7 summarises the main conclusions of this study.

Methodology
In this section, the physical model tests that are used later to validate the numerical model are introduced first.Following that, a detailed description of the numerical model setup is given.the wave overtopping discharge at dikes.The selected physical model tests include five campaigns over five dike configurations (Fig. 1 and 2): a smooth straight slope (SS), a smooth slope with a berm (SB), a straight slope covered by protruding blocks (PbS), a bermed slope covered by protruding blocks (PbB) and a bermed slope with only the upper slope covered by protruding blocks (PbB_up).The protruding block revetment is considered as impermeable.A slope of 1:3 was used for the straight slope and for the down and upper slopes of structures that have a berm.The core of all the structures was impermeable and was made of concrete.The berm position and berm width (0.2 m), see Fig. 2b, remained unchanged in the physical model tests.In Fig. 2c,d & e, the size of the cubes was 5 cm × 5 cm × 5 cm.The chessboard of protruding block revetment was made by installing a concrete tile of 1 cm thick underneath the cubes.Thus, the protrusion height was 1 cm.In practice, installation of these artificial blocks is conducted mechanically.About 16-18 units forming a 1 by 1 m panel are installed at once (Capel, 2015) and Fig. 3 shows the chessboard pattern block revetment.The modelled structures were placed at a distance of 11 m from the wave board.Three  W. Chen et al. wave gauges were placed near the toe of the tested structure to measure wave conditions including the significant wave height H m0 and the wave period T m− 1,0 .Incident and reflected waves were separated making use of the method developed by Mansard and Funke (1980).The wave height H m0 , wave steepness s m− 1,0 (s m− 1,0 = 2πHm0 gT 2 m− 1,0 ), freeboard R c (which refers to the distance between the crest and the still water level, SWL), and berm level (d h ) relative to SWL were varied in the physical model tests.Note that the freeboard R c and berm level d h were varied in the ranges of (0.12 m, 0.18 m) and (− 0.03 m, 0.03 m) respectively by varying the water depth between 0.57 m and 0.63 m with the positions of crest and berm fixed.The overtopped water was collected using an overtopping box in which a wave gauge was installed to detect the water surface variation.The average overtopping discharge was measured for about 1000 waves for each test run.

Numerical setup
The numerical method in OpenFOAM® is based on a finite volume discretisation with a collocated variable arrangement on grids.In this study, the Navier-Stokes equations are solved for the two-phase flow (air and water).The interface of water and air is captured using the iso-Advector developed by Roenby et al. (2016) as suggested by Larsen et al. (2019).The isoAdvector algorithm is a VOF based interface advection method.It is applicable for arbitrary meshes aiming to keep the accuracy of geometric schemes and obtain acceptable calculation time making the geometric operations at a minimum.
The 2DV numerical wave flume was set up to mimic the experimental layout as illustrated in Fig. 4 with a total length of 17 m and a height of 1.2 m.The dike was placed 11 m away from the inlet boundary, which is the same as the distance between the dike and waveboard in the experiments of Chen et al. (2020b).Three wave gauges used for incident wave analysis were defined at the same locations as those in the laboratory experiments.The modelled impermeable structures and bottom were set as non-slip conditions.A constant pressure was enforced at the atmosphere boundary which allowed the air to flow in and out while the water could only flow out.Relaxation zones were applied at the inlet and outlet boundaries to generate and absorb waves, which is described in more detail in Section 2.2.2.

Mesh
For smooth slopes, including SS and SB, the numerical mesh was created by using blockMesh and snappyHexMesh.Both blockMesh and snappyHexMesh are mesh-generation tools implemented in the Open-FOAM® package.The base mesh was created using blockMesh taking the structure into account.The mesh from the inlet boundary to the toe of the dike was orthogonal and conformal with the cell size of 0.02 m both in X direction (horizontal) and in Y direction (vertical).The overtopping discharge is closely related to the wave run-up and Wroniszewski et al. (2014) suggested that the quadrilateral grids parallel with the slope surface could improve the accuracy of the simulation of the wave run-up on a slope.Thus, this type of mesh near the slope surface was used in this study as shown in Fig. 5a snappyHexMesh was then applied to refine the mesh around the free surface in front of the modelled structure with one refinement level, resulting in the cell size of 0.01 m in X and Y.There are nearly 12 cells to resolve per wave height, which is sufficient for modelling the wave propagation (Larsen et al., 2019).
The protruding blocks installed on the waterside slope were in chessboard pattern as shown in Fig. 2c which is impossible to model in a 2DV numerical model.Therefore, the protruding blocks were simplified by only modelling the protrusions shown in Fig. 5b.The mesh for the slopes covered by protruding blocks was created by using GMSH and snappyHexMesh.GMSH is a third-party software that is capable of generating smooth boundary conforming mesh near the wall with geometric transitions.Protrusions of block revetment were modelled explicitly by creating fine grid cells around the protrusions.GMSH was applied to create a base mesh similar to the mesh for smooth structures taking the dike and protrusions into account.The mesh around protrusions was refined based on the base mesh using snappyHexMesh (see Fig. 5b) with the refined grid cell size of 0.0025 m in the direction perpendicular to the slope and 0.0031 m in the direction parallel with the slope.The influence of the grid size around the protrusions has been analysed for the overtopping discharge as presented in Appendix A. Further refining the grid did not result in more accurate model results of average overtopping discharge.Refinement was also applied around the free surface by using snappyHexMesh with one refinement level.The aspect ratio of the computational cells in all simulations was close to 1 as suggested by Jacobsen et al. (2012) for the simulation of wave propagation and wave breaking.Table 1 presents the mesh resolution in different regions of the numerical domain.

Wave generation
OceanWave3D and waves2Foam were used in order to generate the consistent time series of irregular waves with experiments.Ocean-Wave3D is a robust and efficient flexible-order finite difference model based on a fully nonlinear and dispersive potential flow model (Engsig-Karup et al., 2009).OceanWave3D has already been coupled into OpenFOAM®/waves2Foam through the interface described in Paulsen et al. (2014).The one-way coupling makes use of the relaxation zones provided by the waves2Foam utility.The potential flow solver provides the target solution ψ target in the relaxation zones also named coupling zones where the velocity field u and the water volume fraction α are updated at each time step according to where ψ target is the target solution in time and space given by the potential flow solver; ψ computed are the computed quantities by the Open-FOAM model; χ ∈ [0, 1] is a weighting factor.The coupling aimed to provide a fully nonlinear kinematics in the numerical model thereby improving the accuracy of modelled hydraulic loading on coastal structures (Jacobsen et al., 2018).
In this study, waves were first generated with the application of OceanWave3D by inputting the wave paddle signal files that were used in the experiments.Then, OceanWave3D provided target solution for waves2Foam from the starting time defined in the OpenFOAM model through the inlet relaxation zone which is referred to as coupling zone (shown in Fig. 4) and then the two solvers ran side-by-side.In this way the consistent time series of water surface variation with those in the physical model tests were simulated in the numerical model.The relaxation zone implemented at the inlet also absorbed reflected waves from structures.Another relaxation zone was applied at the outlet to avoid waves reflected.The length of the inlet relaxation zone was 5 m which was about one wave length in this study.Since the wave motion behind the structure was not the focus of this study, the length of the outlet relaxation zone was only 2 m to save computational time.

Turbulence modelling
Breaking waves at dikes will induce turbulence which has an effect on the flow velocity along the waterside slope surface, thereby influencing wave overtopping discharge at dikes.The turbulence was modelled by applying the stabilised k − ω turbulence model developed by Larsen and Fuhrman (2018).The conventional k− ω turbulence model (Wilcox, 2006) results in an exponential flow with finite strain, which gives rise to severe over-estimation of turbulence levels.Therefore, Larsen and Fuhrman (2018) proposed a solution to fix this problem by introducing a new stress limiter λ 2 in the calculation of eddy viscosity ν t without alteration of any fundamental closure coefficients.Herein, λ 2 defines the effective potential flow threshold.The added limiter will become active only in a region of nearly potential flow.Larsen and Fuhrman (2018) mentioned that λ 2 should be small, but also large enough to work for practical applications.λ 2 = 0.05 as suggested by Larsen and Fuhrman (2018) was used in all simulations in this study.Wall functions were applied in the boundary layer instead of resolving the boundary, which saves computational time.Application of the wall functions required that the first grid near the wall should be located in the log layer.Therefore, the mesh near the wall was refined as shown in Fig. 5a, resulting in the grid size near the wall of Δy ≈ 0.0022 m.A Nikuradse's roughness height of 0.001 m was adopted to represent the smooth concrete surface used in the experiments.All simulations were solved by using waveIsoFoam which is a solver implemented with the isoAdvection scheme in waves2Foam.

Post processing of overtopping discharge
Overtopping discharge was obtained by utilising the overtopping function in waves2Foam.A set of cell faces at the inner edge of the crest were selected to calculate the volume of overtopped water.The volume flux over the selected faces S (Jacobsen, 2017) is calculated by in which q is the volume flux [m 3 /s] and S f is the non-unit normal vector to the face.
where t start and t end are the starting time and end time of the calculation of the total overtopping volume V tot [m 3 ]; Δt [s] is the time step.

Validation of the model
In this section, the validation of the OpenFOAM® model for reproducing the incident irregular waves and predicting the overtopping discharges at impermeable dikes is presented.28 cases in total were simulated using the numerical model.The experimental wave conditions of these 28 cases are given in Table B.1 in Appendix B. The simulation time of the numerical model was set as 500 s for all simulated cases, corresponding to between 250 and 350 waves depending on the wave period.This duration of the simulation is determined by comparing the measured overtopping using partial time series with that using the whole overtopping time series.It was found that 500 s of time series gives a factor of 1-1.15 of the overtopping rates based on the entire time series.Thus, 500 s was adopted to comprise between the computational effort and the accuracy of the results.Model results are compared with the data from the physical model tests of Chen et al. (2020b).The simulated time series of free surface elevation, wave energy spectra and wave properties are first compared to the measured results obtained based on the initial 500 s from laboratory experiments to show the ability of the numerical model to reproduce the incident irregular waves.The model results of wave overtopping discharges are then compared to the measured discharges to evaluate the performance of the model for predicting overtopping discharges at dikes with different types of configurations.

Validation of incident irregular waves
Wave conditions for all 28 tests were validated by comparing the time series of free surface elevation of OpenFOAM® to those of the experiments.Only two wave conditions, i.e. case SS1 and case SS7 were selected for sake of brevity to show the ability of the OpenFOAM® model in reproducing irregular waves with different wave steepnesses.and WG3 using the method proposed by Mansard and Funke (1980).The agreement between numerical and physical model results for both cases is generally good with a root mean square error (RMSE) of 0.0184 m for case SS1 and a RMSE of 0.015 m for case SS7.Similar good agreement is found for all other simulated cases, which demonstrates that the numerical model can well reproduce the wave propagation.
An accurate reproduction of the wave energy spectra of incident waves in front of the structure is essential for an accurate prediction of the overtopping discharge at dikes.Fig. 7 shows the comparisons between modelled and experimental wave energy spectra of incident waves for case SS1 and case SS7.The experimental wave spectra are obtained based on the measured time series of incident wave surface elevation of 0-500 s corresponding to the simulation time in the numerical model.The shape of the wave energy spectra of measured incident waves from physical model is generally well predicted using the numerical model.Based on the wave energy spectra, the significant wave height H m0 and spectral wave period T m− 1,0 are calculated.Fig. 8 presents the ratio of  8b shows that the wave period T m− 1,0 is slightly overestimated by the OpenFOAM® model.This might be caused by the slight overestimation of wave period produced by OceanWave3D which provides input of wave signal for waves2Foam.In the OceanWave3D model, the reflected waves are not totally absorbed by the sponge layer and therefore interact with the incident waves resulting in larger wave period.Nevertheless, the mean absolute percentage error (MAPE) for T m− 1,0 is only 5.3%, which has limited influence on the average overtopping discharge.Thus, the difference is regarded as being acceptable.
The measured and modelled incident wave height H m0 and wave period T m− 1,0 for all simulated cases can be found in Table B.1.
Overall, the results presented in this section indicate that the OpenFOAM® model is capable of accurately reproducing the incident irregular waves in terms of wave propagation and wave properties including the significant wave height H m0 and the spectral wave period T m− 1,0 .

Validation of wave overtopping discharge
The capability of OpenFOAM® model to predict overtopping discharge at dikes is evaluated by comparing the computed cumulative overtopping discharge by the numerical model to the measured cumulative overtopping discharge from the experiments.The overtopping volume in the physical model tests was measured using a wave gauge where N is the total number of simulated cases; log(q * experiment ) is the mean value of log(q * experiment ).Bias indicates the tendency of the numerical model to overestimate or underestimate the average overtopping discharge.NSE measures the correlation between the experimental and numerical overtopping discharges and it varies between − ∞ and 1. NSE = 1 represents a perfect agreement of numerical data to the experimental data and NSE = 0 means that the numerical values are as accurate as the mean of the experimental data.NSE < 0 indicates that the experimental mean value is a better predictor than the numerical model.Therefore, the closer the NSE is to 1, the more accurate the numerical model is.
Fig. 10 shows a good agreement between the computed and measured dimensionless average overtopping discharges at smooth straight slopes within a wide range of (0.0001, 0.04).An overall slight overestimation for smooth bermed slopes might be because the air entrainment produced by the berm in experiments cannot be accounted for by the numerical model leading to less energy dissipation.For the slopes (partly) covered by protruding blocks, the average overtopping discharges are overestimated by the OpenFOAM® model, resulting in the overall Bias of 0.19.One of the causes for this overestimation might be that the turbulence generated by the protrusions is underestimated in the 2D model as the blocks are in chessboard pattern in the experiments, which cannot be modelled in a 2D model.The NSE value of 0.78 shows an overall good match of the numerical average discharge with the experimental results.It can be concluded that the OpenFOAM® model is capable of predicting the overtopping discharges at dikes with simple and complex configurations with the accuracy within a factor of 1-3 of the experimental overtopping discharges.EurOtop (2018) stated that overtopping rates estimated by empirically derived equations are within, at best, a factor of 1-3 of the actual overtopping rates.Besides, Fig. 15 in Chen et al. (2020b) showed that the empirical equations provided by TAW (2002) and EurOtop (2018) gave a factor of 1-8 of the measured overtopping rates from the experiments.Therefore, the accuracy of the numerical model in this study is equivalent to, or even better, than that of empirical equations.2002) by means of small-scale physical model tests.However, in those experiments, the width of the berm (B) was fixed at 0.2 m and the berm level relative to SWL (d h ) was varied in a very limited range of (− 0.03 m, 0.03 m), which may limit the applicability of the empirical berm equation.Therefore, the validated OpenFOAM® model is applied to generate a numerical dataset by varying the berm width and berm level in a larger range.The numerical dataset is then analysed to investigate the influence of berm width and berm level on the average overtopping discharge at dikes.

Numerical experiments
The same configuration as shown in Fig. 2b was considered in the numerical model, except that the berm width and berm level were varied.The berm width was varied in the range of 0.0 m-0.6 m (on the scale of the small-scale experiments) combined with different values of the wave steepness in the numerical model.The berm level was varied between − 0.15 m and 0.125 m by changing the position of the berm with the still water level (SWL) fixed at 0.6 m.The berm width was chosen as 0.2 m and 0.4 m. 28 cases (Table B.2) were simulated in total.Each case was simulated for around 500 s.

Model results
Model results of average overtopping discharges are presented in Table B.2. Fig. 11 shows the relationship between the relative berm width (B/H m0 ) and the ratio of the modelled dimensionless average overtopping discharge (q * ) and the modelled discharge without a berm (q * no berm ).There is a clear trend of the dimensionless mean overtopping discharge decreasing with increasing relative berm width, which is in line with the results from previous studies (e.g.Fig. 1 in Van Gent, 2013).It is worth noting that the application of a berm at a dike can significantly reduce the average overtopping discharge.Applying a berm with B/H m0 ≈ 2 can lead to 40% less than the average overtopping discharge over smooth straight slope according to the model results.For the relative berm width B/H m0 > 4, further extension of berm width may not lead to a significant reduction of the average overtopping discharge.The berm influence factor γ b is usually used to represent the reductive influence of a berm on overtopping discharge (e.g.TAW, 2002;EurOtop, 2018;Chen et al., 2020b).Values of the berm factor γ b can be obtained by solving the overtopping equation (C.1) with the overtopping discharge q substituted by the modelled average overtopping discharges.Fig. 12 shows that the berm influence factor γ b decreases as the relative berm width B/L berm (where L berm is the characteristic berm length) increases.Here, B/L berm is used since this makes it easier to compare with the empirical equations.Open markers in Fig. 12 represent some of the experimental data on which the empirical berm equation was derived.A good agreement between the experimental results and empirical equations can be seen.Note that the B/L berm from the experimental data varied due to variation in the wave height but not due to variations in the berm width since the berm width in the experiments was fixed as 0.2 m.The empirical equations show a similar trend with the model results while some differences between the model results and empirical equations can be noticed.
The numerical modelled relationship between the berm influence factor γ b and relative berm level is shown in Fig. 13 in which R u2% is the wave run-up height that is exceeded by 2% of the number of incoming where r B can be calculated using Eq.(C.3) and r dh is calculated using Eq.(C.4).The values of R 2 and RMSE for the fitting equation ( 6) are 0.86 and 0.03 respectively.Even though Eq. ( 6) is not recommended for applications without further validation, it indicates that the empirical berm equation may require recalibration of the empirical coefficient b 0 when applying to cases outside of the experimental range.

Roughness influence
It is time-consuming to change the applied locations of roughness elements in the physical model tests.Therefore, in our earlier physical model tests (Chen et al., 2020b), we only considered the reductive influence of protruding blocks on the upper slope, blocks on the upper slope and berm, as well as blocks on the entire waterside slope surface.Chen et al. (2020b) concluded that the protruding blocks on the upper slope contribute the most to the total roughness influence factor.However, it remains unclear how much influence the protruding blocks have if only applied on a part of the upper slope.The numerical simulations are performed aiming to study the influence of coverage length of the protruding blocks on the partial upper slope on the total roughness factor of the waterside slope.

Numerical experiments
Three configurations, including one protrusion, two protrusions and three protrusions on the upper slope, are modelled in the OpenFOAM® model.Fig. 15 shows the meshes around protrusions for these three configurations.The water depth was kept at 0.6 m, which also means that the freeboard and berm level remained constant in all simulations.The berm width was fixed at 0.2 m.Four wave conditions which are given in Table 2 were considered.Combining the wave condition and dike configuration results in 12 simulations.Thus, for each dike configuration, four wave conditions were applied.

Model results
Model results of four series of simulations are reported in Table B.3.As mentioned in Section 3.2, the average overtopping discharges are overestimated by the OpenFOAM® model for the configurations with protruding blocks.Fig. 16 shows the comparison between the measured and modelled dimensionless average overtopping discharges at slopes with upper slope covered by blocks (PbB_up) and slopes entirely covered by blocks (PbB).The OpenFOAM® model tends to give a larger overestimation for smaller overtopping discharge.A relationship between numerical and experimental results is found as shown in Fig. 16 by fitting the data using the least square method.To analyse the trend of the computational results of the roughness influence factor of protruding blocks, the average overtopping discharge predicted by the numerical model were modified by using the following equation where q * OpenFOAM Mod is the modified dimensionless average overtopping discharge and q * OpenFOAM represents the original overtopping discharge given by the OpenFOAM® model.
Table 3 lists the values of the roughness influence factors for all simulations including the simulations for the bermed slope covered by protruding blocks (PbB) and for the bermed slope with the upper slope    covered by protruding blocks (PbB_up).In the table, L 1 represents one protrusion on the upper slope, L 2 represents two protrusions on the upper slope, L 3 corresponds to three protrusions on the upper slope, L 4 represents the protruding blocks on the upper slope and L 5 means the protruding blocks applied on the entire slope surface.The calculated roughness influence factors for all simulations are also shown in Fig. 17 with the relative coverage length of protruding blocks 'L/ L 5 ' on the horizontal axis and the roughness factor on the vertical axis.The coverage length, L, is the effective length of applied protruding blocks along the slope as shown in Fig. 15c.The roughness elements located below 0.25R u2%, smooth under the still water level have little or no effect on the total roughness factor according to TAW (2002).Thus, for the protruding blocks on the down slope, only the blocks lying above 0.25R u2%, smooth under the still water level are taken into account to calculate the coverage length.R u2%, smooth is the wave run-up level on a smooth slope which can be calculated using Eqs.(C.5 & C.6). Six points from left to right in each scattered diagram in Fig. 16 represent the roughness factors based on the model results of smooth bermed slope, one protrusion, two protrusions, three protrusions, slopes with upper slope covered by blocks and slopes entirely covered by blocks.In general, the roughness influence factor decreases, which means that the reductive influence of roughness increases, as the coverage length of protruding blocks increases.The average contribution of protruding blocks for four wave conditions to the total roughness factor is estimated using Eq. ( 8).
where M is the number of wave conditions (M = 4); γ f− totali [− ] is the calculated overall roughness factor of the waterside slope based on the model results; γ f− pb i [− ] is the roughness factor of the bermed slope entirely covered by protruding blocks.
The application of only one protrusion can contribute 44% to the total roughness influence.It is also worth noting that the roughness factors for three protrusions are almost the same as those for the protruding blocks on the entire upper slope (four protrusions).On average, the blocks on the half upper slope (three protrusions) contribute 76% to the total roughness influence, while five additional protrusions (blocks on the entire upper slope) also lead to a contribution of 76% to the total roughness influence (Table 4).This indicates that the protruding blocks on the lower part of the upper slope do not play a role in reducing the average overtopping discharge.

Applicability and limitations of the model
This study provides a systematic validation of the OpenFOAM® model for predicting the wave overtopping discharge resulting from irregular waves at dikes with simple and complex configurations.The model validation shows that the model is capable of predicting the overtopping at smooth straight and smooth bermed slope accurately.However, it overestimates the average overtopping discharges at slopes (partly) covered by protruding blocks about twice of the measured discharge from the experiments.This overestimation might be partly caused by the underestimation of the turbulence produced by protrusions in the 2D model considering the protruding blocks were placed in a chessboard pattern in the experiments.The width of the model in Fig. 17.Comparison of roughness factors for from left to right the smooth bermed slope, one protrusion on the upper slope, two protrusions on the upper slope, three protrusions on the upper slope, protruding blocks on the entire upper slope and protruding blocks on the whole waterside slope under four wave conditions (W1, W2, W3 and W4).

Table 3
Roughness factors based on model results for all simulations, where W represents the different wave conditions and L the different protrusion lengths with L 1 for one protrusion, L 2 for two protrusions, L 3 for three protrusions, L 4 for protruding blocks on the entire upper slope and L 5 for protruding blocks on the whole waterside slope.the experiments was 1 m.Since the slope surface in a 2DV numerical model was assumed uniform in z direction, the protrusions would act as parallel ribs across the model width instead of as protruding blocks in chessboard pattern imagining that the protrusions in the 2D numerical model were extended by 1 m along z direction.The estimated roughness influence factor of ribs (optimal) is about 0.86 using the equation (5.24) in EurOtop ( 2018) while the mean roughness influence factor is around 0.72 from Chen et al. (2020b).This indicates that the ribs might lead to less reductive influence on the overtopping and therefore larger overtopping discharge than protruding blocks in chessboard pattern.Besides, the numerical model cannot take the aeration caused by protrusions into account, which could also lead to overestimation of overtopping rate.
We applied the validated OpenFOAM® model to investigate the berm and roughness influence.Results of the berm influence show that the OpenFOAM® model is capable of computing the reductive influence of a berm on the average overtopping discharges at dikes reasonably.Comparison between numerical results and the empirical berm equation indicates similar trends of the berm factor changing with the berm width and the berm level.However, there is still some difference between the model results and the empirical berm equation derived in Chen et al. (2020b), which indicates that the empirical equation might require modification for applications outside of the experimental ranges.To extend the range of application, a modified berm equation was developed based on the model results by recalibrating the empirical coefficient b 0 in the empirical berm equation.We emphasize that this new berm equation needs further validation before it can be recommended for application in practice.Nevertheless, the numerical model extended the insights into the dependencies of parameters and insights into the applicability of empirical equations.
The influence of the coverage length of protruding blocks on part of the upper slope on the total roughness factor was investigated with the OpenFOAM® model.The numerical results are in line with the findings in Chen et al. (2020b) that the roughness elements on the higher part of the waterside slope play a major role in reducing the wave overtopping (see also Capel, 2015).Model results also show that the total roughness factor is sensitive to the location of the protruding blocks applied on the half upper slope.Applying one row of protruding blocks on the top of the upper slope can effectively increase the reductive roughness influence on the overtopping discharge.Protruding blocks on the upper half of the upper slope lead to almost the same total roughness factor as those on the entire upper slope.
The model results for berm and roughness influence demonstrate the added value of the OpenFOAM® model for predicting the average overtopping discharges at dikes with complex configurations.The numerical model can be used in combination with empirical equations to predict the overtopping discharge during a functional design of a dike, since the numerical model can deal with a wider range of wave conditions and dike configurations than that in the experimental setup.Besides, when applying berm and/or roughness elements at a dike is necessary, this numerical model can provide insights into the optimal berm width and the optimal locations and coverage length of roughness elements along the waterside slope for reducing the wave overtopping discharge.For example, model results for the roughness influence presented in Section 5.2 show that the roughness elements applied on the upper half of the upper slope can lead to approximately the same reductive influence of those applied on the entire upper slope, which suggests a possibility of reducing the costs related to dike construction.
Although the numerical model is capable of predicting the discharge with a good accuracy, it still has some limitations.One limitation of the 2DV numerical model developed in this study is that it assumes the coastal structures and waves are alongshore-homogeneous and therefore the directional and directional spreading effects on wave overtopping and wave transformation cannot be taken into account in a 2DV approach.The numerical model was developed for impermeable structures considering dikes are often impermeable in the practical engineering.Therefore, this OpenFOAM® model requires adaptation for the cases that permeable armour units or permeable core material are applied at the coastal structures.

Computational efficiency of the model
In the present study, one simulation for wave overtopping at a smooth straight slope for 500 s took about 5 days parallelizing the case into 3 processors (3.6 GHz) which is quite a small number of processors.For the wave overtopping at slopes covered by protruding blocks, it took 10 days to finish one simulation of 500 s with 3 processors (3.6 GHz).The OpenFOAM® model is not the most computationally efficient model.Models based on NLSW are more computationally efficient but less accurate than the models based on RANS equations.The computational time of the numerical model developed in this study is acceptable considering the good accuracy in reproducing the incident waves and in predicting the overtopping discharge.For example, during the functional design or reinforcement of a dike, an accurate prediction of average overtopping discharge is essential.In this case, it is beneficial to have a model like this 2DV OpenFOAM® model with good accuracy and reasonable computational cost.

Conclusions
In this paper, an investigation of the capability of the 2DV Open-FOAM® model in predicting wave overtopping for impermeable structures as well as the applicability of the numerical model in simulating the berm and roughness influence is presented.
Twenty-eight tests for five configurations from the experiments in Chen et al. (2020b) have been used to validate the 2DV OpenFOAM® in order to estimate the wave overtopping.Model results show that the incident waves that are used in the experiments can be well reproduced in the numerical model with the specific settings in this study in terms of wave propagation and wave properties.The comparison between the numerical and experimental overtopping discharges demonstrates that the numerical model is capable of predicting the average overtopping discharge at simple and complex configurations within a factor of 1-3 of the experimental discharges and with a NSE of 0.78.Besides, similar trends of modelled and experimental time series of cumulative overtopping volume indicate that the OpenFOAM® model can also capture the individual overtopping events.
The validated OpenFOAM® model was then applied to study the berm and roughness influence on the average overtopping discharge.The numerical model showed reasonable behaviour of the variation of average overtopping discharge when taking the berm and roughness elements into account.The model results indicate a potential improvement of the empirical berm equation developed by Chen et al. (2020b) for the cases outside the experimental conditions.Results for roughness influence show that the roughness influence factor is very sensitive to the coverage length of roughness elements along the upper half of the upper slope.We found that protruding blocks on the upper half of the upper slope resulted in similar roughness factor with those applied on the entire upper slope.
We recommend to analyse the influence of the berm and roughness also in combination with oblique waves, with and without directional spreading.
The results of model validation and numerical experiments show potential applications of the 2DV OpenFOAM® with the specific settings used in this study in the dike design and safety assessment.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Table B.3
Simulated cases and model results of average overtopping discharge (q) for roughness influence for the four wave conditions (W1-W4) and varying coverage length of the protruding blocks on the upper slope (L1-L5).Where r B and r dh are calculated using the following equations.
Fig. 1.Smooth slopes with and without a berm tested in the experiments with SWL varied between 0.57 m and 0.63 m (unit in graphs: m).

Fig. 4 .Fig. 5 .
Fig. 4. Layout of the 2D numerical wave flume (unit in graphs: m) which is inside the domain (colored in light yellow) of OCW3D model and the OCW3D model is coupled with the OpenFOAM model through the inlet relaxation zone referred to as coupling zone.(For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) φ F [m 3 /s] is the flux of fluid across a face multiplied with the indicator function; φ ρ [m 3 /s] is the flux of fluid across a face multiplied with the density of the fluid; φ [m 3 /s] is the flux of fluid across a face.The cumulative overtopping volume V [m 3 ] can be obtained based on the overtopping volume flux by using the equation:

Fig. 6
Fig.6presents the comparison of incident free surface elevation time series given by the OpenFOAM® model and the experiments for cases SS1 and SS7.The incident water surface elevation time series were obtained based on the time series of free surface elevation at WG1, WG2 and WG3 using the method proposed byMansard and Funke (1980).The agreement between numerical and physical model results for both cases is generally good with a root mean square error (RMSE) of 0.0184 m for case SS1 and a RMSE of 0.015 m for case SS7.Similar good agreement is found for all other simulated cases, which demonstrates that the numerical model can well reproduce the wave propagation.An accurate reproduction of the wave energy spectra of incident waves in front of the structure is essential for an accurate prediction of the overtopping discharge at dikes.Fig.7shows the comparisons between modelled and experimental wave energy spectra of incident waves for case SS1 and case SS7.The experimental wave spectra are obtained based on the measured time series of incident wave surface elevation of 0-500 s corresponding to the simulation time in the numerical model.The shape of the wave energy spectra of measured incident waves from physical model is generally well predicted using the numerical model.Based on the wave energy spectra, the significant wave height H m0 and spectral wave period T m− 1,0 are calculated.Fig.8 and the ratio of Tm− 1,0− OpenFOAM Tm− 1,0− experiment for all 28 simulated cases.Scattered data points closely around line of y = 1 in Fig. 8a demonstrate that the significant wave height H m0 can be well reproduced by the OpenFOAM® model.The comparison of the spectral wave period in Fig.

Fig. 6 .
Fig. 6.Comparisons of time series of free surface elevation (η) given by the numerical model and experiments for the cases SS1 and SS7.
Chen et al. (2020b) developed an empirical equation (C.7) as shown in Appendix C to account for the berm influence on the average overtopping discharge based on the berm equation (C.2) provided by TAW(

Fig. 10 .
Fig. 10.Comparisons of numerical and experimental dimensionless average overtopping discharges q * at smooth straight slope (SS), smooth slope with berm (SB), straight slope covered by protruding blocks (PbS), bermed slope covered by protruding blocks (PbB) and bermed slope with only upper slope covered by protrduding blocks (PbB_up).
waves at the toe of the structure.R u2% can be calculated using equations (C.5) and (C.6).The empirical equation matches well with the experimental data points marked with open squares.It can be seen from the model results that the berm has the most reductive influence on the average overtopping discharge when the berm is located near the still water level while the reductive influence decreases as the distance between the berm and SWL increases.This is in accordance with the expressions for the influence of the berm as described inChen et al. (2020b).For the berm width of 0.2 m (corresponding to B/ H m0 = 1.6) used in the experiments, the berm equation works reasonably.Nevertheless, there are still some discrepancies between the numerical model results and the empirical equations, especially for the comparison between the computed berm factors and the berm equation for larger berm width (B/H m0 = 3.2 with B = 0.4 m), which is outside of the experimental ranges.Figs.12 and 13 demonstrate that the numerical model is capable of simulating the behaviour of berm influence on the average overtopping discharge at dikes to a reasonable extent.Additionally, differences between the modelled results and the empirical equations show a potential modification ofChen et al. (2020b) berm equation for a wider range of berm width and berm level.The berm equation (C.7) indicates a linear relationship between γ b and rB(1− r dh ) ̅̅̅̅̅̅̅̅̅ sm− 1,0 √ .Therefore, we plotted rB(1− r dh ) ̅̅̅̅̅̅̅̅̅ sm− 1,0 √ on the horizontal axis, where r B is a parameter taking the berm width into account and r dh is a parameter that represents the effect of the berm level on the berm influence factor, and the berm influence factor based on model results on the vertical axis (Fig. 14).This shows a linear relationship between γ b and rB(1− r dh ) ̅̅̅̅̅̅̅̅̅ sm− 1,0 √ but with a different value of the empirical coefficient b 0 .By applying the least square method, the empirical coefficient b 0 was recalibrated as 0.14 based on the numerical data for wider berms and larger berm levels that are outside of the ranges of experimental data.The best-fit formula based on the model results is defined as follows.

Fig. 12 .Fig. 13 .
Fig. 12.Effect of berm width on the berm influence factor given by OpenFOAM® model in comparison with empirical equations.

Fig. 15 .
Fig. 15.Mesh for a) one protrusion on the upper slope, b) two protrusions on the upper slope and c) three protrusions on the upper slope about half of the upper slope).

Fig. 16 .
Fig. 16.Comparison between experimental and numerical dimensionless average overtopping discharges over slopes with the upper slope covered by protruding blocks and slopes entirely covered by protruding blocks.
W.Chen et al.

Appendix C .
Empirical equations for estimating the average overtopping discharge at dikes Chen et al. (2020b) recalibrated the TAW (2002) overtopping equation for breaking wave condition by recalibrating one empirical coefficient in the TAW (2002) overtopping equation.The recalibrated overtopping equation (C.1) was then used as a reference formula to investigate the influence of a berm and roughness on wave overtopping discharges at dikes.m− 1,0 γ b γ f γ β γ v ) (C.1) in which q [m 3 /s/m] is the average overtopping discharge; α is the angle of waterside slope; ξ m− 1,0 (ξ m− 1,0 = tan α/ ̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ 2πH m0 /(gT 2 m− 1,0 ) √ ) is the breaker parameter; R c [m] is the crest level relative to the still water level; γ b [− ] is the influence factor for berms; γ f [− ] is the influence factor for roughness; γ β [− ] is the influence factor for oblique waves; γ v [− ] is the influence factor for vertical walls.TAW (2002) and EurOtop (2018) provide equations for estimating the berm influence factor γ b .γ b = 1 − r B (1 − r dh ) if 0.6 ≤ γ b ≤ 1.0 (C.2) where L berm[m]  is the characteristic berm length; d h[m]  is the berm level relative to the SWL; R u2%[m]  is the wave run-up height that is exceeded by 2% of the number of incoming waves at the toe of the structure which can be calculated by using the following equations(TAW, 2002):R u2% H m0 = 1.65γ b γ f γ β ξ m− 1,0(C.5) .(2020b) improved the TAW (2002) berm equation (C.2) by taking the wave steepness into account based on the analysis of experimental data.b 0 is an empirical coefficient and b 0 = 0.21 for an impermeable berm; r B and r dh are calculated using Eqs.(C.3) & (C.4).CRediT author statement Weiqiu Chen: Writing-original draft preparation, Software, Formal analysis, Conceptualization ; Jord J. Warmink: Conceptualization, Methodology, Writing-Review & Editing, Supervision; Marcel R.A. van Gent: Conceptualization, Resources, Writing-Review & Editing, Supervision; Suzanne J.M.H. Hulscher: Conceptualization, Writing-Review & Editing, Supervision, Project administration.

Table 1
Mesh resolution for different regions.

Table 2
Wave conditions simulated in the OpenFOAM® model.
W.Chen et al.

Table 4
Average contributions of protruding blocks for four wave conditions to the total roughness influence.