Numerical and Physical Modelling of Wave Overtopping on a Smooth Impermeable Dike with Promenade under Strong Incident Waves

Abstract

This paper presents a study of run-up/overtopping over a smooth impermeable dike with promenade using 2D and 3D mesh-based and mesh-free numerical models and results from 2D physical modelling for strong energetic incident waves. These waves induce plunging wave breaking and a complex water/air mixture turbulent flow before overtopped the dike, a challenging configuration for numerical models. The analysis is structured in two phases: (i) evaluates the results of 2D numerical and physical models for run-up and overtopping; (ii) compares qualitatively the results of 3D numerical models for overtopping over a dike with promenade between groins located in front of a slope beach. The results indicate that the main differences obtained in run-up and overtopping are due to differences in wave generation and active absorption systems used in physical and numerical models and in turbulent models used by the numerical models. These differences lead to changes on incident wave height and on wave breaking and, consequently, on reflection, run-up and overtopping over the structure. For 3D simulation, even if larger discrepancies were found on overtopping along the dike, mean wave overtopping discharge and water flow height at the crest of the groin head show a similar order of magnitude.

1. Introduction

Flooding risk in coastal urban fronts protected by very shallow foreshores dikes is mainly dominated by strong wave action and high-water level. These wave conditions are challenging for numerical models, frequently breaking as plunging and inducing strong water/air mixture turbulent flow.

Physical modelling (e.g., [1,2,3,4]) is a good tool to study these cases, but it is time-consuming and expensive, especially if a 3D model is needed. Many authors presented physical model studies on overtopping at a dike with promenade [5,6,7,8,9,10]. However, only few studies are available in literature on wave overtopping prediction for dikes with very shallow, long, and gentle foreshores [11], especially for waves frequently breaking as plunging and inducing strong water/air mixture flow.

Several numerical models (e.g., [12,13,14,15,16]) are improving lately but need hard calibration/validation. Particularly for sloping beaches, the validation of the wave breaking dissipation, a very complex, non-linear, and energetic phenomenon, is essential to correctly obtain the run-up and overtopping over the structure. For the case of a beach protected by a dike constructed between groins, the phenomena involved are essentially 3D and, consequently, more difficult to correctly model numerically. However, it lacks analysis of numerical results of run-up and overtopping for strong wave conditions breaking as plunging. For example, IH2VOF [17,18,19,20], solves the 2DV Reynolds-Averaged Navier–Stokes (RANS) equations and the free surface flow is defined by the VoF (Volume-of-Fluid) method. IH2VOF was extensively validated for both surf zone hydrodynamics and the stability and functionality of conventional or non-conventional coastal structures (e.g., [21]). FLUENT [22] solves the 2D and 3D RANS-VoF equations using a Volume Finite method on conform/non-conform and structured/unstructured meshes and can be used for modelling wave-structure interactions. FLUENT was applied and validated for 2D and 3D applications, such as submarine outfalls [23], onshore and offshore wave energy converters [24,25,26], and wave overtopping of coastal structures [27]. DualSPHysics [28] is a Lagrangian method based on the Smoothed Particle Hydrodynamics (SPH), where the fluid is described as a set of discrete particles [29,30]. This method is considered a versatile method, naturally dealing with unsteady and nonlinear flows, extreme deformations, and complex topological evolutions, such as a breaking free surface, with sufficient accuracy to provide meaningful solutions to engineering problems [31]. DualSPHysics has some limitations such as solve complex boundary conditions and 3D cases, due to the large time consumption [32]. However, in literature, there is no application of any of these three numerical models to study run-up and overtopping over a dike due to wave breaking, as plunging and induced strong water/air mixture turbulent flow.

This work is focused on modelling strong waves breaking as plunging and the consequent run-up/overtopping over a smooth impermeable dike with promenade between groins located in front of a slope beach, using physical and numerical models. An impermeable structure is used to glean insight into the accuracy of the models avoiding the influence in the results of possible inaccuracies coming from the model/method used to simulate porous medium in the numerical models. The still water level used in the study, which leads to a large freeboard (almost the height of the smooth dike), and the 1:10 slope before the smooth dike contribute to the occurrence of plunging breaking far from the toe of the dike, generating a complex turbulent water/air mixture flow above the ramp. Two mesh-based models, IH2VOF and FLUENT, and a mesh-free model, DualSPHysics, are used.

The methodology has the following main components: (i) 2D physical model tests for obtaining breaking characteristics and measuring run-up and overtopping on a smooth impermeable dike for scenarios of sea level rise; (ii) verification of performance and accuracy of the three 2D numerical models comparing results with the physical model data; and (iii) 3D modelling, using the RANS-VoF and mesh-free models, of run-up/overtopping in a 1:15 slope beach protected by a smooth impermeable dike located between two impermeable groins. This is the first step for numerical modelling of a dike between groins, an important and urgent tool to assist in the perdition of flooding risk at coastal cities protected by structures and in studying alternative solutions to mitigate these risks.

2. Physical Model Tests

The physical tests were carried out in the Atmosphere-Ocean Interaction Flume (CIAO) of the Andalusian Institute for Earth System Research (IISTA)-University of Granada (Figure 1a). This wave flume is dedicated to the combined study of marine and atmospheric processes. The facility is equipped with two opposed piston-type generating paddles, with an active system to absorb reflection, and a closed circuit for wind generation. The wave flume has a width of 1 m, length of 16 m, and a range of generating waves of periods from 1 to 5 s and wave heights up to 25 cm (see Figure 1 and Figure 2).

Tests were performed with a VTI controller. AwaSys software package was used to generate waves with the simultaneously active absorption of reflected waves [33]. Active absorption system uses real-time measurements of water level in front of the paddles to obtain a desired reflective condition of regular and irregular waves. In these experiments, regular waves were generated using paddle 2 (Figure 2).

Physical model tests were conducted in a wave flume for a schematic dike constructed on land with a beach and promenade corresponding to common values on the Atlantic coast of Spain and Portugal. Figure 2 shows a scheme of the physical model, consisting of an impermeable ramp made of wood, with three slope angles of 1:10, 1:12, and 1:2 for the sea wall. The crest and the toe of the schematic dike are 78 cm and 58 cm above the floor, respectively. Three water depths, h, have been considered, of 0.50, 0.58, and 0.62 m, respectively corresponding to three sea level rise (SLR) scenarios: current situation, SLR = 0.45, and SLR = 0.90 m. A scale factor of 34.5 was selected.

Free surface elevation was measured using four UltraLab ULS 80D acoustic wave gauges placed along the flume (Figure 2, wave gauges WG1, WG2, WG3, and WG4). Overtopping was identified and measured through two UltraLab ULS 80D acoustic wave gauges (WG5, WG6). WG5 was used to measure overtopping height in the sea wall, and WG6 was used to measure overtopping volume and mean wave overtopping discharge. By measuring the depth of water inside the overtopping tank after each test, using WG6, and considering the dimensions of the tank, it was possible to evaluate the overtopping volume and the flow rate of each test. These gauges have a maximum repetition rate of 75 Hz, a space resolution of 0.18 mm, a working range of 350 mm and a reproducibility of ±0.15%.

Wave run-up was obtained by analysing video images from inside the wave flume, and from outside (Figure 3). For this purpose, tape measures were placed both on the side of the flume and along the slope (Figure 1b and Figure 3), which allowed, in the subsequent viewing of the videos, the evaluation of the maximum run-up reached by the waves. The accuracy of the measurement is that corresponding to the discretisation of the tape measure, 0.01 m.

2.1. Test Programme

Regular waves were simulated and defined by its wave height, H, and period, T, considering typical storm wave characteristics of Atlantic Spanish and Portuguese coasts, considering a scale of 1:34.5.

Table 1. Wave characteristics of the physical tests.

Test ID	h (m)	T (s)	H (m)	Duration (s)
1	0.50	2.04	0.116	204
2	0.50	2.04	0.23	204
3	0.50	1.50	0.23	150
4	0.58	2.04	0.116	204
5	0.58	2.04	0.23	204
6	0.58	1.50	0.23	150
7	0.62	2.04	0.116	204
8	0.62	2.04	0.23	204
9	0.62	1.50	0.23	150

shows the target wave parameters in each configuration and used in this paper.

Following [34], the results from regular waves can be extended to irregular waves via the hypothesis of equivalence introduced by Saville [35], empirically proven by Bruun and Gunbak [36] for run-up on rough permeable slopes.

2.2. Experimental Results

For each Test ID, three replicates were performed to analyse the repeatability of the test in the wave flume. As expected, small deviations were observed between tests, more significant in those cases where the generation system works closer to its operating limits. Thus, calculating the percentage difference between maximum and minimum values measured for the mean wave height, following the formula proposed by [37], differences between 0.4% in the most favourable case and 7.6% in the most unfavourable case are observed, within similar values to those found in the aforementioned work.

Free surface was measured on three points close to the paddle (Figure 4), which makes possible the separation between incident and reflected wave trains, and to calculate the value of the reflection coefficient.

The proper behaviour of wave generation can be evaluated comparing target wave height, H_target, versus the incident wave height, H_inc, obtained by applying the Baquerizo Method ([38], based on [39]), to separate incident and reflected wave based on the measured time series of surface elevation on the three wave gauges (Figure 5).

Maximum differences of 0.044 m were found between target and estimated wave height and the maximum difference reaches 19.9%. Maximum difference is observed in the case of h = 0.62 m (test IDs 8 and 9) and may be due to the fact that for those cases, the generation system works almost at the operating limit, so that the generation behaviour deviates slightly from what is expected, which is not the case when working within the optimal generation limits.

Maximum wave run-up, Ru_max, was measured in all the tests. Finally, overtopping volume was measured and mean wave overtopping discharge can be calculated. In the experimental tests, only the tests 4, 6, 7, and 8 showed measurable values of overtopping volumes.

Table 2. Maximum run-up and mean wave overtopping discharge obtained in the tests.

Test ID	h (m)	T (s)	H (m)	Rumax (m)	Q (L/s/m)
1	0.5	2.04	0.116	0.21	-
2	0.5	2.04	0.23	0.23	-
3	0.5	1.5	0.23	0.22	-
4	0.58	2.04	0.116	-	0.02
5	0.58	2.04	0.23	-	0.35
6	0.58	1.5	0.23	-	0.05
7	0.62	2.04	0.116	-	0.35
8	0.62	2.04	0.23	-	2.84
9	0.62	1.5	0.23	-	0.98

presents the maximum run-up measured in the cases where no overtopping has occurred, corresponding with the cases of h = 0.50 m (tests 1, 2, and 3) together with mean wave overtopping discharge, Q, calculated from the measured data of overtopping volume. Larger mean wave overtopping discharge occurs for test ID 8 reaching 2.84 L/s/m.

3. Model Characteristics

In this paper three, different models have been used: DualSPHysics [28], FLUENT [22], and IH2VOF [19,21]. A summary of the model characteristics is presented hereafter.

DualSPHysics is a hardware accelerated Smoothed Particle Hydrodynamics (SPH) open-source code [28]. In SPH, the fluid is described as a set of discrete particles where any physical property can be computed as an interpolation of the values of the nearest neighbouring particles. The contribution of the neighbouring particles is weighted according to their distance from a target particle using a kernel function and a smoothing length [29,30]. The kernel function is expressed in a discrete form, where the approximation of any physical property is interpolated at a given location and a summation is performed over all the particles within the region of compact support of the kernel. In this paper, a Quintic kernel suggested by [40] was used.

In DualSPHysics, the fluid is treated as Weakly Compressible SPH (WCSPH) and the equation of state is used to determine fluid pressure based on particle density [41,42]. The compressibility is adjusted so that speed of sound can be artificially lowered; this means that the size of time step taken at any one moment (which is determined according to a Courant–Friedrich–Lewy (CFL) condition, based on the currently calculated speed of sound for all particles) can be maintained at a reasonable value. Such adjustment, however, restricts the speed of sound to be at least ten times faster than the maximum fluid velocity, keeping density variations to within less than 1%, and therefore, not introducing major deviations from an incompressible approach. In this paper, the governing equations are integrated in time using Symplectic scheme. This scheme is time reversible in the absence of friction or viscous effects [43]. It can also preserve geometric features, such as the energy time-reversal symmetry present in the equations of motion, leading to improved resolution of long-term solution behaviour. The scheme used here is an explicit second-order Symplectic scheme with an accuracy in time of second order and involves a predictor and corrector stage. With explicit time integration schemes, the time step is dependent on the CFL condition, the force terms, and the viscous diffusion term, so that a variable time step is calculated according to [44].

The Dynamic Boundary Condition (DBC) is the default method provided by DualSPHysics [45]. The boundaries are described by a set of particles that satisfy the same equations as the fluid particles (mass and momentum conservation equations); however, they do not move according to the forces exerted upon them. An interesting advantage of these particles is their computational simplicity, since they can be calculated inside the same loops as fluid particles. However, they remain either fixed in position (e.g., fixed boundaries) or move according to an imposed/assigned motion function (i.e., moving bodies such as wave-makers) [46]. Validations with dam-break flows and sloshing tanks have been published with good results and comparing these boundary conditions with other approaches [46]. In addition, DBC have also been shown to be suitable to reproduce complex geometries [31] and moving bodies [32]. A piston-type wave-maker was used to generate waves following Stokes II. An active wave absorption system (AWAS) is used to absorb the reflected waves at the piston to mimic the behaviour of an open sea where reflected waves propagate outside the computational domain.

Numerical model FLUENT [22] is employed to solve RANS-VoF (Reynolds-Averaged Navier–Stokes Volume-of-Fluid) equations using a Finite Volume technique, where variables are defined in the centre of each control volume (CV).

RANS-VoF equations, based on the decomposition of the instantaneous velocity and pressure fields of the Navier–Stokes equations into mean and fluctuating components, and the subsequent time-averaging of the set of equations, are used. This process introduces Reynolds stress terms associated with the turbulence. The k-ω SST turbulence model is used for relating Reynolds stresses to mean flow variables and close the equations, since it seems more adequate for modelling wave-breaking [27].

Free surface flow is defined by the VoF method [47], which is based on the transport equation of the volume fraction, a scalar that takes values 0 in the air, 1 in the water, and 0.5 at the free surface.

In the FLUENT numerical model, the solver scheme SIMPLEC (with standard under-relaxation factors) is used and the scheme PRESTO! is recommended for discretizing pressure. The momentum is discretised by the third-order scheme MUSCL, and the turbulence kinetic energy and specific dissipation rate are discretised by the second order upwind scheme. The geo-reconstruct scheme, well adapted for modelling the complex shape of free surface flows, such as wave breaking and overtopping, is used for the VoF equation, compatible with the first-order time integration scheme and variable time step that depends on the CFL condition.

The RANS-VoF equations are solved in a structured/unstructured and conform/non-conform mesh. However, good accuracy of free surface is obtained using a structured mesh in the zone of free surface deformation and a regular mesh is recommended in the zone of run-up, wave breaking, and overtopping [27].

Wave generation is performed using a static numerical wave-maker (composed of several independent numerical paddles for 3D wave tank applications), imposing the velocity component profiles to the wave-maker boundary and the corresponding free surface position. An active absorption technique is imposed at the wave generation to eliminate reflected waves [24,48].

The non-slip condition is imposed on the structure and the bottom of the wave flume. An exit boundary is applied at the top of the breakwater. Free surface level at rest, null velocity components, hydrostatic pressure on the water, and atmospheric pressure on the air are the initial conditions.

IH2VOF solves the 2DV RANS equations based on the decomposition of the instantaneous velocity and pressure fields into mean and turbulent components at the clear-fluid region (outside the porous media) and inside the porous media by the resolution of the Volume-Averaged Reynolds-Averaged Navier–Stokes (VARANS) equations based on the decomposition of the velocity and pressure fields into mean and turbulent components using a k-ε turbulent model on a 2D vertical domain [17,21]. As in FLUENT, free surface flow is defined by the VoF method.

The RANS-VoF equations are solved in a rectangular structured mesh. Pressure, p, and velocity, u, variables are defined in a staggered scheme. The model uses a cutting cell method first introduced by [49] in order to consider obstacles contained within the numerical domain. The basic idea behind this technique is that the obstacle can be modelled as a special case of the flow with an infinite density [17].

Wave generation is performed here using a flat numerical wave-maker, imposing the velocity component profiles to the wave maker boundary and the corresponding free surface position. The numerical wave-maker is equipped with an active wave absorption, using the methodology proposed by [50]. IH2VOF is one of the most advanced RANS-VoF models thanks to its capabilities, robustness, and extensive validation for both surf zone hydrodynamics and the stability and functionality of conventional or non-conventional coastal structures.

Summing up, all models resolve the Navier–Stokes equations. However, DualSPHysics is a Lagrangian model that solves Navier–Stokes equations and FLUENT and IH2VOF are Eulerian models that solve RANS equations. IH2VOF is a 2D model where FLUENT and DualSPHysics have 2D and 3D versions. DualSPHysics has a dynamic paddle with active absorption, FLUENT has a static paddle (multiple static paddles for 3D wave tank) and IH2VOF can be used with a dynamic or static paddle, both models with an active absorption implemented. Meshes are also different in FLUENT (structured/unstructured and conform/non-conform mesh) and IH2VOF (rectangular structured mesh).

4. Sensitivity Analysis

Two cases are chosen to perform a sensitive analysis of each model: one with no overtopping, where the influence of the discretisation on the run-up is analysed (test 1), and another where the influence of the discretisation on the overtopping is analysed (test 5).

To select the discretisation of the models that allows us to obtain the better solution, a sensitivity analysis of the influence of the mesh discretisation on the results was carried out for FLUENT and IH2VOF models, and of the influence of number of particles on the results were acquired for DualSPHysics. For the three models, the results presented correspond to an analysis of 25 waves, i.e., between 29 and 80 s of simulation, in order to analyse exactly the same waves.

4.1. FLUENT

For FLUENT, a non-conform mesh is used since it allows a better optimisation of mesh refinement in the proximity of the breakwater and inside the wave breaking zone (wave breaking occurs around 9.0 m). Both regular incident waves analysed in this study have high nonlinear characteristics; therefore, the Fourier wave theory, which provides a good accurate for a wide range of cases, is used [51]. Seven meshes are used, from 40316 to 98533 control volumes (CV). For the first and second configuration, meshes 4 and 7 are used for mesh resolution sensitivity analysis, respectively.

Table 3. Mesh characteristics zone by zone with total number of control volumes and number of control volumes in x and y discretisation used in FLUENT.

Mesh	Number of CVs	Propagation Nx × Ny = Nº CVs dx × dy	Wave Breaking Nx × Ny = Nº CVs dx × dy	Beach Nx × Ny = Nº CVs dx × dy	Breakwater Nx × Ny = Nº CVs dx × dy
1	40316	271 × 76 = 20,596 50 × 6 mm	280 × 40 = 12,320 7 × 7 mm	140 × 40 = 5600 7 × 7 mm	72 × 25 = 1800 7 × 7 mm
2	41684	271 × 76 = 20,596 50 × 6 mm	280 × 40 = 12,320 7 × 7 mm	140 × 40 = 5600 7 × 7 mm	72 × 44 = 3168 7 × 3.5 mm
3	44852	271 × 76 = 20,596 50 × 6 mm	280 × 40 = 12,320 7 × 7 mm	140 × 40 = 5600 7 × 7 mm	144 × 44 = 6336 3..5 × 3.5 mm
4	45604	271 × 76 = 20,596 50 × 6 mm	280 × 40 = 12,320 7 × 7 mm	140 × 68 = 9520 7 × 3.5 mm	72 × 44 = 3168 7 × 3.5 mm
5	90757	309 × 76 = 23,484 50 × 6 mm	571 × 80 = 45,680 3.5 × 3.5 mm	275 × 67 = 18,425 3.5 × 3.5 mm	72 × 44 = 3168 7 × 3.5 mm
6	93925	309 × 76 = 23,484 50 × 6 mm	571 × 80 = 45,680 3.5 × 3.5 mm	275 × 67 = 18,425 3.5 × 3.5 mm	144 × 44 = 6336 3.5 × 3.5 mm
7	98533	309 × 76 = 23,484 50 × 6 mm	571 × 80 = 45,680 3.5 × 3.5 mm	275 × 67 = 18,425 3.5 × 3.5 mm	144 × 76 = 10,944 3.5 × 1.7 mm

shows the characteristics of the seven meshes, indicating the number of volume controls and the discretisation in horizontal (Nx) and vertical (Ny) directions zone by zone: wave propagation (0 < x < 9.1 m), wave breaking (9.1 m < x < 11.1 m), beach (11.1 m < x < 12.06 m) and breakwater zone (slope and crest) (12.06 m < x < 12.5 m). For wave propagation zone, mesh resolution in horizontal and vertical direction is dx = 0.05 m and dy = 0.006 m, respectively, which corresponds to a minimum of 94 control volume per wavelength and 20 control volumes per wave height and allows good accuracy of wave propagation [25,48,52]. The mesh resolution in this zone is the same for all meshes. In the wave breaking zone, two mesh resolutions are used with dx = dy = 0.007 m for meshes 1 to 4 and dx = dy = 0.0035 m for meshes 5 to 7, with special care to conserve a relatively regular mesh more adapted for wave breaking modelling. In beach and breakwater zones, combinations of 0.007 and 0.0035 m control volume dimensions are used as presented in Table 3. Mesh 7 shows the finer mesh with dy = 0.0017 m in breakwater zone.

Run-up case, test 1, is analysed for mesh 1, 4, 6, and 7.

Table 4. Run-up case: mean wave height at gauges located at x = 3.0, 8.0, 11.1, and 12.06 m and mean and maximum run-up obtained in FLUENT.

Mesh	Wave Height (m)				Ru (m)
	x = 3.0 m	x = 8.0 m	x = 11.1 m	x = 12.06 m	Mean	Max
1	0.1245	0.1370	0.0367	0.0082	0.098	0.116
4	0.1256	0.1366	0.0324	0.0090	0.100	0.124
6	0.1251	0.1438	0.0436	0.0151	0.104	0.139
7	0.1274	0.1411	0.0457	0.0185	0.124	0.171

resumes the mean wave height at four gauges located at x = 3.0 and x = 8.0 m, before the wave breaking, x = 11.1 m, after the wave breaking, and x = 12.06 m at the toe of the dike and the mean and maximum run-up (see Figure 6). Two gauges are also used at x = 2.8 and 3.3 m for calculating the reflection coefficient.

Wave height before wave breaking, at gauge x = 3.0 and 8.0 m, is about the same for the four meshes and the maximum difference is around 3%. After the wave breaking, at gauge x = 11.0 m, differences are 20% and 29% for mesh 1 and 4, respectively, and wave height for mesh 6 and 7 is similar with only 4.6%. Difference is more important between the coarse meshes 1 and 4 and the finer one 7 at gauge x = 12.06 m, around 55%. The complex flow above the beach zone is sensitive to the mesh refinement (meshes 6 and 7 have a more refine mesh in the x direction). Nevertheless, difference between mesh 6 and 7 is 19%—this is probably due to the mesh refinement on the dike in the y direction which allows modelling the small accumulation of water at the toe of the structure. Run-up reaches and oversteps the toe of the dike for the four meshes. Maximum differences occur between mesh 1 and 7, the coarser and finer mesh, with 21% and 44.3% for mean and maximum Ru, respectively. The sensitivity of run-up with mesh refinement in the beach and dike zones is clear and indicates than mesh refinement in the x direction is important for water flow even if refinement in the y direction is also important for capturing free surface position. Differences between mesh 6 and 7 for mean and maximum Ru are 15.7% and 25.3%, respectively, due to the meh resolution in vertical direction on the dike (dy = 3.5 and 1.7 mm for mesh 6 and 7, respectively). The mesh 7 seems to be the most appropriated for the simulation of run-up case.

The overtopping case, test 5, is analysed for mesh 1 to 7.

Table 5. Overtopping case: mean wave height at gauges located at x = 3.0, 10.0, 12.06, and 12.42 m and the mean wave overtopping discharge obtained in FLUENT.

Mesh	Mean Wave Height (m)				Q (10−3 m3/s/m)
	x = 3.0 m	x = 10.0 m	x = 12.06 m	x = 12.42 m
1	0.2554	0.1804	0.0997	0.0197	2.70
2	0.2552	0.1859	0.0919	0.0176	2.43
3	0.2570	0.1861	0.0974	0.0187	2.26
4	0.2538	0.1802	0.0982	0.0192	2.54
5	0.2564	0.1938	0.0988	0.0203	2.01
6	0.2486	0.1916	0.1040	0.0215	2.71
7	0.2515	0.1950	0.1029	0.0207	2.26

presents the mean wave height at four gauges located at x = 3.0 and 10.0 m, before the wave breaking, x = 12.06 m, at the toe of the dike, and 12.42 m, at the top of the dike, and the mean wave overtopping discharge. Wave height at gauge x = 3.0 m is almost equal for all meshes with only a difference smaller than 2% compared to results obtained for the finer mesh 7. For gauge x = 10.0 m, differences of wave height are slightly larger with 5% to 7% for meshes 1 to 4 and smaller than 1% for meshes 5 and 6. Wave height at the toe of the dike, x = 12.06 m, is almost equal for all meshes, with differences of 1% to 5%, even if flow motion after the plunging wave breaking shows complex structures (secondary wave breaking, convected vortex structures, reverse wave breaking, splash-up), except for mesh 2, which reaches 10.6% of difference. It seems indicating that mesh refinement in the wave breaking zone used for mesh 1 to 4 is enough to model correctly the wave breaking and flow motion above the beach zone. Water height at the crest of the dike shows a difference between 5% to 15% for meshes 1 to 4 and only 2% and 4% for meshes 5 and 6, respectively. Wave overtopping discharge was calculated for each wave based on the water height at the crest of the dike and the horizontal component of water velocity. Mean wave overtopping discharge was significantly over-estimated for meshes 1 and 6, with a difference around 19.5% compared to finer mesh 7, and under-estimated for mesh 5 with a difference of 11.3%. Results obtained for mesh 3 and 7 were equal even if mesh 7 was significantly more refined in breaking wave, beach, and breakwater zones than mesh 3. The same observation holds true for results obtained from mesh 1 and 6, with mesh resolution two times smaller for mesh 6 than mesh 1 in breaking wave, beach, and breakwater zones. Convergence of mean wave overtopping discharge with mesh refinement was not observed and can be explained by: (i) overtopping, and subsequently the water flow on dike, was very sensitive to mesh refinement in both directions and depends on water height and water velocity at the dike crest; and (ii) wave breaking presents a chaotic behaviour with a large variability of water height and flow field before and over the dike. Mesh 7, since it was the finer mesh, seems to be the best mesh for the simulation of overtopping case.

Based on this sensitivity analysis, all simulations in Section 5 are conducted using mesh 7.

4.2. IH2VOF

For IH2VOF, four different rectangular meshes were used for mesh resolution sensitivity analysis, with a minimum discretisation in two directions varying from 6 to 3 mm. Both regular incident waves were generated using Stokes III. The domain was divided into three horizontal zones and two vertical zones. In the horizontal direction, regular meshes were used for wave propagation (0 < x < 6 m), and beach (10 m < x < 13 m) zones, where an irregular mesh was used in the wave breaking (6 m < x < 10 m) zone. In the vertical direction, an irregular mesh was used in the lower part (0 < x < 0.43 m) and a regular mesh was used in the upper part (0.432 m < x < xfinal). The xfinal changes from grid to grid to guarantee that the water depth is coincident with the mesh as well as the top of the beach and is equal to 1.003 for mesh 2 and 4, 1.004 for mesh 3, and 1.006 for mesh 1. The grid characteristics of the different meshes are presented in

Table 6. IH2VOF mesh characteristics.

Mesh	Total Number of Cells	Propagation Nx (dx)	Wave Breaking Nx (dx)	Beach Nx (dx)	Lower Zone Ny (dy)	Upper Zone Ny (dy)
1	1601 × 162	600 (dx = 10 mm)	501 (10 mm < dx < 6 mm)	500 (dx = 6 mm)	72 (dy = 6 mm)	90 (6 mm < dy < 7 mm)
2	2041 × 214	800 (dx = 7 mm)	641 (7 mm < dx < 5 mm)	600 (dx = 5 mm)	96 (dy = 4 mm)	118 (4 mm < dy < 5.2 mm)
3	2551 × 244	1000 (dx = 6 mm)	801 (6 mm < dx < 4 mm)	750 (dx = 4 mm)	144 (dy = 4 mm)	100 (4 mm < dy < 5 mm)
4	3201 × 321	1200 (dx = 5 mm)	1001 (5 mm < dx < 3 mm)	1000 (dx = 3 mm)	168 (dy = 3 mm)	153 (3 mm < dy < 4 mm)

. The total number of cells together with the number of cells and their dimension in horizontal (Nx, dx) and vertical (Ny, dy) directions are presented for each zone of the domain.

Results are obtained in 11 numerical gauges, with at four in the first zone, three in the second zone, and four in the third zone.

For the run-up case, test 1, the time series of free surface elevation obtained at the numerical gauges are analysed for mesh 1 to 4.

Table 7. Run-up case: mean wave height at gauges located at x = 3.0 8.0, 11.1 and 12.06 m and the mean and maximum run-up obtained in IH2VOF.

Mesh	Wave Height (m)				Ru (m)
	x = 3.0 m	x = 8.0 m	x = 11.1 m	x = 12.06 m	Mean	Max
1	0.1190	0.1202	0.0359	0.0054	0.070	0.114
2	0.1173	0.1250	0.0493	0.0127	0.090	0.165
3	0.1191	0.1172	0.0521	0.0175	0.090	0.134
4	0.1174	0.1281	0.0534	0.0135	0.090	0.152

resumes the mean wave height at four gauges located at x = 3.0 and x = 8.0 m, before the wave breaking, x = 11.1 m, after the wave breaking, and x = 12.06 m at the toe of the dike and the mean run-up.

Wave height before wave breaking, at gauge x = 3.0 and 8.0 m, is about the same for the four meshes and the maximum difference is 6% for mesh 1 and reduced to 2% for meshes 2 and 3, when compared with mesh 4. As for FLUENT, the difference is more important between the meshes 1 and 4. After the wave breaking, at gauge x = 11.0 m, differences are clearly with the mesh dimension and are 33% between meshes 1 and 4, 8% between meshes 2 to 4, and reduced to 2% between meshes 3 and 4. At gauge x = 12.06 m, the wave thickness becomes very small and the reduction in differences with the refinement of the mesh are not as expected, since differences to the mesh 4 are reduced from mesh 1 (60%) to 2 (5%) but increase from mesh 2 to 3 (0.004 m, corresponding to 30%). However, this difference is, in fact, smaller than the finer mesh discretisation for mesh 3, that is 5 mm. Moreover, the wave height in this zone is very much influenced by small changes in the wave breaking dissipation, leading to differences in run-up at the slope of the structure. In fact, looking to the maximum run-up in these 25 waves, the differences follow the same pattern as differences in wave height at x = 12.06 m: the maximum difference between each mesh and mesh 4 varies from 25% for mesh 1, 9% for mesh 2, and 12% for mesh 3. However, mean run-up is only different for mesh 1, with a value 22% lower than the meshes 2 to 4, that present the same value of mean run-up. Taking all the results for this case, mesh 2 and 3 are appropriated for the simulation of run-up case.

For the case with overtopping, test 5, results obtained with meshes 1 to 4, the same 25 wave periods were analysed.

Table 8. Overtopping case: mean wave height at gauges located at x = 3.0, 10.0, 12.06, and 12.42 m and mean overtopping obtained in IH2VOF.

Mesh	Mean Wave Height (m)				Q (10−3 m3/s/m)
	x = 3.0 m	x = 10.0 m	x = 12.06 m	x = 12.42 m
1	0.1758	0.2651	0.0633	0.0117	0.35
2	0.1864	0.2863	0.0694	0.0127	0.94
3	0.1844	0.2801	0.0668	0.0161	1.80
4	0.1893	0.2777	0.0793	0.0174	1.27

presents the mean wave height at four gauges located at x = 3.0 and 10.0 m, before the wave breaking, x = 12.06 m, at the toe of the dike, and 12.42 m, at the top of the dike, and the mean wave overtopping discharge.

Mean wave height at gauge x = 3.0 m and 10 m is almost equal for all meshes with only a difference smaller than 7% for mesh 1 and smaller than 3% for meshes 2 and 3, compared with results obtained for the finer mesh, mesh 4. Wave height at the toe of the dike, x = 12.06 m, presents higher differences of 20% for mesh 1 and less than 16% for meshes 2 and 3, possibly due to the flow motion after plunging wave breaking, as suggested before. It seems indicating that mesh refinement in the wave breaking zone used for mesh 2 and 3 is enough to model correctly the wave breaking and the flow motion above the beach zone. Water height at the crest of the dike shows differences of 55% (mesh 1), 27% (mesh 2), and 7% (mesh 3). Here, only mesh 3 seems to obtain similar results to the finer mesh.

Mean wave overtopping discharge is significantly over-estimated for mesh 1 with difference of 72%, compared with finer mesh. For mesh 2 and 3, the difference with mesh 4 is 26% and 42%, respectively. These results seem to indicate that overtopping, and subsequently the water flow on dike, is sensitive to mesh refinement. Mean wave overtopping discharge depends on water height at the crest of the dike but also on water velocity in water layer. Mesh 4 has a Q between the values obtained by mesh 2 and mesh 3, showing that a convergence is starting after mesh 2. Since mesh 3 gives similar values to mesh 4 for run-up (differences of 2%) and mean wave overtopping discharge convergence is starting after mesh 2, all simulations in Section 5 are conducted using mesh 3.

4.3. DualSPHysics

For DualSPHysics, the sensitivity analysis was performed by changing the initial inter-particle distance, d_p. In this analysis, seven different resolutions were used, with initial inter-particle distance ranging from 1 mm to 20 mm. An active wave absorption system (AWAS) was used to absorb the reflected waves at the piston [31]. Waves were generated following Stokes II and using a piston-type wave-maker. The initial inter-particle distance, total number of particles, and number of particles per wave height are presented in

Table 9. DualSPHysics domain resolutions (with total number of particles per H).

Resolution	dp (m)	Particles	Particles/H
1	0.02	11166	6
2	0.015	19398	8
3	0.01	43745	12
4	0.0075	90093	16
5	0.005	173483	23
6	0.002	1087453	58
7	0.001	4326777	116

.

Run-up case is analysed for all resolutions.

Table 10. Run-up case: mean wave height at gauges located at x = 3.0, 8.0, 11.1, and 12.06 m and the mean and maximum run-up obtained in DualSPHysics.

Resolution	Wave Height (m)			Ru (m)
	x = 3.0 m	x = 8.0 m	x = 11.1 m	Mean	Max
1	0.1086	0.1074	-	0.024	0.032
2	0.1086	0.1165	0.0178	0.039	0.048
3	0.1150	0.1161	0.0459	0.055	0.061
4	0.1144	0.1208	0.0437	0.065	0.069
5	0.1140	0.1220	0.0455	0.074	0.076
6	0.1160	0.1210	0.0510	0.077	0.082
7	0.1125	0.1259	0.0580	0.080	0.083

presents the mean wave height at four gauges located at x = 3.0 and 8.0 m, before the wave breaking, x = 11.1 m, after the wave breaking, and x = 12.06 m at the toe of the dike and the mean and maximum run-up.

For more than 10 particles per wave height, before the wave breaking at gauge x = 3.0 and 8.0 m, the wave height is about the same for the resolution 3 to 7, with maximum difference of about 3% between resolutions 3 and 7. This is in line with [32], who proposed that more than 10 particles per wave height allows for the modelling of wave propagation with good accuracy. However, after the wave breaking at gauge x = 11.1 m, the difference between resolutions 3 and 7 is of 20%. This is caused by the small number of particles in the beach zone, where a noncomplete kernel is inevitably present. Therefore, to capture the complex flow above the beach zone, it is necessary to increase the resolution for 23 particles per wave height or more. For small resolution, the run-up does not reach the gauge x = 11.1 m. As expected for mean and maximum run-up, the results show the convergence of the numerical model when increasing the resolution for 23 particles per wave height or more (resolutions 5 to 7). Considering the low computational runtime, 23 particles per wave height seems to be the more appropriate for the simulation of the run-up case.

For the overtopping case,

Table 11. Overtopping case: mean wave height at gauges located at x = 3.0, 10.0, 12.06 and 12.42 m and the mean wave overtopping discharge obtained in DualSPHysics.

Resolution	Mean Wave Height (m)				Q (10−3 m3/s/m)
	x = 3.0 m	x = 10.0 m	x = 12.06 m	x = 12.42 m
1	0.2250	0.1932	0.0562	-	-
2	0.2259	0.1942	0.0627	-	-
3	0.2386	0.1979	0.089	0.0120	0.76
4	0.2445	0.1992	0.0922	0.0124	0.81
5	0.2458	0.2042	0.0955	0.0150	0.87
6	0.2466	0.2058	0.0978	0.0155	0.88
7	0.2470	0.2058	0.098	0.0160	0.87

presents the mean wave height at four gauges located at x = 3.0 and 10.0 m, before the wave breaking, x = 12.06 m, at the toe of the dike, and 12.42 m, at the top of the dike, and the mean wave overtopping discharge.

The wave height at gauge x = 3.0 and 10.0 m is almost equal for more than 10 particles per wave height, with differences of about 3% comparing resolutions 3 and 7. For gauge x = 12.06 m, these differences of wave height are slightly larger of about 9%. This shows that even 10 particles per wave height allows us to accurately simulate the wave height before the breaking wave, but it is not enough to model correctly the wave breaking and flow motion above the beach zone. For small resolutions (1 and 2), the wave does not reach the gauge x = 12.42 m. Mean wave overtopping discharge depends on water height at the crest of the dike but also on water velocity in water layer. With the increasing of spatial resolution, the flow field is captured more accurately. This confirms that increasing the resolution will improve the accuracy of the wave overtopping discharge. In addition, more than 23 particles (resolutions 5 to 7) per wave height seems to be the more appropriated for the simulation of overtopping case. Based on this sensitivity analysis, all simulations in Section 5 are conducted using 23 particles per wave height (resolution 5).

5. Results and Discussion

Firstly, 2D results obtained by the three numerical models for the nine physical experimental tests are compared, with a deeper analysis for test 1 and 5 used before for sensitivity analysis of each model. The same experimental structure geometry, wave conditions, and water levels were reproduced by the numerical models.

Secondly, 3D versions of DualSPHysics and FLUENT are applied for modelling run-up and overtopping in a 1:15 slope beach protected by a smooth impermeable dike with promenade located between two impermeable groins and the results of both models are compared.

5.1. 2D Modelling

Physical modelling results (see Section 2.2) are here compared with the results of each numerical model for the discretisation choose in the sensitivity analysis (see Section 4).

Results of two models are compared using differences in percentage, calculated as ((mod₁ − mod₂)/(0.5 (mod₁ + mod₂)), where mod₁ and mod₂ refer to results from model 1 and 2, respectively. Results of physical and all numerical models are compared using normalised RMSE, calculated as RMSE = (mean (exp − num)²)^1/2/exp, where num and exp refer to numerical and experimental results, respectively.

The same H and T were given to the wave paddle in the physical model and in the numerical models. However, wave generation varies between models.

A deep analysis is conducted for results of wave height along the flume for the two cases studied in the sensitivity analysis: test 1, for h = 0.50 m, T = 2.04 s, and H = 0.116 m and test 5, with h = 0.58 m, T = 2.04 m, and H = 0.23 s. For the three models and for the physical model, the results presented here correspond to an analysis of 25 waves, i.e., between 29 and 80 s of simulation for cases with T = 2.04 s and between 21.5 and 59 s of simulation for cases with T = 1.5 s. In both tests, wave breaking is plunging type.

Figure 7 presents mean wave height results obtained with the physical model and with the three numerical models for test 1. Wave gauge located in the generation zone is replicated by the numerical models (gauge at x = 3 m). For this first case, with no overtopping, the numerical models have very similar results for x < 10 m, with differences varying from values less than 6% between FLUENT and IH2VOF and values up to 15% between FLUENT and DualSPHysics. Physical model presents a normalised root-mean-square error (RMSE) of 15% for 2 < x < 4 m. At x = 10 m is the one with highest differences between the models, due to breaking behaviours and wave breaking point differences: waves in IH2VOF break after that point and for FLUENT and DualSPHysics waves are breaking. After that, the results become more similar, especially for x = 12 m, at the toe of the dike.

For test 5, Figure 8, where wave overtopping occurs, the numerical models have higher differences in results for x < 8 m than for the test 1, with differences varying from values less than 7% between FLUENT and DualSPHysics and values up to 25% between IH2VOF and DualSPHysics. Normalised RMSE is 36%. IH2VOF presents higher values of RMSE than IH2VOF and lower than FLUENT and DualSPHysics. As for test 1, x = 10 m is the one with higher differences between the models, due to breaking differences: waves in IH2VOF break after that point and for FLUENT and DualSPHysics, waves broke before. After that, the results become more similar, especially for x = 12 m at the toe of the dike. FLUENT and DualSPHysics present very similar results all around the flume.

To analyse the differences observed in the wave height, the time series of surface elevation at section x = 3.0 m are compared in Figure 9 (test 1) and Figure 10 (test 5) and the values of the spectrum for the first, second, and third harmonics are presented in Figure 11 for both tests.

As it can be seen in Figure 9 and Figure 10, the shape of the waves is different in the physical model and in the numerical models. Looking to the spectrum for the three first harmonics, Figure 11, these differences appear clearly. For test 1 and x = 3.0 m, for the first harmonic, the spectrum is higher for DualSPHysics, followed by the physical model when for the second harmonic is the FLUENT that presents highest values, with the others presenting similar values and for the third harmonic, it is the physical models that have highest values. Normalised RMSE is of 21%, 87%, and 89% for first, second, and third harmonic, respectively. For test 5, the behaviour is different. For the first harmonic, physical model has the highest value at x = 3.0 m, IH2VOF the smallest (2 times smaller than experimental value). The normalised RMSE is of 44%, 61%, and 71% for first, second, and third harmonic, respectively. FLUENT has the lowest values for the third harmonic in all cases. IH2VOF has the lowest values for the first harmonic in all cases. Wave characteristics are more different for test 5 than test 1, which seems expected since wave of test 5 presents higher nonlinear characteristics than wave of test 1.

For run-up (Figure 12), mean and maximum values are underestimated by all numerical models. FLUENT presents the higher mean and maximum Ru with differences of 49% and 21%, respectively, compared with physical models. IH2VOF and DualSPHysics show similar smaller values of mean Ru with a difference of 19%. DualSPHysics has the smaller estimation of the maximum Ru.

For wave overtopping discharge (Figure 12), DualSPHysics and FLUENT, which have similar pattern in the wave height along the flume, have differences in mean and maximum wave overtopping discharge greater than the differences present between FLUENT and IH2VOF. Maximum wave overtopping discharge was very similar between FLUENT and IH2VOF, with only 4% of difference, but mean wave overtopping discharge has a larger difference, 23%, and was up to 6.5 times larger than the experimental one. Mean and maximum wave overtopping discharge obtained by FLUENT and IH2VOF was up to 2.7 times greater than those of DualSPHysics, which overestimated the experimental mean wave overtopping discharge by a factor 2.5. These large differences can be due to a different interaction between the water/air mixture flux after breaking and the structure slope due to differences in the breaking position and the type of wave generation. Numerically, wave breaking and overtopping were very sensitive on turbulence models, particularly for plunging wave breaking, water/air mixture flows, and overtopping. For FLUENT, the standard k-ω SST turbulence model, which seems less diffusive than the standard k-ε model [27], was used, and IH2VOF used a k-ε turbulent model on a 2D vertical domain [17,21]. DualSPHysics does not use turbulence model. Energy dissipation due to turbulence during wave breaking can be significantly different and modify the flow dynamics after breaking and the interaction between the water flow and the slope structure.

Analysing the nine tests, Table 1, the results of incident wave height, reflection coefficient, maximum run-up, and mean wave overtopping discharge obtained on the physical model tests are here compared with the results of the numerical models.

The incident wave height values, H_inc, estimated using the method of [38] for all cases, are presented in Figure 13. As can be seen, due to differences in the wave generation, the incident wave height differs from the target one, with the physical model giving higher values and the numerical models, in general, given lower values than the target.

For the physical model, the differences are between 10% and 4%, being lower than the target for all cases 1 to 5 and higher for cases 6 to 9. Numerical models give lower values than the target, except for case 3 for DualSPHysics and case 6 and 7 for FLUENT. The differences are less than 5.8% for DualSPHysics and less than 12.3% for FLUENT and IH2VOF. This is expected since wave generation in DualSPHysics is conducted using a piston-type paddle, as in a physical model.

Figure 14 presents the reflection coefficient, Kr, obtained by each model. Since the structure is impermeable and has the same slope in all cases, the differences in the wave height and Kr could be related with the different wave generation and with the turbulence models used in the numerical models, that lead to differences in the breaking wave and consequent differences in wave dissipation. In general, the values of Kr present significant differences in tests with overtopping. This may be due to differences in predicted overtopping between models, leading to different energy dissipated through this phenomenon and consequently, different reflection. For tests 4, 5, 6, and 9, FLUENT presents Kr higher than all the others, possibly due to less dissipation on the breaking process or differences in wave transmitted by overtopping. Nevertheless, for test 7, FLUENT presents a good accordance with physical model but IH2VOF and DualSPHysics underestimate Kr, possibly due to higher dissipation. FLUENT presents globally the higher estimation of Kr while DualSPHysics and IH2VOF have values of Kr of the same order of magnitude.

Maximum wave run-up and mean wave overtopping discharge obtained in the physical model are compared with those obtained by the numerical model for the nine cases, shown in

Table 12. Maximum run-up calculated by physical and numerical models.

Test ID	h (m)	T (s)	H (m)	Rumax (m)
				Physical Model	DualSPHysics	FLUENT	IH2VOF
1	0.50	2.04	0.116	0.212	0.076	0.171	0.134
3	0.50	1.50	0.230	0.217	0.105	0.179	0.144
4	0.58	2.04	0.116	-	0.178	-	0.152

and

Table 13. Mean wave overtopping discharge calculated by physical and numerical models.

Test ID	h (m)	T (s)	H (m)	Q (L/s/m)
				Physical Model	DualSPHysics	FLUENT	IH2VOF
1	0.50	2.04	0.116	0.00	0.00	0.00	0.00
2	0.50	2.04	0.230	0.00	0.01	0.13	0.01
3	0.50	1.50	0.230	0.00	0.00	0.00	0.00
4	0.58	2.04	0.116	0.02	0.00	0.20	0.00
5	0.58	2.04	0.230	0.35	0.87	2.26	1.80
6	0.58	1.50	0.230	0.05	0.02	0.14	0.05
7	0.62	2.04	0.116	0.35	0.45	0.98	0.51
8	0.62	2.04	0.230	2.84	2.41	5.30	3.95
9	0.62	1.50	0.230	0.98	0.70	1.55	1.58

, for a time interval corresponding to 25 wave period for analysis.

Figure 15 presents the dimensionless maximum run-up, Ru_max/H_inc, obtained for each model and for the physical model tests as function of Iribarren number, ξ, defined as ξ = tan α/(H_inc/L)^1/2, where α is the slope of the smooth dike and L is the wavelength, for the cases 1, 3, and 4. Figure 16 presents the dimensionless mean wave overtopping discharge, Q/(g H_inc³)^0.5, as a function of the freeboard, Rc, divided by H_inc.

Higher values of maximum run-up were obtained for tests 1 to 3 on the physical model and the minimum values for the DualSPHysics model: for test 1, Ru_max varies from 0.08 m to 0.21 m, and for test 3, from 0.11 m to 0.22 m. For test 1, higher Ru^max occurred for the higher H_inc obtained at the physical model. It can be noted that even if the H_inc of the free numerical models is similar, 0.11 m, Ru_max exhibits large variations: 0.08 m for DualSPHysics, 0.13 m for IH2VOF and 0.17 m for FLUENT. A normalised RMSE between 37% and 44% was obtained. Kr for the three models is similar, only slightly higher for FLUENT. For test 3, the maximum value of H_inc was obtained with DualSPHysics and the physical model, were the lower value of Ru_max was obtained at the DualSPHysics and higher value of Ru_max at the physical model. For this test, DualSPHysics presents a slightly higher Kr. Test 4 presents similar behaviour with DualSPHysics and IH2VOF presenting lower values of run-up leading to no overtopping but FLUENT and physical model present higher values of Ru_max with waves overtopping the structure. However, for this case, maximum run-up is higher for IH2VOF than for DualSPHysics. For this test, these two models present almost the same Kr and incident wave height, both lower than FLUENT and physical model values, which can explain the lower run-up. FLUENT presents Kr higher than the other models and a higher transmission by overtopping what might mean that the breaking process was less dissipative.

For mean wave overtopping discharge, as referred before, DualSPHysics and FLUENT, which have similar pattern in the wave height along the flume, have, in general, differences in mean wave overtopping discharge greater than the differences between FLUENT and IH2VOF. FLUENT have generally the higher results and DualSPHysics and physical model have the lower ones. Higher differences between models were found for tests 5 and 8. A normalised RMSE between 51% and 405% for non-zero mean wave overtopping discharge was obtained.

The reasons for the differences found in run-up and mean wave overtopping discharge can be summarised as follows: differences in wave generation and active wave absorption; differences in turbulence models used in the numerical models; numerical diffusion on DualSPHysics due to the noncomplete kernel in the run-up zone where the number of particles is small; viscous dissipation effects of SPH method in solving small water sheet; and the DBC does not include a specific value to define wall friction to correctly solve the flow near the solid walls.

5.2. 3D Modelling

A 3D configuration of a beach protected by a dike constructed between groins is analysed at scale prototype for a regular incident wave (T = 12 s and H = 8.0 m) and wave propagation normal to the coast. Both dike and groins are considered impermeable and located on an impermeable 1:15 foreshore uniform slope, Figure 17. The crest of the 1:2 slope dike is 6.21 m above the still water level, which is at the toe of the structure. The foreshore slope is developed until a water depth 19.22 m and the wave-maker is located at the end of a 127.77 m horizontal bottom. The distance between the wave-maker and the toe of the dike is 416.07 m. The groin is characterised by a 1:2 slope and a length of 105 m between its head (x = 323.49 m) and the crest of the dike (x = 428.49 m). The crest of the groin with 6 m width is composed of a first 1:40 slope between the head and x = 403.49 m and a 1:25 slope between x = 403.49 m and the crest of the dike. The distance between the centreline of two groins is 200 m. The toe of the groin head is at −7.45 m. The wave incidence, normal to the dike, allows limiting the computational domain by two symmetric planes: the centreline plane of groins and the middle plane between two groins.

Within the objective to analyse the effects of the groins in wave overtopping along the dike, the same configuration of dike and foreshore beach without the groins is considered. A 2D computational domain with the same mesh characteristics of the 3D mesh on the symmetric plane of the wave tank is used for this configuration.

Figure 18 shows the position and name of the different gauges used for flow monitorisation. Free surface elevation is monitored using local gauges inside the wave tank: G11 to G14 are located at the beginning of the 1:15 ramp (G10 for 2D configuration); G22, G23, and G24 are placed at the beginning of the groin (G20 for 2D configuration); G32 to G44 are between the head of the groin and the toe of the breakwater (G30 and G40 for 2D configuration); G52, G53, and G54 are at the toe the breakwater (G50 for 2D configuration). Water height and overtopping volume are computed at the head groin (gauge C01) and along the breakwater crest (gauges B01 to B10) using gauges of 3 m wide on the head groin and 10 m wide on the breakwater crest.

Simulation is carried out using RANS-VoF FLUENT and DualSPHysics numerical models each one with its specific characteristics presented in Section 3. The differences between each model are hereby presented and discussed.

In FLUENT, a non-conform mesh is constructed based on hexaedrics control volumes. This approach allows a better optimisation of mesh refinement and a control of the number of control volumes. The mesh refinement criterion respects the results of the sensitive analysis performed in Section 4 at model scale 1:34.5, allowing a good representation of wave propagation, wave breaking, and overtopping in 3D. Far from the groin, the mesh is slightly stretched considering that, in 3D, principal phenomena occur in wave propagation direction.

Figure 19 shows the mesh resolution in several planes the free surface at rest and on the dike and groin structures. The vertical (z direction) resolution in the proximity of the structure is defined as follows. Control volumes with 0.07 m height are used at the dike and groin crest; control volume with 0.28 m height is used at the dike toe and varying from 0.28 m to 0.07 m along the dike slope; control volume height varies along the toe of the groin from 0.28 m, near the dike, to 0.40 m, near the head groin. The horizontal resolution is defined as follows. Control volumes wide along the dike crest, in the y direction, is 0.28 m near the junction between the dike and the groin increasing to 2.0 m near the symmetric plan; control volume length along the groin, in the x direction, is 0.28 m near the junction between the dike and the groin and increases to 0.40 m near the head groin. Mesh resolution in the x direction inside the wave tank is defined by the criterion established in Section 3 for wave propagation, beach, wave breaking, and breakwater zones. Mesh resolution in transversal direction y varies from 0.28 to 2.0 m in the zone between the groins. Mesh refinement in vertical direction in the vicinity of free surface is around 20 control volumes of 0.4 m height between wave through and the crest which is enough in wave propagation zone. Near the dike and the groin, mesh resolution in vertical direction around the free surface is between 0.28 and 0.4 m. The total number of control volumes is 1,454,953.

A Stokes II incident wave is generated using a numerical wave-maker composed of 10 paddles (each with 10 m length) for a better wave absorption of reflected waves by active absorption technique.

CPU time on PC Intel^® Core™ i7-3820 CPU @ 3.60GHz using four cores for parallel computing is 17 h per wave period.

In DualSPHysics, the 3D simulation was carried out with a resolution of 23 particles per wave height, leading to an initial interparticle distance of 0.4 m and total number of particles in the computational domain of 9,907,470 (being 6,910,103 fluid particles and 2,997,367 bound particles). This resolution allows a good representation of dike and groins as well as 3D wave propagation, wave breaking, and overtopping.

The simulation was performed for 600 s of physical time, which means that around 50 waves are generated, on a Nvidia GTX 2080 alongside with an Intel Xeon E5, allowing for large number of particles to be stored. The computational time is approximately 170 h, i.e., 3.4 h per wave period. The AWAS is used to absorb the reflected waves at the piston to mimic the behaviour of an open sea where reflected waves propagate outside the computational domain. Waves were generated following Stokes II and using a piston-type wave-maker.

Figure 20 shows the minimum, mean, and maximum wave overtopping discharge over the dike for both FLUENT and DualSPHysics for configuration with groins (3D case) and without groins (2D case). Figure 21 shows the mean wave overtopping discharge along the dike, in each 10 m wide section, for FLUENT and DualSPHysics. The wave overtopping discharge and the water flow height at the crest of the groin head was also presented for both numerical models in Figure 22, where a large wave overtopping occurs.

Some discrepancies can be observed between FLUENT and DualSPHysics results for 2D and 3D simulations, featuring a much smaller difference of mean wave overtopping discharge in 2D and 3D. Due to the numerical diffusion and viscous dissipation effects of SPH method in solving small water sheet, in some wave periods, DualSPHysics does not present any wave overtopping discharge in 2D and 3D (Figure 20). Nevertheless, the minimum wave overtopping discharge obtain by FLUENT in 2D was very small, only 0.004 m³/s/m, which is not the case in 3D, since FLUENT shows a minimum wave overtopping discharge of 0.020 m³/s/m. For mean wave overtopping discharge, a similar value was obtained by FLUENT and DualSPHysics in 2D, 0.074 and 0.068 m³/s/m, respectively, and maximum wave overtopping discharge exhibits the same order of magnitude. For 3D simulation, similar mean wave overtopping discharge was also obtained by FLUENT and DualSPHysics, 0.036 and 0.033 m³/s/m, respectively, but maximum wave overtopping discharge for FLUENT is half that of DualSPHysics.

It can be noted that both models show the same tendency for mean wave overtopping discharge, which was smaller in 3D than 2D, probably due to the effect of the groin in 3D configuration.

For the 2D model scale simulation (1:34.5) presented Section 4.1, Figure 12, with same incident wave characteristics and slope structure, and only small differences in the bathymetry, FLUENT maximum and mean wave overtopping discharge were 2.6 times larger than DualSPHysics. These results suggest a better agreement between FLUENT and DualSPHysics for prototype scale than scale models, perhaps due to a better accuracy for modelling dissipative structures in turbulent flow and energy dissipation.

The analysis of mean and maximum wave overtopping discharge along the dike, Figure 21, shows large variations along the structure for each numerical model and large discrepancies between them. In section B01, in the axe of the groin, the larger discrepancy between both models can be noted: both mean and maximum wave overtopping discharge obtained by DualSPHysics present a smaller value compared to FLUENT, with 0.003 and 0.046 m³/s/m for mean wave overtopping discharge, respectively. Nevertheless, in the adjacent section B02, both models show a similar trend and the mean wave overtopping discharge has of same order of magnitude. The same observation can be found in sections B09 and B10, were both models show a smaller mean and maximum wave overtopping discharge than sections B03 to B08. In these sections B03 to B08, DualSPHysics presents a similar behaviour for mean wave overtopping discharge, which varies between 0.036 and 0.042 m³/s/m. The same can be noted for maximum wave overtopping discharge. FLUENT shows a more variable behaviour of mean wave overtopping discharge than DualSPHysics, which varies from 0.012 to 0.063 m³/s/m. Nevertheless, as observed before, the global maximum and mean wave overtopping discharge present a similar behaviour even if large differences occur along the dike.

At the groin head, it can be observed a large wave overtopping for both models (Figure 22, Figure 23b and Figure 24b). FLUENT and DualSPHysics present the same order of magnitude of the mean and maximum wave overtopping discharge and of the water flow height at the crest of the groin head, for this very complex free surface flow. Minimum values for wave overtopping discharge and water flow height show larger discrepancies.

Figure 23 and Figure 24 show the velocity intensity at the free surface flow for FLUENT and DualSPHysics, respectively, and the complex 3D interaction between the incident wave and the structure. In general, the 3D wave field in the vicinity of the sea wall shows important differences with a large spatial variation in free surface elevation and velocity intensity. This behaviour is caused by the different nonlinearities of the interactions of the wave with the bottom and by the interactions between the wave and the groin. In FLUENT, the overtopping over the groin becomes more relevant than in the DualSPHysics. It can also be noted that active absorption at the wave-maker was a complex task particularly due to the oblique reflected waves that reach at the wave-maker.

6. Conclusions

This paper presents a study of a run-up and overtopping over an impermeable dike using two mesh-based models, IH2VOF and FLUENT, a mesh-free numerical model, DualSPHysics, and results from 2D physical modelling for a strong energetic incident wave. These waves induce plunging breaking and a chaotic water/air mixture turbulent flow on a small ramp before the water/air mixture run-up/overtopped the smooth dike. These characteristics are challenging for the numerical models and render the differences in wave generation, absorption, and turbulence models more important in the results. In addition, 3D modelling is presented using FLUENT and DualSPHysics for a 1:15 impermeable slope beach protected by a smooth impermeable dike located between two impermeable groins.

For a regular wave target height and period for the same water depth, wave height normalised RMSE of 15% was found. These differences are related to wave generation (experimental and numerical) and to the turbulence models used in the numerical models. Differences in the way the models account for the free surface—the VoF method for FLUENT and IH2VOF and only the water flow for DualSPHysics—lead to differences in breaking position and dissipation. Consequently, differences in run-up, from 94% to DualSPHysics to 21% to FLUENT, and on mean wave overtopping discharge, from 146% to FLUENT (for smaller mean wave overtopping discharge) to 16% to DualSPHysics, were found when compared with the physical model results. A normalised RMSE between 37% and 44% was obtained for maximum run-up and between 51% and 405% for mean wave overtopping discharge. The complexity of the case tested, i.e., plunging wave with intense turbulent flow and a chaotic mixture flow interacting with a dike, highlights the sensibility and limitations of the numerical models for this kind of applications.

Three-dimensional simulation was performed by FLUENT and DualSPHysics for a beach protected by a smooth impermeable dike constructed between impermeable groins. Similar maximum and mean wave overtopping discharges were obtained by FLUENT and DualSPHysics for both configurations. The effect of the groins on overtopping of the dike leads to a significant reduction in mean wave overtopping discharge of about twofold when compared with the same configuration without the groins. Large discrepancies were found on overtopping along the dike, particularly near the junction between the groins and the dike. Nevertheless, mean wave overtopping discharge and water flow height at the crest of the groin head show a similar order of magnitude.

In FLUENT, the overtopping over the groin becomes more relevant than in the DualSPHysics. Besides the differences pointed out in the 2D simulations, 3D active absorption at the wave-maker is a complex task, particularly due to the oblique reflected waves that reach at the wave-maker, and can be one of the causes of the differences found on the results.

The results presented show that even for numerical models that have proved to be in agreement with the experimental data for some wave conditions, as the three models presented here, there are some cases where this is not verified. In this paper, the case analysed is a very complex wave-structure interaction involving plunging breaking followed by an important and intense turbulent air/water mixture flow. Here, the models reach their application limits, resulting in very large differences between the numerical models and the experimental tests.