Hydrodynamic Modelling of Wave Overtopping over a Block-Covered Flood Defence

Wave overtopping can cause erosion on the landward slope due to high flow velocities and turbulence that cause high stresses on the cover. Innovative block revetments such as Grassblocks protect the subsoil of the dike against erosion. The blocks are permeable, which reduces the flow velocity and the pressures along the landward slope. The performance of these blocks is assessed in physical tests, which provides insights into the stability of the blocks. However, such experiments are expensive and accurate measurements are difficult due to highly turbulent conditions. Therefore, the goal of this study is to determine the hydrodynamic conditions at the dike cover caused by the wave run-up on the seaward slope and by the overtopping flow over the crest and landward slope. The geometry and wave conditions from the physical test at the Deltares Delta flume are implemented in an OpenFOAM® numerical model. Using the porousWaveFoam solver, a porous layer on the crest and landward slope is implemented, where the flow resistance of this porous layer largely depends on the resistance coefficients α [-] and β [-]. The numerical model is calibrated based on resistance coefficients as introduced earlier in the literature, which showed that the resistance coefficients of α = 500 and β = 2.0 performed best for the peak flow velocities and the peak pressures. The numerical model is evaluated by using these resistance coefficients in other time series of the physical tests. The evaluated model is then used to determine the hydrodynamic conditions on the landward slope, which showed that the pressure was the most influential hydrodynamic condition at the time of failure. Finally, the model showed that a porosity of n = 0.6 and the porous layer thickness η = 36 mm reduced the peak pressure the most.


Introduction
Wave overtopping is one of the main failure mechanisms of dikes. The wave overtopping process consists of three different phases: wave impact, wave run-up, and wave overtopping. The waves can overtop over the dikes during storm events. When the water reaches the seaward slope, crest, or landward slope of the dike, it can cause erosion due to high flow velocities and turbulence [1][2][3][4]. Dike cover erosion due to wave loads by overtopping waves could eventually lead to dike breaching, which can result in severe losses of properties and lives. Experiments [5][6][7][8][9][10] and numerical models [11,12] have shown that erosion on dikes is mainly caused by wave loads, such as pressure, normal stress, shear stress and flow velocities. Therefore, accurate estimations of these wave loads are important for dike designs and safety assessments.
Innovative covers strengthen the earthen dike cover against wave run-up and overtopping flow. When the grass cover erodes during extreme storm events, the blocks will still protect the underlying subsoil of the dike against erosion by overtopping waves (Figure 1). Blocks such as Grassblocks ( Figure 2) have an open volume between the blocks to absorb the up-and down-rushing water and dissipate the wave energy. This open volume per square meter defines the porosity of the system. Physical overtopping tests were performed to assess the performance of these blocks [13]. These large-scale physical tests provide insights into the stability of the blocks, but such tests are expensive. Furthermore, it is difficult to measure the flow conditions, such as flow velocity and pressure due to the highly turbulent conditions and the limited number of measuring locations. Additionally, it is time consuming to change the dike configurations in the physical tests.  [13]. (b) Grassblocks installed on crest and landward slope in the Deltares Delta flume [13].
Numerical models can provide the flow characteristics on any location on the dike cover and, therefore, can provide insight into the flow conditions on the dike and the limits of the dike cover. Previous research, e.g., Jacobsen et al. [15], Van Bergeijk et al. [16], and Chen et al. [17], has shown the capability of hydrodynamic numerical modelling in simulating the wave overtopping flow. However, these models were applied to simulate the flow on the waterside of the flood defence [15,17,18] or on the landward side [11,16]. The models for the landward slope use empirical volume distributions at the dike crest as the input for modelling the wave load along the landward slope, which can lead to inaccurate estimations of the wave load. Furthermore, the connection to the design conditions such as the water level and wave height is missing in these landward slope models. A hydrodynamic numerical model that has been validated for accurately simulating the hydrodynamic conditions due to wave overtopping over the waterside slope, crest, and landward slope of dikes, including permeable block revetments, does not exist yet. Therefore, the main goal of this paper is to make a reliable model for predicting the wave overtopping load along the dike cross-section and to use this model to quantify the physical stresses during the wave overtopping process. Using this new numerical model, we will access the strongest configuration of porous layer thickness and porosity of the block revetment in order to protect the dike against wave overtopping erosion.
For that purpose, a 2DV numerical model was set up using OpenFOAM®, in which a porous layer was implemented to represent the porosity of the Grassblock revetment. The OpenFOAM® model was calibrated and evaluated by comparing the numerically modelled flow velocities and pressures with the experimental data of the large-scale physical tests of Van Steeg [13]. The evaluated numerical model was then applied to determine the optimal configuration of the block revetment for reducing the overtopping forces.
The paper is structured as follows: Section 2 describes the physical model set-up, whereas Section 3 shows the numerical model set-up, calibration, and evaluation. In Section 4, the results are presented followed by the discussion and the conclusions in Sections 5 and 6, respectively.

Physical Model Tests
Two flume tests were performed in the Deltares Delta flume by Van Steeg [13] to determine the strength of the Grassblock against the wave forces. The waves were generated using a wave paddle that can generate regular and irregular waves with a significant wave height H s up to 2.0 m and is equipped with an active wave absorption system that prevents wave reflection [19]. Figure 3 shows the test set-up with the location of the wave generator paddle, the dike, the wave height meters, and the flow velocity and pressure sensors along the dike crest and landward slope. The toe of the dike structure was located at a cross-flume distance of X = 164.19 m. The lower waterside slope of the dike consists of concrete, followed by a section of Grassblocks (X = 190.29 m, Z = 8.95 m) on the waterside slope covering the crest until the landward slope (X = 201.79 m, Z = 6.20 m). The Grassblocks that were used in this study are permeable and were installed on a layer of geotextile with a clay layer of 0.4 m underneath [13]. The free surface elevation of the waves was measured by three wave height meters (also referred to as wave gauges) located at 108.5 m, 114.5 m, and 117.5 m, respectively, from the wave paddle.  The flow velocity was measured with a frequency of 100 Hz by two paddle wheels (PW) which were installed on the landward slope of the dike (Figure 3b). The centres of the paddle wheels were located at 27 mm above the top of the Grassblocks. The pressure was measured on five locations on the dike by pressure sensors (PS) with a frequency of 100 Hz. The pressures at PS2-PS5 were measured at 8 mm above the top of the Grassblocks and at 30 cm from the wall of the flume.
The test program consisted of two different tests in which the wave load was increased by increasing the wave height and reducing the wave steepness, see Table 1. The irregular waves were generated based on the JONSWAP spectrum. The forces on the blocks were only caused by water flowing over the blocks, because the overtopping water did not separate at the landward crest line and therefore the blocks did not experience any impact forces (Figure 4a). After Test A, no damage was observed on the Grassblocks. After one critical overtopping wave at 300 s in Test B the Grassblocks collapsed after which this test was aborted (Figure 4b). The collapse started at the toe of the landward slope, directly followed by the collapse of the block formation from the landward toe towards the crest.

Method
The measurements of the physical wave overtopping tests by Van Steeg [13] are used to set-up, calibrate, and evaluate the numerical model for a range of wave load conditions. After the evaluation, the hydraulic loads on the blocks were determined using the numerical model to determine the block revetment specifications leading to the smallest hydraulic forces. This section gives a description of the model set-up (Section 3.1), model calibration (Section 3.2), model evaluation (Section 3.3), and hydraulic loads on the blocks (Section 3.4).

Numerical Model Set-Up
OpenFOAM® 1712 was used in this study to set up the numerical model. Ocean-Wave3D coupled with the waves2Foam toolbox developed by Jacobsen et al. [15] was used to generate the same waves in the numerical model as those in the physical model tests. Detailed information about the coupling approach can be found in Paulsen et al. [20]. The movement of the wave paddle at the Deltares Delta flume is controlled by a steering file that describes the movement of the wave paddle to generate the wave time series. The steering files were used to translate the movement of the wave paddle into the OpenFOAM® model using OceanWave3D [21]. Since the Grassblocks that were used in this study are permeable, the Volume-averaged Reynolds-averaged Navier-Stokes (VRANS) equations were solved using the porousWaveFoam solver [18], which is included in the waves2Foam toolbox using the resistance coefficients α and β (see Appendix A).
The dike geometry was implemented into the OpenFOAM® model using blockMesh ( Figure 5). The flume in the numerical model starts at 42.5 m away from the wave generator up to 206 m with a height of 12.1 m. The relaxation zone avoids internally reflecting waves at the wave generating boundary of the numerical model. The length of the relaxation zone was set as close to one wavelength [22]. According to Jacobsen et al. [15] and Pedersen et al. [23], maintaining an aspect ratio of ∆x/∆z = 1 throughout the computational domain improves the accuracy and numerical stability of the model. The mesh size is 30 × 30 cm, leading to 40 cells in the vertical direction. The mesh size on the crest and landward slope is 7 × 7 cm. Similar to the research by Chen et al. [17], a wall function was used in the boundary layer near the bottom surface of the flume. The mesh was refined near the bottom along the entire domain using 10 cells with a height of 1 cm to simulate the vertical velocity profile accurately.
Windt et al. [24] have shown that 10 layers of cells per wave height are needed to resolve the wave height accurately. The significant wave height varies from 1.48 m up to 1.69 m. With a minimum of 10 vertical layers per wave height, this leads to a mesh size of 15 × 15 cm. Therefore, the grid was locally refined near the water surface from Z = 4.9-8.9 m in vertical direction and from the start of the flume until the start of the slope (X = 42.5-164.19 m) in horizontal direction ( Figure 6). The maximum Courant number was set to 0.25 as Gruwez et al. [25] showed that C = 0.25 led to a good balance between accuracy and computational costs. Van Gent [26] proposed the resistance coefficients α = 1000 and β = 1.1 based on measurements of the permeability during experiments with porous flows through rubble-mound material. Losada et al. [27] showed that the resistance coefficients depend on parameters such as the Reynolds number, the shape of the stones, the grade of the porous material, and the permeability and the flow characteristics. However, the precise descriptions of the α and β coefficients are still not fully understood for oscillatory flows in wave propagation and breaking over slopes [27]. For these conditions, the values as proposed in literature may not be valid since the experimental conditions for obtaining those formulae were not considering these effects. Based on a comparison of experimental data and numerical modelling results of wave overtopping processes (COBRAS-UC), Losada et al. [27] proposed α = 200 and β = 0.8 as the best-fit parameters in his research. Furthermore, Jensen et al. [18] compared experimental data with numerical outcomes. Based on research by Burcharth and Andersen [28], Jensen et al. [18] proposed that the Reynolds number in porous media Re p can be determined using Equation (1).
Based on this Reynolds number, a distinction was made by Jensen et al. [18] for different flow regimes. The first is a non-linear flow regime, also described as the Forchheimer flow regime (10 < Re p < 150). The second is an unsteady laminar flow regime, also called the transitional flow regime (150 < Re p < 300). The last one is the fully turbulent flow regime, with Re p > 300. According to Jensen et al. [18], the resistance coefficients α = 500 and β = 2.0 performed best considering all flow regimes.
Three different sets of resistance coefficients were selected based on Van Gent [26], Losada et al. [27], Jensen et al. [18] to calibrate the numerical model, as displayed in Table 2.  [27] 200 0.8 Since wave breaking was observed in the physical experiments, the turbulence should be accounted for in the numerical model. The stabilised k − ω turbulence model developed by Larsen and Fuhrman [29] was used in this study, where k is the turbulent kinetic energy [m 2 /s 2 ] and ω the specific rate of dissipation of turbulent kinetic energy. The k − ω turbulence model has a stress limiter λ 1 and an effective potential flow threshold λ 2 . Larsen and Fuhrman [29] suggested a value of λ 2 = 0.05 and for λ 1 either 0.2 and 0.875 is mentioned. Chen et al. [17] showed that λ 1 = 0.2 performed better in a k − ω model for wave overtopping. The stabilised k − ω model with λ 1 = 0.2 and λ 2 = 0.05 was used in this study.
The generated waves were calibrated following the approach of Chen et al. [17] to ensure that the modelled waves accurately represent the waves from the physical experiment. Wave height meters were defined in the numerical model at the same location as used in the physical experiment. Using the wave separation method of Mansard and Funke [30], the spectral density of the free surface elevation was determined. Based on the m 0 and m −1 moment of the spectral density, the significant wave height H m0 = 4 √ m 0 and spectral wave period T m−1,0 = m −1 /m 0 were obtained. The waves were then calibrated by comparing H m0 and T m−1,0 of the numerical model simulated waves and the physical test waves.

Calibration of the Resistance Coefficients
The pressures and flow velocities were calculated in OpenFOAM® using probes at the same locations as the pressure sensor and paddle wheels in the physical test. The first half of Test B (t = 0-170 s) was used for calibration. The peak flow velocity and peak pressure of each overtopping wave were assessed, because the peak values of the flow velocity and pressure have the most influence on the stability of the dike cover. An overtopping wave was defined as the exceedance of a minimum pressure of 0.2 kPa and a minimum flow velocity of 0.5 m/s. The modelled peak pressures and velocities of each run with resistance coefficients from Van Gent [26], Jensen et al. [18], and Losada et al. [27] were compared with the peak pressures and velocities as measured in the physical tests. The goodness of the obtained fits for each set of resistance coefficients are scored using the Nash-Sutcliffe Efficiency (NSE) and the root-mean-square error (RMSE), which are often used in previous studies (e.g., Chen et al. [31], Van Bergeijk et al. [4], and Van Gent [32]). A perfect fit between the model results and the formulation results in a NSE of 1 and a NSE of 0 indicates that the formulation is as accurate as the mean of the model results [33]. The RMSE shows the average deviation of the modelled values compared to the observed values. The resistance coefficient set with the highest NSE and smallest RMSE is selected for the further simulations in this study.

Evaluation of the Resistance Coefficients
The evaluation of the numerical model was performed using the test data of Test A (t = 0-500 s) and the second half of Test B (t = 170-350 s). Although the block revetments collapsed after 300 s in Test B, the test was stopped after 350 s. Therefore, the evaluation also partly includes the hydrodynamic conditions on the collapsed revetment. Similar to the calibration phase, an overtopping wave is defined as the exceedance of a minimum pressure of 0.2 kPa and a minimum flow velocity of 0.5 m/s. The resistance coefficients set of Jensen et al. [18] using α = 500 and β = 2.0 has been evaluated by comparing the modelled peak flow velocities and peak pressures of each overtopping wave with the peak pressures and velocities of each overtopping wave as measured in the physical tests.
The modelled peak flow velocities and pressures are then assessed based on the NSE and RMSE values.

Hydraulic Loads on Blocks
The evaluated numerical model was used to determine the forces that occurred on the landward slope during Test B and caused the block revetments to collapse after 300 s in Test B. The location of failure during the test was at the toe of the landward slope of the dike. Around this moment the hydrodynamic conditions pressure P [kPa], flow thickness h [m] and flow velocity u [m/s] were determined at the toe of the landward slope using the numerical model.
Subsequently, the evaluated numerical model was applied to investigate the sensitivity of the hydrodynamic conditions for the porosity and thickness of the permeable layer. Multiple simulations were completed where the peak pressure P, peak flow velocity u, and peak flow thickness h are calculated at the toe of the landward slope around the moment of failure (t = 290-310 s). As displayed in Table 3, the porosity n was varied between 0.2 and 0.6 with increments of 0.1. Separately, the porous layer thickness η was varied between 0.016 m and 0.056 m with increments of 0.01 m. The relative change in peak values is calculated using the difference of the calculated peak values and the default peak values for porosity n = 0.4 and porous layer thickness η = 36 mm, respectively. The results of this sensitivity analysis can be used in practice to determine which block revetment properties reduce the hydraulic loads the most.

Validation of Wave Conditions
The calibration of the free surface elevation was based on the significant wave height and spectral wave period. The OpenFOAM® model shows great performance for modelling the free surface elevation with an inaccuracy of 0.001 m for both tests. The spectral wave period was accurately modelled, only with a slight overestimation of 0.2 s for Test B (Table 4). Figure 7 shows the comparison of the measured and modelled incoming propagation of the free surface elevation on the location of the first wave gauge (X = 108.5 m) between 200 and 250 s of Test B. This signal shows that wave heights were accurately simulated.   Figure 8 shows the measured time-varying flow velocity and the numerically modelled flow velocity using the three different sets of resistance coefficients between t = 150-170 s. It can be seen that the numerical model can capture the timing and the variation of the flow velocity reasonably well. Figure 9 shows the one-to-one comparison of the modelled peak flow velocities and peak pressures per overtopping wave from t = 0-170 s in Test B. The three different resistance coefficients sets perform similarly for the flow velocity as seen in Figures 8 and 9a and resulted in comparable NSE (0.332, 0.315, 0.299) and RMSE (0.874, 0.885, 0.895) in Table 5.

Calibration of the Resistance Coefficients
As can be seen in Figure 9b, the modelled flow velocities are more accurate than the modelled pressures. Table 5 shows that using the Jensen et al. [18] (NSE = 0.266) and Losada et al. [27] resistance coefficients (NSE = 0.225), the model performed reasonably well, whereas the Van Gent [26] resistance coefficients (NSE = −0.024) resulted in less accurate pressures. Overall, these small differences in performance show that the model already performs reasonably well without too much calibration. The resistance coefficient set as suggested by Jensen et al. [18] was selected for the further simulations in this study based on the comparison of the NSE and RMSE values.    Figure 10 shows the one-to-one comparison of the modelled peak flow velocities and peak pressures per overtopping wave during the evaluation from t = 170-350 s in the test using resistance coefficients of Jensen et al. [18] (α = 500, β = 2.0). The modelled peak flow velocities (Table 6) show less scatter for PW1 (NSE = 0.681) than for PW2 (NSE = 0.502), which shows that the flow velocities on the crest are modelled more accurately than lower on the slope. With a total NSE of 0.606 the peak flow velocities were predicted accurately.

Evaluation of the Resistance Coefficients
The difference between the modelled and measured pressure shows more scatter with an NSE of 0.125 (Table 7). Furthermore, there is quite some variation between the pressure sensors with NSE coefficients ranging from −1.192 up to 0.546. The pressure of PS5 has the highest NSE value, which is close to the location of interest (toe of landward slope). The pressures of PS2 and PS4 have below zero NSE values, which might be caused by wall effects (Section 5). In general, the numerical model simulating wave overtopping over a permeable block layer performed reasonably well, whereas the modelled peak flow velocities showed higher accuracy than the modelled peak pressures, which will be discussed in more detail in Section 5.  Table 8 shows the five largest overtopping waves in Test B. The flow velocity, flow thickness, and pressure were simulated at the location where damage occurred, which is at the toe of the landward slope. This shows that especially at 300 s the pressure is larger than in other peaks in Test B, which is caused by the relatively high flow velocity and flow thickness occurring at the same time. The high flow velocities result in high dynamic pressures (p d = 1 2 ρ|u| 2 ) and the large flow thickness in high static pressures (p s = ρgh). These pressures combined result in a relatively large total pressure (p 0 = p s + p d ), which exceeded the limit state of the Grassblock revetments. This is in accordance with the experimental observation that the failure of the Grassblocks occurred around t = 300 s in Test B.  Considering the overtopping wave after t = 300 s in Test B led to failure of the Grassblocks, this specific overtopping wave has been used for the sensitivity analysis. The results of the porosity based sensitivity analysis are displayed in Figure 11a using n = 0.2 up to n = 0.6, resulting in a relative change in porosity of −50% to +50% compared to the reference case of n = 0.4. The peak flow velocity at approximately t = 300 s shows to depend on the porosity, where the peak flow velocity increases as the porosity increases. Furthermore, the peak flow velocity occurs earlier for a higher porosity (300 s for n = 0.6) compared to a lower porosity (302 s for n = 0.2). This can be explained by the water flowing faster through the porous layer at a higher porosity, which leads to the water reaching the toe of the landward slope earlier. The flow thickness decreases when the porosity increases. This can be explained by more water flowing through the porous layer. The only outlier in this case is n = 0.3, which has a larger flow thickness than n = 0.2. The pressure time series shows a large dependency on the porosity.

Hydraulic Loads on Blocks
The results of the sensitivity analysis based on the thickness of the porous layer using 16,26,36,46, and 56 mm resulting in a relative change in porous layer thickness of −55.5% to +55.5% compared to the reference case of η = 36 mm are displayed in Figure 11b. It shows that the peak flow velocity decreases when the thickness of the porous layer increases. Increasing the porous layer thickness leads to more porous volume, which therefore gives more resistance to larger volumes of water, thus decreasing the flow velocity. The opposite effect occurs for the overtopping flow thickness, which increases as the porosity increases. Noticeable in the pressure time series is that the current porous layer thickness results in the smallest peak pressures.

Towards Improved Modelling of Porous Dike Revetments
The results showed that we were able to construct a reliable numerical model in order to simulate the hydrodynamic conditions along the whole cross-dike profile. For computational benefits, the numerical model was set up as an 2D-vertical (2DV) model. In the physical test, the side walls of the Deltares Delta flume create friction, which leads to 3D effects [34]. This means that the flow velocity in the middle of the flume might be larger than near the walls of the flume. The 2DV model represents the middle of the flume, the location with the smallest amount of wall friction. However, the measuring equipment was installed near the side walls of the flume, with a distance of 30 cm. The 2DV numerical model excludes the influence of side walls. This might explain a timing and magnitude difference between the model output and physical test data. In future physical testing it might be desirable to have measurements from the centre part of the flume to minimise the effect of wall friction in future physical tests.
In the numerical model, the composition of the revetment layer below the block revetments was slightly different than in the physical test. In the physical test, the blocks were installed on a geotextile filter layer located on the subsoil consisting of sand, which may cause infiltration of water into the dike core. In the OpenFOAM® model, the infiltration of the overtopping water into the dike was not implemented. The additional calibration of the porosity and resistance coefficients of the filter layer, clay layer, and sand layer would be too complex and computationally expensive for this research goal. Therefore, the porous block revetments in the numerical model were assumed to be located on the same type of surface as of the waterside slope, which is concrete. Although this concrete surface in OpenFOAM® is impermeable, the geotextile filter layer, as used in the physical test, has a low permeability with a hydraulic conductivity of 0.07 m/s. Thus, although there is a difference in permeability between the physical and the modelled dike cover, it is a feasible approximation within this research. For future research, it is recommended to investigate the implementation of this infiltration in the dike to assess the effect of the permeability of the dike. Furthermore, a sensitivity analysis of the resistance coefficients α and β could give insight in which resistance coefficient is most critical for the peak hydrodynamic conditions. Koutrouveli and Dimas [35] found that for low-crested breakwaters the calibration of the α coefficient was much more crucial.

Performance and Limitations of the Model
The simulation of higher frequency waves turned out to be most difficult for the numerical model. This is also visible when we zoom in on the free surface elevation as plotted in Figure 7. Here, the data of the physical test show very small changes in surface elevation within a wave. The numerical model does not simulate these small bumps in surface elevation and has a more general approximation. However, these minor changes in surface elevation are not expected to have a significant influence on the wave propagation and wave overtopping. The peaks and the period of the surface elevation are most important and these are approximated well.
Although most of the overtopping waves in Test A were predicted well, there were some inaccuracies in the timing of the modelled pressures. A probable cause of this is the way the pressures were measured during the test. The pressure sensors have shown to be sensitive to the flow thickness and become less accurate for smaller flow thicknesses. In Test A, there was less wave overtopping, both in frequency and overtopping volume, resulting in less accurate pressure measurements. This is especially proven in Test B, where there was more wave overtopping and the pressure sensors functioned well. The results showed that the model performed best for Test B, which has the largest overtopping volumes and therefore the largest flow thicknesses.
The modelled peak pressures showed less agreement with the measurements compared to the modelled flow velocities. A possible explanation is the difficulties in the measurements of pressure during the physical test. These inaccurate pressure measure-ments might have influenced the performance indicators (NSE and RMSE) of the calibration and evaluation phase of this research. Due to the lesser amount of available data of pressure peaks, the sample size for calibration and evaluation was reduced. However, the numerical model was able to accurately simulate the timing of each wave overtopping peak.
The sensitivity analysis for porosity values showed that the porosity of n = 0.6 led to the most reduction in peak pressure. This coincides well with a sensitivity analysis on the porosity by Ren et al. [36], which showed that the effective value of porosity for a submerged breakwater is in the range from n = 0.4 to 0.8. However, it should be noted that the range of porosity in the sensitivity analysis of this research was from 0.2 to 0.6, not considering porosity values larger than 0.6. A drawback of the current model is that the model is two-dimensional and, therefore, 3D aspects could not be simulated. Another drawback is the computation time, it takes about 22 h to compute one simulation for about 300 s on a computer with 2 processors (2.5 GHz). In order to improve the computational efficiency we recommend to use the newly developed boundary condition by Wellens and Borsboom [37] for wave generation and absorption since this boundary condition eliminates the need of relaxations and therefore reduces the amount of cells, or to parallelise the computation in case more processors are available. However, such computational times are still less than executing physical tests in the laboratory.

Application
The numerical model can be used to assess the different kinds of flow characteristics on each location on the upper side of the waterside slope, the crest and the landward slope over time to provide more insight into the dominant flow characteristics. Furthermore, the dike geometry of the model can be easily adapted compared to physical flume experiments, which can then be used to assess the effect of different waterside slopes, crest widths, or landward slopes of the dike on the hydrodynamic conditions. Additionally, the wave conditions can be adjusted to determine how they affect the flow over the waterside slope, crest, and landward slope. The porous layer can also be adjusted to represent other types of block revetments, although the flow over these blocks depends on the grading and shape of the permeable material [38]. This could, therefore, require a calibration of the resistance coefficients α and β. However, the resistance coefficients of Jensen et al. [18] have proven to be valid for different kinds of flow regimes and could be widely applicable.
Using the model, the hydrodynamic conditions on the dike cover can be determined. These hydrodynamic loads are important for two types of dike failure: erosion of the dike cover and the block stability. The output of the numerical model can be used to determine the hydraulic load resulting in dike cover erosion [39]. Existing erosion models require hydraulic variables such as the flow velocity [4,40], the shear stress [3], or the normal stress [10], which can be calculated using this new numerical model. Furthermore, the model can be used for the design of innovative covers, such as described by Capel [41] and Chen et al. [42] leading to more reduction in the wave overtopping volume (waterside slope) or more erosion resistant covers (landward slope). The designs can be tested in the model and optimal specifications can be determined, which can optimise the block stability and reduce the amount of physical testing and, therefore, is more time-and cost-effective. An example is the sensitivity analysis performed on the grass blocks in this study where the hydraulic load on the blocks was calculated for different values of the porosity and the layer thickness. The characteristics and locations of the block revetment can easily be adapted in the numerical model to find the best design of the block revetment resulting in the smallest wave load.

Conclusions
A VRANS-based hydrodynamic model was built using a coupling of waves2Foam and OpenFOAM® by implementing the geometry and wave conditions from large-scale physical tests conducted at the Delta flume of Deltares. The porous function of the Grassblocks on the crest and landward slope was implemented into the model using a porous layer, where the porous function of this layer largely depends on the resistance coefficients α and β. Based on research by Van Gent [26], Losada et al. [27], and Jensen et al. [18], three sets of resistance coefficients were selected to calibrate the numerical model using a comparison of the modelled and observed peak flow velocities and pressures, which showed that the resistance coefficients of Jensen et al. [18] performed the best. The numerical model was then evaluated by using these resistance coefficients in other time series of the physical tests, which showed good results for the peak flow velocities (NSE = 0.606) and sufficient results for the peak pressures (NSE = 0.154). The calibrated and evaluated model was then used to determine the hydrodynamic conditions that led to failure of the block revetments. Model results show that the pressure was the most influential hydrodynamic condition at the time of failure. The model was also used to determine the influence of the porosity or the porous layer thickness on the hydrodynamic conditions, which showed that a porosity of n = 0.6 and the porous layer thickness of 36 mm reduced the pressure the most.
The hydrodynamic model has proven to be able to simulate the hydrodynamic conditions along the whole cross-dike profile. Although the peak pressures and flow velocities have proven to be difficult to simulate, the evaluated model could especially approximate the peak flow velocity well, whereas the peak pressure have shown to be approximated reasonably well. The model also proved to be a valuable tool in the design phase of block revetments by assessing which hydrodynamic conditions occurred. Using the model we were able to determine the properties of the block that result in the most reduction in the hydrodynamic loads on the dike cover.  is the excess pressure, x r is a reference location (defined at still water level), and µ tot is the total dynamic viscosity [18].
The flow resistance F p is described by the Darcy-Forchheimer resistance equation (Equation (A5)). where a and b are resistance terms as displayed in Equation (A6), defined by Van Gent [26].
Here, ν is the kinematic viscosity and d 50 is the nominal diameter of the porous material. The first term dominates in linear flow regimes, whereas the second term dominates in turbulent regimes. α and β are called the resistance coefficients.
KC is the Keulegan-Carpenter number and describes the relative importance of the drag forces over inertia forces in an oscillatory fluid flow. For small KC numbers the inertia dominates, whereas for large KC number the turbulent drag forces dominate. According to Jacobsen et al. [43], the KC number for revetment blocks is difficult to estimate. Ideally, the KC number should be used with a temporal and spatial distribution with a temporal variation to account for the changing hydrodynamics during an irregular time series. Since this is not feasible in OpenFOAM®, the KC number is determined based on the incident wave field and shallow water wave theory, as shown in Equation (A7).