DV RANS-VOF MODELLING OF DEPTHS AND VELOCITIES OF OVERTOPPING WAVES AT OVERWASHED DIKES

The purpose of this contribution is the representation of real wave overtopping over sea dikes with the Rans-Vof code (IH-2VOF) developed by the University of Cantabria. More specific objectives are: to identify the real capacity of the IH-2VOF model in the prediction of overflow and to determinate the accuracy of these predictions in order to provide designers with a generally applicable methodology to use the code. The model is validated against experimental tests conducted by Hughes at the U.S. Army Engineer Research and Development Center (ERDC). This analysis shows that the model tends to overestimate wave reflection and better results are obtained by introducing a porous layer around the structure as artificial way to reduce reflected wave energy.


INTRODUCTION
Vulnerability and resilience of dikes and sea banks play a key role in the safety of coastal areas.Ideally these defenses should be designed with a crest height and width to prevent flooding for any climate scenarios but the uncertainty in climate conditions, worsened by climate change, and the high costs require that a certain overtopping level has to be accepted.
Dike overtopping produces high-speed currents and turbulence (Schuttrumpf and Oumeraci, 2005) that may induce damage of the protection layers (for instance grass cover, Van der Meer et al., 2010) and expose at erosion the under layers in case of soil dikes.If the overtopping persists for sufficient time, crest lowering and structure breaching may occur (Hughes and Nadal, 2009).Catastrophic consequences of this process happened during the hurricane Katrina in New Orleans (ASCE Hurricane Katrina External Review Panel, 2007).
Accurate estimates of the statistics of overtopping waves in terms of flow depths, duration and especially velocities for a set of climate conditions are needed and have to be combined with consolidated criteria for identifying tolerable overtopping threshold.
The evaluation of wave-structure interaction has been mainly investigated by means of physical model testing in laboratories (a.o.Schuttrumpf and Oumeraci, 2005;Van Gent, 2002) and at prototype scale, achieving in depth knowledge by means of the "Wave Overtopping simulator", Van der Meer et al. (2006).However the theories adopted so far for predicting overtopping depths and velocities along the dike still require analysis and verification (Van der Meer et al., 2010).
Aim of this contribution is to verify the possibility to accurately predict wave overtopping and related fluxes over sea dikes with the RansVof code (IH-2VOF) developed by the University of Cantabria (Lara et al., 2008).The possibility to use a numerical model with high resolution for such purposes would allow obtaining many results in a controlled environment regarding depths and velocities on the dikes, i.e. relevant parameters for a more resilient design.
The paper structure is as follows.First the experimental tests used for comparison with the numerical model are presented.Then the numerical model and its set-up are described, including the validation process to check the model accuracy.The model predictive capacity is then examined by means of comparison with experimental results, such as statistic distribution of wave overtopping depths and velocities (h2%, T2%, u2%) over and landward the dike.A sensitivity analysis is carried out to set the best modeling configuration.Numerical results for the best modeling parameters are finally compared with the existing theoretical approach developed by Schuttumpf et al. (2005).

Description of the tests
The overtopping discharge resulting from combined storm surge overflow and wave overtopping of a levee with a trapezoidal cross section was studied in 1:25 scale in the 0.91-m-wide wave flume at the U.S. Army Engineer Research and Development Center (ERDC), Coastal and Hydraulics Laboratory (CHL) in Vicksburg, MS.
The tested levee cross section replicated in the physical model is shown in Fig. 1 in model-scale units.Fig. 2 shows the relative location of the levee model in the wave flume.The levee crest was approximately 32 m from the wave board (left end of Fig. 2) with the crest elevation 0.61 m above the wave flume bottom.Seaward of the levee model section there was a long 1:100 approach slope and a shorter 1:20 slope transition to the bottom of the flume.Surge and waves that overtopped the levee flowed into a reservoir (right end of Fig. 2), and a pump recirculated the water to the seaward end of the flume.
Four capacitance-type wave gauges were placed at the locations shown in Fig. 2. Gauge 1 was located over the horizontal bottom of the wave flume closest to the wave board.Gauges 2-4 were placed as a three-gauge array near the toe of the levee model berm.Two of the pressure cells were located on the levee crest, and the remaining five instruments were evenly spaced down the protectedside slope.Fig. 3 shows the locations of the pressure cells with spacing given in model-scale units.The purpose of the pressure cells was to measure flow thickness variations over the levee as a function of time.
A fiber-optic laser Doppler velocimeter (LDV) was used to measure the horizontal component of flow velocity directly above pressure cell P2 (see Fig. 3) near the rearward edge of the levee crest.This cross-flume location was directly above pressure cell P2.For each experiment the vertical position of the LDV beam crossing was adjusted to an elevation approximately half of the water depth of the steady overflow.This vertical position was thought to provide a reasonable value for depth-averaged horizontal velocity, but the drawback to this vertical location was the loss of velocity signal during wave troughs when the water level fell below the elevation of the laser beams.

Simulated tests
The numerical analysis is carried out for four different tests that are characterized by different wave height, peak periods and structure submergence (see Table 1).The work is organized into two different phases: one to check the precision of the model and the other to verify the model predictive capacity.

DESCRIPTION OF THE NUMERICAL MODEL
Numerical models of fluid/wave-structure interactions are increasingly becoming a viable tool in furthering our understanding of the complicated phenomena that govern the hydraulic response of breakwaters, including effects of permeability (Losada, 2003).These models are based on different approaches, for instance: Lagrangian models with particle-based approaches such as the Moving Particle Semi-Implicit method (Koshizuka et al., 2005) and Smooth Particle Hydrodynamics (SPH) (Dalrymple et al., 2009), Reynold Averaged Navier Stokes-Volume Of Fluid (RANS-VOF) models such as those developed by Lara et al. (2008) and Shi et al. (2004).
The IH-2VOF model (Lara et al. 2008) has been developed by implementing various extensions to the RIPPLE model (Kothe et al, 1991; originally designed to provide a solution of two-dimensional versions of the Navier-Stokes equations in a vertical plane with a free surface), making it specifically applicable to the study of wave interaction with coastal structures.
The model solves the 2DV Reynolds averaged Navier-Stokes (RANS) equations based on the decomposition of the instantaneous velocity and pressure fields into mean and turbulent components, and the   k equations for the turbulent kinetic energy k and its dissipation rate  .The influence of turbulence fluctuations on the mean flow field is represented by the Reynolds stresses.A nonlinear algebraic Reynolds stress model is used to relate the Reynolds stress tensor and the strain rate of mean flow.The free surface movement is tracked by the volume of fluid (VOF) method, for only one phase, water and void.
In order to replicate solid bodies immersed in the mesh, instead of treating them as sawtooth-shape, the model uses a cutting cell method first presented by Clarke et al. 1986.This technique uses an orthogonal structured mesh in the simulations to save computational cost.This approach defines the openness function  to define the fraction of volume of free space in the cell.According to this definition 0   is a "solid cell" (entirely occupied by the solid), Variables in the cell or on the cell faces are redefined by the product of the openness coefficients times the original variables.To numerically implement this, partial cell coefficients need to be defined at both the cell centers and boundaries: , which correspond to the openness at the left, right, top and bottom, respectively.As a result, the governing RANS equations are redefined as follows: The last term in equation ( 2) is a virtual force.It appears because the IH-2VOF model considers two different numerical techniques to simulate moving bodies within the computational domain.The first one is called "virtual force method" (see Mittal and Iaccarino, 2005) and the other one is the "direct forcing method" (Mohd-Yusof, 1997).

Model set-up
In these preliminary tests the IH-2VOF model is used forcing as offshore boundary condition the same water levels measured by the gauge in front of the structure (at around 30 m from the wavemaker).
 Seaward boundary condition measured water level at the laboratory gauge P2  Landward boundary condition absorption The length of the numerical channel is therefore shorter than the one in the laboratory since the numerical simulation starts from Wave Gauge (WG) 2 on (placed at a distance of 29.7151 m from the WM).
 Flume dimension 5.7 m x 0.8 m  Mesh resolution 0.01 m x 0.005 m

Results
With this type of approach the representation of water levels along the channel is very good.In Fig. 4 three graphs are reported: the first represents the water level at the WG in front of the structure (WG2), the second at the WG sets on the crests of the structure at the offshore edge (P1) and the last at the WG sets at inshore edge (P2).At the WG 2 the numerical water level (red line) reproduces almost perfectly the laboratory one (blue line) because it is precisely the boundary condition imposed.Fig. 5 shows a comparison among measured and simulated overtopping volumes at gauge P1 (gauge at the crest offshore edge).It can be observed that the representation is overall optimal and both experimental and numerical data are very well approximated by a Weibull statistic distribution.
In Table 5 (case "A") the experimental and numerical significant results are compared in terms of the wave height and peak period at the gauges WG1 and WG2 (i.e. the gauges between the wave maker and the structure).The more significant quantitative results (h2%, T2%, u2%) are reported in Table 6 (case "A").The scatter of the numerical results from the experimental data is very modest, suggesting that the RansVof model ensures a good precision in the representation of the overtopping process.

Model prediction capacity
In order to check the predictive capacity of the model, the numerical channel is identical to the experimental one and a Jonswap spectrum characterized by the same significant wave height, peak period and spectrum parameters as in the laboratory is imposed as offshore boundary condition.Therefore the following model settings are adopted: As before, in Table 5 (case "B") it is possible to see the comparison between the experimental and numerical significant wave height and peak period at the gauges WG1 and WG2.
The results of u2%, h2% and T2% for the tests R14, R18 and R20 (characterized by wave height, peak period, water deep and submergence shown in Table 1) are reported in Table 6 (case "B"): all the parameter values derived from the experiments are significantly higher than the numerical values.
It is possible to conclude that in these tests the numerical approach cannot capture the higher overtopping volumes.This poor approximation may be essentially justified by the lower incident wave height simulated in front of the structure (see Table 5, case "B").Several simulations are then conducted by changing the discretization of the domain.The best solution seems to be achieved when the resolution of the numerical mesh equals the 1% of Hs in the horizontal direction and the 0.5% of Hs in the vertical direction.Even though this expedient, the variation of the results is modest.
To improve the quality of the simulations, three different hypotheses can be made: 1. too high dissipation induced by the selected turbulence model; 2. possible wave decay along the numerical flume; 3. too high reflection induced by the dike-obstacle.
The results of a sensitivity analysis to check these three hypotheses and define a method to improve the performance of the numerical model are reported in the next sections.
The simulations are performed first without including the turbulence model, then in a shorter flume (that is around 3 wave lenghts long in front of the structure); finally the structure is modifyied by including a composite porous medium.

Simulations without turbulence
To verify the dissipation induced by the turbulence model used in the RansVof code (   k turbulence model), some simulations are repeated without turbulence.
The same numerical flume described for case "B" is here used again (identical to the experimental one).Also boundary conditions and mesh resolution are kept unchanged.
Without turbulence (see Table 6, case "C") the underestimation of the overtopped volumes and consequent overestimation of the characteristic parameters h2%, T2%% and u2% are reduced but not in a satisfactory manner.In fact the incident wave height simulated in front of the structure still remains lower that the experimental one (see Table 5, case "C").

Simulations with a shorter channel
If the problem is that a wave decay occurs along the channel, a way to reduce this decay could be the use of a numerical channel shorter than the experimental one.Some simulations are repeated in a channel whose length from the wave-maker to the structure approximately equals three times the maximum wave length.The following settings are therefore adopted in this case: The results of tests R14, R18 and R20 (characterized by wave height, peak period, water depth and structure submergence in Table 1) are shown in Table 6 (case "D").The representation of overtopped volumes does not improve significantly, the overestimation of the incident wave height being still high (Table 5, case "D").
Therefore this model setting does not allow to reach a final conclusion regarding the improvement of the model predictive capacity.However it can be concluded that shortening the channel up to 3 times the peak wave length does not induce any relevant change of the results, suggesting that it would be possible to reduce in this way the computational effort without loosing accuracy.

Simulations without the structure
To check if the discrepancy with experiments is due to wave reflection or to an internal numerical dissipation, a simulation without the structure is carried out.
The simulation is therefore performed in the same conditions as in case "B" (i.e.R14 with a numerical channel perfectly identical to the experimental one and without turbulence):  Mesh resolution cell width: 0.01÷0.02m cell height: 0.01 m The results obtained at the gauges set between the numerical wave maker and the structure are reported in Table 2:  As it can be seen in Fig. 7, the computed and measured wave statistics at P1 perfectly match.This means that the problems with mismatching and dissipation should be related to wave reflection from the dike seaward slope or eventually to the way the run-up process is reproduced.

Simulations with a part of porous medium
In order to reduce wave reflection, the representation of the structure is changed from perfectly impermeable to partially permeable.Different configurations are tested (see Fig. 8): 1. the seaward slope is replaced by a homogeneous porous medium (case "E") 2. the seaward slope is replaced by a composite porous medium (case "F") 3. a thin porous layer is wrapped around the structure (case "G") As a first attempt (case "E") a single porous layer at the seaward-side is inserted.Then, as a second attempt (case "F"), the seaward slope is replaced by three layers each characterized by different value of porosity and nominal diameter.This approach allows to create a gradient of porosity: the lowest layer presents a greater porosity than the upper layers, in order to reduce both wave reflection in the lower part and wave percolation in the upper part (i.e. to maintain wave run-up).Finally (case "G") to avoid the discontinuity that alters the flow, the best solution seems to include a thin and homogenous layer of porous medium around the entire structure (seaward side, crest and landward side).The characteristics of each porous layers for all the cases are shown in Table 3.The computed and measured wave statistics at P1 are reported in Fig. 9 and Fig. 10 for the cases "E" and "F" respectively.
In Table 4 it is reported a comparison between u2% obtained at the gauge P2 with the configuration "E" and "F" (Test R14).The use of the porous medium as seaward slope of the structure reduces the wave reflection but accelerates the flow on the crest due to the discontinuity between the permeable and impermeable parts of the structure.The use of a permeable seaward slope generates a discontinuity between the porous seaward side and the impermeable crest that alters the water flow.So in case "G" the porous medium is set around the entire structure.This thin layer of porous medium (about 1/10 Hs), is characterized (see Table 3) by a low porosity in order to limit the water percolation into the dike.

CASE "E" CASE "F" CASE "G"
Hence to avoid structure discontinuity that alters the flow, the best solution seems to include a thin and homogenous layer of porous medium around the entire structure (seaward side, crest and landward side).Keeping constant all the other parameters, the results show a significant improvement in the performance of the model in reproducing the statistical distribution of the overflowing volumes in correspondence of the gauge P1 at the offshore edge.

Verification of the reflection coefficient
The reflection coefficient obtained by the simulations is compared in Fig. 12 with the data for smooth straight slopes that are included in the reflection database by Zanuttigh and Van der Meer (2006).It can be observed that the value (cyan filled-in circle) obtained by forcing the model exactly with the measured levels in front of the structure (case "A") falls perfectly in the range of the experimental values as well as the values (red filled-in circles) derived with the inclusion of the porous layer (case "G") and the use of the target wave spectrum as boundary conditions.The values obtained by the model in predictive set-up when the structure is fully impermeable (case "B") are well above the experimental data, being around two times the measured values for the same Iribarren-Battjes parameter.
Therefore the overestimation of wave reflection with the numerical code derived from the previous analysis is confirmed.

Relevance of profile versus local measurement
Fig. 13 reports, in the same graph, the computed and the measured velocities for the case "G".The figure shows the velocity profile obtained with the introduction of the porous medium.The red, green and blue colours represent respectively the maximum, the minimum and the average velocity.The velocity values (maximum, minimum and averaged) at the measurement point are well approximated.
The model allows to get the velocity profiles over the flow depth in time -usually not measured in the lab or at prototype scale.Peak velocities along the profile differ in some cases of a factor 2 from the u2% at the measurement point, highlighting the relevance of this modelling result for design purposes.14 shows the good agreement among the points extracted from the numerical simulations and the theoretical curve.Table 7 reports the values of rmse and the ratio between the standard deviation relative to the numerical curve and the standard deviation relative to the theoretical curve (see line "evolution of layer thickness along the crest").

Overtopping velocities on the dike crest
The variation of the overtopping velocities on the dike crest is derived from the Navier-Stokes equation and reads = overtopping velocity at the beginning of the dike crest.The friction coefficient must be determined empirically, here the value used in Fig. 15 is 0.2.Table 7 shows the values of rmse and the ratio between the standard deviation relative to the numerical curve and the standard deviation relative to the theoretical curve (see line "evolution of velocity along the crest").

Flow field on the landward slope of the dike
Schuttrumpf developed an approach to describe overtopping velocities on the landward slope and, as for the evolution of the velocity and the layer thickness along the crest, verified this theoretical approach by the small scale model test results.Fig. 16 shows the computed and theoretical evolution of the velocities on the landward slope.Table 7 shows the values of rmse and the ratio between the standard deviation relative to the numerical curve and the standard deviation relative to the theoretical curve (see line "evolution of velocity on the landward side").The model reproduces with high accuracy the physical process of wave overtopping when it is forced with the measured levels obtained in the channel.
When used for predicting wave overtopping based on a Jonswap target wave spectrum, the model does not fully represent the overtopping volume distribution and tends to underestimate higher volumes, especially for low waves.This can be only partially ascribed to the mesh resolution and turbulent dissipation.Specifically best results are obtained with mesh resolution: 1% Hs horizontal and 0.5% Hs in vertical without turbulence model.Actually the discrepancy among numerical and experimental distribution of overtopping volumes and overtopping statistic parameters (h2%, u2%, T2%) can be fully explained by the model overestimation of wave reflection.
The inclusion of a porous layer around the structure reduces wave reflection and therefore significantly improves -keeping constant the other parameters -the performance of the model in reproducing wave overtopping statistics together with flow characteristics at the dike landward slope.To avoid discontinuity in water fluxes and contemporarily reduce wave reflection from the dike, the best solution seems to include a homogenous thin porous layer around the dike characterised by low thickness (around 1/10 Hs) and low porosity (n=0.1 to avoid percolation).
The model represents very well u2% at the measurement point over the crest and at the same time allows to get the velocity profiles over the flow depth in time -usually not measured in the lab or at prototype scale.Peak velocities along the profile differ in some cases of a factor 2 from the u2% at the measurement point, highlighting the relevance of this modelling result for design purposes.
The numerical results also fit very well the theoretical approach developed by Schuttrumpf and Oumeraci (2005) to represent flow depths and velocities over and landward the dike.

Figure 2 .
Figure 2. Profile view of wave flume (model-scale units).The wave-marker is to the left side of the figure.

Figure 3 .
Figure 3. Locations of flush-mounted pressure cell on levee (model-scale units).

Figure 4 .
Figure 4. Trend of three water levels measured respectively at WG 2 (placed in front of the structure), at the pressure gauge P1 (placed on the crest of the structure) and at the pressure gauge P6 (placed on the land side of the structure).


Fig.6shows a comparison among the overtopped volumes measured in the laboratory (blue color) and the numerical ones (red color) for the test R14 at the gauge P1.It is possible observe clearly that the model underestimates the overtopped volumes especially in correspondence of the bigger volumes.As before, in Table5(case "B") it is possible to see the comparison between the experimental and numerical significant wave height and peak period at the gauges WG1 and WG2.The results of u2%, h2% and T2% for the tests R14, R18 and R20 (characterized by wave height, peak period, water deep and submergence shown in Table1) are reported in Table6(case "B"): all the parameter values derived from the experiments are significantly higher than the numerical values.It is possible to conclude that in these tests the numerical approach cannot capture the higher overtopping volumes.This poor approximation may be essentially justified by the lower incident wave height simulated in front of the structure (see Table5, case "B").


seaward boundary condition Hs and Tp (Jonswap spectrum)  landward boundary condition absorption  flume dimension 14 m x 0.8 m  mesh resolution cell width: 0.001÷0.0011m cell height: 0.005 m

Figure 7 .
Figure 7. Test R14 (test without structure): comparison of computed and measured overtopping volumes at gauge P1.

Figure 8 .
Figure 8. Geometrical scheme for the test with a single porous layer.

Figure 12 .
Figure 12.Comparison between the reflection coefficient obtained from test "A" (cyan filled in circle), test "B" (red empty in circles), test "G" (red filled in circles) and Zanuttigh and Van der Meer's database (2006).
Figure 13.Test R14 (case "G"): comparison of computed and measured horizontal velocity at gauge P2 -test with porous medium.
thickness at the beginning of the landward slope; B v = overtopping velocity at the beginning of the landward slope; f = bottom friction coefficient;  = landward slope angle.

Figure 15 .
Figure 15.Evolution of the overtopping velocities on the dike crest.

Figure 16 .
Figure 16.Evolution of the overtopping velocities on the landward side.

AND THEORY (SCHUTTRUMPF) Layer thickness on the crest of the dike
On the dike crest, the layer thickness c h depends on the width of the dike crest B and the x-

Table 7 . Values of rmse and ration between the standard deviation relative to the numerical curve and the standard deviation relative to the theoretical curve for the case "G".
Numerical simulations with the RansVof model were performed and compared with the tests provided by Hughes on impermeable submerged or zero freeboard dike.