Failure of Grass Covered Flood Defences with Roads on Top Due to Wave Overtopping : A Probabilistic Assessment Method

Hard structures, i.e., roads, are commonly found over flood defences, such as dikes, in order to ensure access and connectivity between flood protected areas. Several climate change future scenario studies have concluded that flood defences will be required to withstand more severe storms than the ones used for their original design. Therefore, this paper presents a probabilistic methodology to assess the effect of a road on top of a dike: it gives the failure probability of the grass cover due to wave overtopping over a wide range of design storms. The methodology was developed by building two different dike configurations in computational fluid dynamics Navier–Stokes solution software; one with a road on top and one without a road. Both models were validated with experimental data collected from field-scale experiments. Later, both models were used to produce data sets for training simpler and faster emulators. These emulators were coupled to a simplified erosion model which allowed testing storm scenarios which resulted in local scouring conditioned statistical failure probabilities. From these results it was estimated that the dike with a road has higher probabilities (5 × 10−5 > Pf >1 × 10−4) of failure than a dike without a road (Pf < 1× 10−6) if realistic grass quality spatial distributions were assumed. The coupled emulator-erosion model was able to yield realistic probabilities, given all the uncertainties in the modelling process and it seems to be a promising tool for quantifying grass cover erosion failure.


Introduction
Structural embedment of roads on top of dikes will generate stability effects in the flood defence during normal operation and during a flood event [1].Yet, these effects are not explicitly considered in the current probabilistic assessment methods.It is expected that climate change increases sea level rise and that it will affect the wave climate around the world [2,3].For rivers, extreme water levels occur more frequently as shown for the River Rhine basin [4], for example.Such climate change effects increase the risk of flooding of low lying, highly populated areas around world [5,6].In the last decades, flood defence design has moved towards a risk-based approach [7,8] in which it is possible to express the design parameters uncertainty as probabilistic distributions [8].Such methods allow inclusion of the expected future effects due to climate change as modifications in the parameter statistical distributions.Nevertheless, there is a lack of methods for assessing non-convectional flood defences so that their interaction effects are also captured and reflected in the estimated failure probabilities [1].
In case of grass covered dikes the failure mechanism known as wave overtopping is influenced by the presence of non-water retaining functions [1].This failure mechanism consists in the intermittent surpass of water excess volumes over the dike derived from the 'attacking' waves during a storm.Only the waves that surpass the crest height will overtop and 'attack' the landward slope of the dike, potentially eroding the cover until an eventual dike breach.
For dike cover erosion assessment, reliability methods based on partial safety factors validated by a fully probabilistic Monte Carlo method are available [9].These methods are based on the classical overtopping approach of defining dike safety based on an allowable overtopping discharge q max [m 3 /s/m].Other studies concluded that overtopping erosion is better estimated from individual wave overtopping volumes V [m 3 /m] [10][11][12][13].In Ref. [14], three different erosion predictors (velocity excess, shear stress excess, and work excess) were tested for calculating the scour depths per wave volume and they found that work excess was the best predictor.Simultaneously, the cumulative overload method [15] was developed which relies on the shear stress principle, using the total overtopping duration per wave as input while acknowledging that the "real" erosion time per wave may be shorter.Additional research [16,17] allows inclusion of the effects of obstacles and transitions in both the hydraulic load and resistance terms in the form of calibration coefficients for the cumulative overload method from Ref. [15].
With respect to the grass cover resistance, a conceptual geo-mechanical approach [18] describes the physical process of grass cover erosion based on the vertical and horizontal stresses that act upon a turf element model.In addition, two main studies [19,20] recommended grass and soil quality values expressed in terms of their erodibility rates.Based on these studies and concepts, Ref. [21] aimed to determine the grass cover failure of conventional dikes in terms of probabilities, while considering the spatial distribution of grass quality.The method used in this last study is fully probabilistic and could be adapted to more complex dike profiles.Nevertheless, it relies on empirical formulas for estimating the water depths and flow velocities on top of the crest and along the landside slope.These empirical formulas [22] simplified the estimation of water depth and flow velocity for the design of sea dikes by proving a set of equations.More holistic approaches [23,24] are available for including the effects of bottom irregularities in the flow.Both studies used one dimensional (1D) models, which allowed to simulate the time-dependent erosion process of a grass covered coastal dike.Yet, both models estimate the exerted bottom shear stresses as a function of uniform flow models, such as Chezy or Manning, which are only applicable for smooth bed slope transitions.If abrupt bottom slope variations occur, the vertical velocity distribution tends to have a larger variability [25], thereby violating the uniform flow assumption.This variability might result in more turbulent flows, higher energy dissipation and hence less accurate bed shear stress estimations.The work presented in Ref. [26] used a more detailed computational fluid dynamics (CFD) model and showed that the presence of a road leads to higher scour depths in places with little to non-existent grass cover.In this paper, we extend the deterministic CFD approach of Ref. [26] towards a fully probabilistic assessment method for failure of the grass cover due to wave overtopping erosion.
Emulation of detailed models is as an efficient way to reduce computational times, for implementation in probabilistic failure studies.Emulation consists of building a computationally inexpensive model and training it with data sets generated from the input-output data sets obtained from more complex models [27,28].Hydrodynamic model emulation has proven to be a powerful tool for water level predictions while improving the calculation speed significantly [29,30].Also, for the reliability assessment of flood defences, emulation has already been implemented for other failure mechanisms, such as backward piping erosion [31][32][33].For wave overtopping, [34] used the results of 10,000 physical model tests of different types of coastal defences to train a neural network.
The aim of this study was to develop an emulator-based dike failure assessment methodology for estimating the effects of a road over the crest including the spatial variation of the grass quality in the failure probability of the grass cover.This failure is estimated as a function of a storm surge wave overtopping event composed of a set of single overtopped water volumes due to wind generated gravity waves, which flow over the dike profile.The methodology was based on two different CFD transient models; one with a road on top [35] and one fully covered with grass.The research question we are answering in this paper is: What is the effect of a road on the grass cover failure probability of a dike due to overtopping waves?
The article is composed of 6 main sections.Section 2 explains the theoretical background of the shear stress excess (T).Section 3 explains the experiment in Millingen from which the collected data was used for construction and calibration of the models.Section 4 presents a detailed explanation of the steps of the developed methodology and its correspondent assumptions.Sections 5 and 6 present the results and discussion and Section 7 gives the conclusions to highlight the major findings.

Theoretical Background
Wave overtopping erosion is a transient process in which tensile and shear stresses are exerted along the dike grass cover over time because of the overtopped water volume's momentum.While the erosion is generated by one single wave, it has been observed that failure of the dike cover is mostly achieved by the overtopping of numerous waves during a storm event [12].The time span in which erosion caused by a single wave takes place at a specific location is in the order of seconds [26] and it is termed total wave overtopping time (t total ).This duration differs from wave to wave depending on its volume, flow depth, and flow velocity and it is spatially influenced by the surface irregularities and bottom slope.If the generated stresses exceed the mechanical resistance of the grass cover elements (e.g., grass leaves, roots, substrate), erosion occurs by pulling, scouring, and flushing the eroded material downstream.

Shear Stress Excess (T)
Most of the overtopping scour erosion methods are based on definitions of critical threshold values for average overtopping discharges in steady flow conditions [14].For combined grass and soil covers, the scour threshold can also be defined by a critical shear stress τ c , which depends on the grass cover quality [19].This value only represents the threshold which defines if scouring occurs or not.To include two other important scour characteristics (erosion rate and real erosion duration), the time dependent shear stress excess [N•s/m 2 ] concept from [36] was adopted in the present study.It represents the surplus of shear stress during the period of time in which erosion occurs.Only during this period does erosion take place.This implies that for larger values of τ c , the difference between the total overtopping duration (t total ) and the real erosion time span (e.g., f (τ c2 ) = t 4 − t 2 in Figure 1) is significant.The value of T is calculated for a given location along the profile as the integral of the bottom shear stress function τ(t) minus the critical erosion stress threshold τ c , over the erosion time span (e.g., area A for time span t 4 − t 2 and τ c2 in Figure 1).A lower value of τ c represents lower resistance to erosion and implies a longer time span of erosion and a greater.
research question we are answering in this paper is: What is the effect of a road on the grass cover failure probability of a dike due to overtopping waves?
The article is composed of 6 main sections.Section 2 explains the theoretical background of the shear stress excess ( ).Section 3 explains the experiment in Millingen from which the collected data was used for construction and calibration of the models.Section 4 presents a detailed explanation of the steps of the developed methodology and its correspondent assumptions.Sections 5 and 6 present the results and discussion and Section 7 gives the conclusions to highlight the major findings.

Theoretical Background
Wave overtopping erosion is a transient process in which tensile and shear stresses are exerted along the dike grass cover over time because of the overtopped water volume's momentum.While the erosion is generated by one single wave, it has been observed that failure of the dike cover is mostly achieved by the overtopping of numerous waves during a storm event [12].The time span in which erosion caused by a single wave takes place at a specific location is in the order of seconds [26] and it is termed total wave overtopping time ( ).This duration differs from wave to wave depending on its volume, flow depth, and flow velocity and it is spatially influenced by the surface irregularities and bottom slope.If the generated stresses exceed the mechanical resistance of the grass cover elements (e.g., grass leaves, roots, substrate), erosion occurs by pulling, scouring, and flushing the eroded material downstream.

Shear Stress Excess ( )
Most of the overtopping scour erosion methods are based on definitions of critical threshold values for average overtopping discharges in steady flow conditions [14].For combined grass and soil covers, the scour threshold can also be defined by a critical shear stress , which depends on the grass cover quality [19].This value only represents the threshold which defines if scouring occurs or not.To include two other important scour characteristics (erosion rate and real erosion duration), the time dependent shear stress excess [N•s/m 2 ] concept from [36] was adopted in the present study.It represents the surplus of shear stress during the period of time in which erosion occurs.Only during this period does erosion take place.This implies that for larger values of , the difference between the total overtopping duration ( ) and the real erosion time span (e.g., ( ) = − in Figure 1) is significant.The value of is calculated for a given location along the profile as the integral of the bottom shear stress function ( ) minus the critical erosion stress threshold , over the erosion time span (e.g., area A for time span − and in Figure 1).A lower value of represents lower resistance to erosion and implies a longer time span of erosion and a greater.The fact that the erosion time span (ti+1 − ti) and the total wave overtopping duration time ( ) differs may have a significant influence on the estimated erosion was acknowledged in Ref. [37], who included the "real" erosion excess time by assuming a triangular shape of the wave overtopping The fact that the erosion time span (t i+1 − t i ) and the total wave overtopping duration time (t total ) differs may have a significant influence on the estimated erosion was acknowledged in Ref. [37], who included the "real" erosion excess time by assuming a triangular shape of the wave overtopping discharge hydrograph.However, this assumption may be too general for representing the shear stress excess over a highly irregular bottom profile.The overtopping integration bounds (t i+1 and t i ) vary between locations due to the presence of bottom irregularities which may either accelerate or decelerate the overtopping flows.In addition, their values are also determined by the value of τ c (See Figure 1).

Erosion Model as a Function of T
The erosion model used in the present study is derived from the time dependent mass erosion rate equation for cohesive soils presented in Ref. [38] as: In which M p [kg/(m 2 s)] corresponds to the characteristic soil mass transport coefficient, dm dt corresponds to the mass transport rate in time, τ [N/m 2 ] corresponds to the exerted bottom shear stress, and τ c [N/m 2 ] corresponds to the critical shear stress threshold.
According to [36], an overall erosional strength parameter C E [m/s −1 ] for soil and grass together was derived [19] as: Note that the C E parameter is directly proportional to M p and inversely proportional to τ c .Based on this equivalence and by dividing Equation ( 1) by the dike cover relative density per unit of width (ρ cover [kg/m 2 ]), the erosion rate at any specific location for a unitary width of dike can be expressed as: Since C E or ρ cover are independent of time, the scour depth (ε) in a particular location can be calculated by integrating Equation (3) as: The integral on the right side of Equation ( 4) represents the T value for a given τ c for a single overtopped wave.The integral bounds of the right-side integral are defined by the specific moments in time where τ(t) ≥ τ c , as erosion will occur during this condition (Figure 1).The integrand bound x i represents the initial bottom elevation before the scouring process and x i+1 the resultant bottom elevation after the scouring process.This model includes the time-dependent variability of the shear stress, but it only estimates the erosion process as a function of the bottom shear stresses.The vertical tensile resistance of grass and roots of the grass cover [36] are accounted for in the parameterization by the C E coefficient.These tensile stresses are described by the turf element method [36] but were too complex to be included in the model.Therefore, the shear stress excess erosion concept presented in Section 2.1 is selected as it requires less input data and allows to use the erodibility value C E [19].This simplification makes the implementation more feasible for a probabilistic model as it represents the grass quality and resistance uncertainty in one single stochastic variable.

Millingen aan de Rijn Wave Overtopping Experiment with a Road on Top
Experiments for wave overtopping were conducted on top of a riverine dike along the river Waal nearby Millingen aan de Rijn (The Netherlands) in February and March of 2013 [39].Random wave volumes V [m 3 /m] were released by the wave overtopping simulator (WOS) [15] during a period of time equivalent to the duration of a storm.The WOS was located on top of the crest of the dike close to the riverside vertex (see Figure 2).For the present study, a storm condition is defined by a crest height ( [m]), a water level and the total number of waves [-] generated in the foreshore, with significant wave height, [m], during the storm duration period [s].In addition, a storm condition is often characterized by its mean overtopping discharge per unit of width [m 3 /s/m].This discharge is calculated as the total overtopped volume of water per storm divided by storm duration as shown in Equation ( 5): Note that the same storm event may generate a different depending on the crest level of the dike as only a limited number of wave volumes will overtop the dike [-].Hence, for this study, a combination of factors (e.g., wave conditions, dike geometry, and roughness) that originate a certain is referred to a storm condition.In the WOS experiment, 10 storm conditions were simulated consecutively.After each simulated storm condition, the dike was scanned using a 3D laser scanner to estimate the scour depths and spatial scouring patterns.This experiment included a road on top of the dike to analyze the effect of transitions on the scouring process due to wave overtopping.The experiment was divided into two parts performed in two adjacent locations.Each test location was 4 m wide with lateral walls that prevented leakage during the tests.

Millingen Experiment Part I: Scour Measurements
During the first part of the experiment, the WOS was located near the riverside edge of the dike (zone A, Figure 2).From this location, a series of volumes were released during a period of six hours, which is equivalent to the duration of the design storm.The released volumes were randomly sampled following a Weibull distribution which represents the stochastic nature of the waves that overtop the dike during a storm.The former Dutch safety standards, for example [11], allow storm conditions of maximum overtopping discharge between 0.001 m 3 /s/m and 0.01 m 3 /s/m depending on the dike location, cross sectional design, and materials [12].The first part of the experiment [39], was conducted for different consecutive storm conditions ( ) as presented in Table 1.For the present study, a storm condition is defined by a crest height (R c [m]), a water level and the total number of waves N w [-] generated in the foreshore, with significant wave height, H s [m], during the storm duration period D [s].In addition, a storm condition is often characterized by its mean overtopping discharge per unit of width q m [m 3 /s/m].This discharge is calculated as the total overtopped volume of water per storm divided by storm duration as shown in Equation ( 5): Note that the same storm event may generate a different q m depending on the crest level of the dike as only a limited number of wave volumes will overtop the dike N ow [-].Hence, for this study, a combination of factors (e.g., wave conditions, dike geometry, and roughness) that originate a certain q m is referred to a storm condition.In the WOS experiment, 10 storm conditions were simulated consecutively.After each simulated storm condition, the dike was scanned using a 3D laser scanner to estimate the scour depths and spatial scouring patterns.This experiment included a road on top of the dike to analyze the effect of transitions on the scouring process due to wave overtopping.The experiment was divided into two parts performed in two adjacent locations.Each test location was 4 m wide with lateral walls that prevented leakage during the tests.

Millingen Experiment Part I: Scour Measurements
During the first part of the experiment, the WOS was located near the riverside edge of the dike (zone A, Figure 2).From this location, a series of volumes were released during a period of six hours, which is equivalent to the duration of the design storm.The released volumes were randomly sampled following a Weibull distribution which represents the stochastic nature of the waves that overtop the dike during a storm.The former Dutch safety standards, for example [11], allow storm conditions of maximum overtopping discharge q max between 0.001 m 3 /s/m and 0.01 m 3 /s/m depending on the dike location, cross sectional design, and materials [12].The first part of the experiment [39], was conducted for different consecutive storm conditions (q m ) as presented in Table 1.While the experiment's durations range from 0 to 3 h for each mean discharge, dikes are often designed for storms that range between 2 and 6 h depending on their boundary conditions.High resolution 3D laser scan images were taken (Figure 2) every 2 h with an accuracy of 2 mm.For some storm conditions, it was necessary to accelerate or decelerate some of the tests due to the WOS pumping constraints and release capacity restrictions.The different storm average discharges are conditioned by the relative position of the crest of the dike with respect to the still water level, the dike geometry, dike cover, and revetment type.
In the initial state profile, already significant damage due to traffic was observed in the transition gap (zone C, Figure 2).This gap increased and propagated landward with each test, as bare soil erodes faster than grass covered soil.Besides this zone, no significant erosion was observed in the grass cover for the 0.001 and 0.010 m 3 /s/m tests.The scan labels presented in Table 1 include the average overtopping discharge and the scan number of that experiment.

Millingen Experiment Part II: Flow Depths and Velocity Measurements
The second part of the experiment consisted of measuring flow depth time series and velocity time series of individual overtopping wave volumes along the dike profile.Paddle wheel devices (PW) were used for measuring flow velocities and surfboard meters (SB) for measuring flow depths in 8 locations along the crest and landward slope (see [39] for more details).To exclude the effects on the measurements of the landside road transition and roughness change, a geotextile cover was placed over the transition between the asphalt road and the remaining crest part (Figure 2, zone C).Additionally, the WOS was located directly over the road (Figure 2, zone B), to avoid the induced effects on the flow measurements from the riverside transition.Velocities and depths time series were measured at least two times for wave volumes (V) of 0.4, 0.6, 0.8, 1.0, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, and 5.5 m 3 /m.

Methodology
The complete methodology is divided in 4 steps as shown in the flow chart presented in Figure 3.The first step consists of constructing and validating two CFD models built for producing accurate and detailed time series of bottom shear stresses (τ) along the dike profile as a function of different wave volumes (V).The first model includes an asphalt road at the crest of the dike and corresponds to the profile of the Millingen dike initial state scan.This CFD model is referred to as "Road on crest dike" model (RCD).The second model represents the same dike but the road on top was replaced by a smooth grass cover.This CFD model is referred to as "Grass crest dike" model (GCD).
In the second step (Section 4.2) the bottom stress time series (τ) are integrated at each location for different ranges of critical shear stress values (τ c ).The obtained set of values represent the shear stress excess integral T, which are used to build two computationally inexpensive models (emulators) which imitate both CFD models.These emulators are capable of estimating T values as a function of combinations of V and τ c values.
In the third step (Section 4.3), the "grass cover quality" is determined for each location and represented by a C E function.These functions are calibrated based on the previously built emulators and the scanned profiles of the Millingen experiment part I.One of the main assumptions of the methodology is that the estimated curves are representative for both dike conditions.
Finally, in the fourth step (Section 4.4) the emulators of the RCD and GCD are used to predict the erosion depth at each location along the dike profile, using stochastic input for V and τ c and the calibrated C E functions from step 3.This yields a random set of erosion profiles from which the failure probabilities in each location for both the road and grass covered dike are derived for each storm condition (q m ).Note that the use of emulators allowed us to reduce the computational times so that crude Monte Carlo simulations could be used that captured the high level of detail of the hydrodynamics instead of using approximate reliability methods such as first order reliability method (FORM) or second order reliability method (SORM) [31].
In the second step (Section 4.2) the bottom stress time series ( ) are integrated at each location for different ranges of critical shear stress values ( ).The obtained set of values represent the shear stress excess integral , which are used to build two computationally inexpensive models (emulators) which imitate both CFD models.These emulators are capable of estimating values as a function of combinations of and values.In the third step (Section 4.3), the "grass cover quality" is determined for each location and represented by a function.These functions are calibrated based on the previously built emulators and the scanned profiles of the Millingen experiment part I.One of the main assumptions of the methodology is that the estimated curves are representative for both dike conditions.
Finally, in the fourth step (Section 4.4) the emulators of the RCD and GCD are used to predict the erosion depth at each location along the dike profile, using stochastic input for and and the calibrated functions from step 3.This yields a random set of erosion profiles from which the failure probabilities in each location for both the road and grass covered dike are derived for each storm condition ( ).Note that the use of emulators allowed us to reduce the computational times so that crude Monte Carlo simulations could be used that captured the high level of detail of the hydrodynamics instead of using approximate reliability methods such as first order reliability method (FORM) or second order reliability method (SORM) [31].

CFD Models
Detailed two dimensional (2D) CFD Reynolds averaged Navier-Stokes equations (RANS) with k-turbulence closure solution models are considered better approximations of flow models where transient bed shear stresses are required with respect to uniform flow-based methods.They simulate turbulent flows for which rapid dissipation of kinematic energy is converted into internal energy in the form of turbulence.The inclusion of turbulence derived effects in the bed shear stresses becomes more important when abrupt geometric bottom changes and local surface roughness variations are present (e.g., RCD bottom profile, see Figure 4).
The RCD and GCD models were built using the COMSOL Multiphysics ® software (COMSOL Inc., Stockholm, Sweden) for solving the 2D RANS k-equations for a two-phase field domain (air/water) based on a two equation (k-) turbulence closure model.This model accounts for turbulence effects by including two additional unknowns; the turbulent kinetic energy ( ) and the  The RCD and GCD models were built using the COMSOL Multiphysics ® software (COMSOL Inc., Stockholm, Sweden) for solving the 2D RANS k-ε equations for a two-phase field domain (air/water) based on a two equation (k-ε) turbulence closure model.This model accounts for turbulence effects by including two additional unknowns; the turbulent kinetic energy (k) and the rate of dissipation of kinetic energy (ε).Both are used for representing the transfer of momentum inside the flow due to turbulent eddies represented by the eddie viscosity.In addition, a phase field two equation transport approach is used for the interface tracking.For further information of the numerical solution see Ref. [40,41].Details of the construction of the grid are given in Ref. [35].
The first CFD model (RCD), originally built in Ref. [26,35], was based on the scanned profile obtained as a final condition after the 0.01 m 3 /s/m experiment (q10_sc2: where q implies the mean overtopping discharge in l/s/m and sc2 refers to scan 2, the first profile scan after the q10 event, see Table 1).This profile includes an asphalt road of 3.1 m width with adjacent eroded transition gaps on each side of approximately 0.5 m each (Figure 4).Additionally, 2.5 m of the landside slope (approx.1:3) is included in the model domain.For a detailed description of the RCD model setup, the accuracy of the mesh refinement and selected turbulence model we refer to Ref. [26].
A second model (GCD) was built to recreate a dike which could be tested with the same hydrodynamic conditions, but without including the asphalt road and transition gaps present in the original model.This model was built based on the original Millingen experiment dike profile (q10_sc2, Section 3.1), but the transition gaps were removed with a line smoothing procedure.The asphalt cover roughness was replaced by the roughness of grass (Figure 4).The rest of the profile remained unmodified in terms of the original elevation and slopes with respect to the RCD model.
The simulation principle consists of defining a water domain inside the WOS equivalent to the wave volume to be routed and run the model by letting the water to flow free from the domain during 15 s to ensure the whole volume has left the system (besides water retained in concave zones).
rate of dissipation of kinetic energy ( ).Both are used for representing the transfer of momentum inside the flow due to turbulent eddies represented by the eddie viscosity.In addition, a phase field two equation transport approach is used for the interface tracking.For further information of the numerical solution see Ref. [40,41].Details of the construction of the grid are given in Ref. [35].
The first CFD model (RCD), originally built in Ref. [26,35], was based on the scanned profile obtained as a final condition after the 0.01 m 3 /s/m experiment (q10_sc2: where q implies the mean overtopping discharge in l/s/m and sc2 refers to scan 2, the first profile scan after the q10 event, see Table 1).This profile includes an asphalt road of 3.1 m width with adjacent eroded transition gaps on each side of approximately 0.5 m each (Figure 4).Additionally, 2.5 m of the landside slope (approx.1:3) is included in the model domain.For a detailed description of the RCD model setup, the accuracy of the mesh refinement and selected turbulence model we refer to Ref. [26].
A second model (GCD) was built to recreate a dike which could be tested with the same hydrodynamic conditions, but without including the asphalt road and transition gaps present in the original model.This model was built based on the original Millingen experiment dike profile (q10_sc2, Section 3.1), but the transition gaps were removed with a line smoothing procedure.The asphalt cover roughness was replaced by the roughness of grass (Figure 4).The rest of the profile remained unmodified in terms of the original elevation and slopes with respect to the RCD model.
The simulation principle consists of defining a water domain inside the WOS equivalent to the wave volume to be routed and run the model by letting the water to flow free from the domain during 15 s to ensure the whole volume has left the system (besides water retained in concave zones).Both GCD and RCD models are used for producing the bottom shear stress time series ( ( )) along the bottom profile by releasing a certain volume contained in the WOS (see Figure 4).The bottom shear stress is calculated as the product of the flow density times the boundary shear velocity ( * ) squared.This shear velocity is calculated parallel to the bottom surface which allows to account for the fact that the erosion model was developed assuming a flat bottom.This model also allows one to estimate the tensile stresses in the normal direction of the bottom, but for the present method it is assumed that the erosion is a product of only the tangential shear stresses.The Manning's roughness coefficients used for both RCD and GCD models were obtained from the measured values reported in Ref. [39].The steel roughness was used as the calibration value for achieving the measured velocities in the outlet of the WOS.Further information about the calibration may be found in Ref.Both GCD and RCD models are used for producing the bottom shear stress time series (τ(t)) along the bottom profile by releasing a certain volume contained in the WOS (see Figure 4).The bottom shear stress is calculated as the product of the flow density times the boundary shear velocity (u * ) squared.This shear velocity is calculated parallel to the bottom surface which allows to account for the fact that the erosion model was developed assuming a flat bottom.This model also allows one to estimate the tensile stresses in the normal direction of the bottom, but for the present method it is assumed that the erosion is a product of only the tangential shear stresses.The Manning's roughness coefficients used for both RCD and GCD models were obtained from the measured values reported in Ref. [39].The steel roughness was used as the calibration value for achieving the measured velocities in the outlet of the WOS.Further information about the calibration may be found in Ref. [1].The transformation from Manning values (n) to equivalent sand roughness (K s ) values required as input by the software was done with Ref. [42]: The values used for the different boundary material roughness are presented in Table 2. 1 This value was not used as input in any of the two CFD models, because its value is similar to grass.
Simulations for 0.15, 0.4, 0.7, 1.0, 1.5, 2.5, 3.2, and 4.5 m 3 /m wave volumes (V) were performed in both RCD and GCD models.These values were chosen such that different orders of magnitude could be covered in the emulator training set while also including wave volumes which were measured in the Millingen experiment (0.15, 0.4, 1.0, and 2.5 m 3 /m) for the CFD models validation (presented in appendices A, B, and C).From each simulation, bottom shear stress time series were generated in 33 different locations along the profiles of both models (Figure 5).From now on these locations are referred to as STP1 until STP33.
J. Mar.Sci.Eng.2018, 6, x FOR PEER REVIEW 9 of 28 [1].The transformation from Manning values ( ) to equivalent sand roughness ( ) values required as input by the software was done with Ref. [42]: The values used for the different boundary material roughness are presented in Table 2. Simulations for 0.15, 0.4, 0.7, 1.0, 1.5, 2.5, 3.2, and 4.5 m 3 /m wave volumes ( ) were performed in both RCD and GCD models.These values were chosen such that different orders of magnitude could be covered in the emulator training set while also including wave volumes which were measured in the Millingen experiment (0.15, 0.4, 1.0, and 2.5 m 3 /m) for the CFD models validation (presented in appendices A, B, and C).From each simulation, bottom shear stress time series were generated in 33 different locations along the profiles of both models (Figure 5).From now on these locations are referred to as STP1 until STP33.These 33 study points (numbered dots in Figure 5) were selected for estimating the erosion depths at the locations where the bottom slope line segment changed significantly with respect to the previous line segment.All 33 STPs are located in the same horizontal position in both RCD and GCD models.The manual smoothing procedure of the RCD profile for generating the GCD profile consisted of placing straight lines between STPs 5 and 7, 9 to 20, and 23 to 26 (Figure 5).According to the field measurements [39,45], the average thickness of the grass cover was around 10 cm (dashed lower line in Figure 5).Points that are below this line in the RCD are located inside the clay zone (e.g., STPS 10 to 13 in the RCD, Figure 5).These 33 study points (numbered dots in Figure 5) were selected for estimating the erosion depths at the locations where the bottom slope line segment changed significantly with respect to the previous line segment.All 33 STPs are located in the same horizontal position in both RCD and GCD models.The manual smoothing procedure of the RCD profile for generating the GCD profile consisted of placing straight lines between STPs 5 and 7, 9 to 20, and 23 to 26 (Figure 5).According to the field measurements [39,45], the average thickness of the grass cover was around 10 cm (dashed lower line in Figure 5).Points that are below this line in the RCD are located inside the clay zone (e.g., STPS 10 to 13 in the RCD, Figure 5).
The RCD model and its validation are presented in Ref. [35] and Appendixs A-C of the present study.As the Millingen experiment included a road at the crest, no experimental measurements were available for validating the GCD model.Yet, maximum value data from two other WOS experiments performed in the Vechtdijk in Zwolle [46] and the Tholen dike near Nieuw-Strijen [39] were used for validating the GCD model for which good agreement was obtained (See Appendixs A and B).

Emulator Surfaces Construction
Bed shear stress time series produced for each STP location from each CFD model and each wave volume (0.15, 0.4, 0.7, 1.0, 1.5, 2.5, 3.2, and 4.5 m 3 /m) are used to generate the training data of the emulators.To do that, each of these time series was integrated in time for one given value of τ c in steps of 1 N/m 2 until 300 N/m 2 .The τ c range was defined based on the recommended values for grass and soil presented in Ref. [36].As a result, lists of training data sets of three columns (V, τ c , and T) for each STP location were obtained and used to train the emulators.
A data-based emulation approach was selected [27], as we are only interested in the final T values for estimating the erosion depth in each location.A data-based emulator reproduces input-output relations from the input-output datasets produced by the CFD models.A 3D linear interpolation surface was used to build the emulators [47].Each emulator on each location is defined as Ω R n and Ω G n according to the referred dike model.The Greek letter Ω denotes emulator, letters R or G denote either "RCD" or "GCD" and the sub index n denotes the STP location.These emulators were built based on the Matlab ® (MathWorks Inc., Natick, MA, U.S.) interpolation scientific package and were designed as 3D surfaces for estimating the T G n (V, τ c ) and T R n (V, τ c ) functions in each STP.The emulators on STP8 and STP28 (Ω R 8 , Ω G 8 , Ω R 28 , Ω G 28 ) are plotted as an example in Figure 6 to show what the surfaces look like.In total, 66 emulators were built to predict T as a function of V, τ c ; 33 for the RCD and 33 for the GCD.These emulators can calculate the set of T values of 1000 storms in less than 10 s on a standard personal computer.Note that the surface that corresponds to the emulator of a dike with a road on top (Ω R 28 ) presents a less smooth behaviour along the slope with respect to the grass emulation surface in the same location (Ω G 28 ), which demonstrates that there is indeed a significant effect derived from the presence of the road upstream.This also shows that the overtopping bed shear stress process is non-linear and may be under-or over-estimated when irregularities are not explicitly considered.
The RCD model and its validation are presented in Ref. [35] and appendices A, B, and C of the present study.As the Millingen experiment included a road at the crest, no experimental measurements were available for validating the GCD model.Yet, maximum value data from two other WOS experiments performed in the Vechtdijk in Zwolle [46] and the Tholen dike near Nieuw-Strijen [39] were used for validating the GCD model for which good agreement was obtained (See Appendices A and B).

Emulator Surfaces Construction
Bed shear stress time series produced for each STP location from each CFD model and each wave volume (0.15, 0.4, 0.7, 1.0, 1.5, 2.5, 3.2, and 4.5 m 3 /m) are used to generate the training data of the emulators.To do that, each of these time series was integrated in time for one given value of in steps of 1 N/m 2 until 300 N/m 2 .The range was defined based on the recommended values for grass and soil presented in Ref. [36].As a result, lists of training data sets of three columns ( , , and ) for each STP location were obtained and used to train the emulators.
A data-based emulation approach was selected [27], as we are only interested in the final values for estimating the erosion depth in each location.A data-based emulator reproduces inputoutput relations from the input-output datasets produced by the CFD models.A 3D linear interpolation surface was used to build the emulators [47].Each emulator on each location is defined as Ω and Ω according to the referred dike model.The Greek letter Ω denotes emulator, letters R or G denote either "RCD" or "GCD" and the sub index denotes the STP location.These emulators were built based on the Matlab ® (MathWorks Inc., Natick, MA, U.S.) interpolation scientific package and were designed as 3D surfaces for estimating the ( , ) and ( , ) functions in each STP.The emulators on STP8 and STP28 (Ω , Ω , Ω , Ω ) are plotted as an example in Figure 6 to show what the surfaces look like.In total, 66 emulators were built to predict as a function of , ; 33 for the RCD and 33 for the GCD.These emulators can calculate the set of values of 1000 storms in less than 10 s on a standard personal computer.Note that the surface that corresponds to the emulator of a dike with a road on top (Ω ) presents a less smooth behaviour along the slope with respect to the grass emulation surface in the same location (Ω ), which demonstrates that there is indeed a significant effect derived from the presence of the road upstream.This also shows that the overtopping bed shear stress process is non-linear and may be under-or over-estimated when irregularities are not explicitly considered.For estimating the scour depth per wave, the emulators become very powerful as they are equivalent to the solution of the excess shear stress integrals present in the right side of Equation (4) (Figure 1), as presented in Equations ( 7) and ( 8) for both GCD and RCD cases, respectively.For estimating the scour depth per wave, the emulators become very powerful as they are equivalent to the solution of the excess shear stress integrals present in the right side of Equation (4) (Figure 1), as presented in Equations ( 7) and ( 8) for both GCD and RCD cases, respectively.
These emulators are valid for estimating the T values of any given wave volume (V) value between 0.150 and 3.2 m 3 /m and for any given τ c value between 1 and 300 N/m 2 .The τ c value reflects the grass quality and hence when erosion occurs.At the same time, the resistance is also affected by the coefficient of erodibility C E (see Equation ( 4)) which represents the resistance of the cover to be eroded which is referred in this manuscript as the grass quality.

C E Functions
The C E coefficient determines the resistance of the grass to be eroded as a function of the τ c in the selected erosion model as shown in Equation ( 2) [19].In the present methodology, it is assumed that C E is a function of τ c instead of a constant value which ensures that the same grass quality is reflected in both the estimation of T and the erosion depth estimation so that Equation (2) still remains valid for any τ c .These C E functions were built using the pre-trained emulators, the WOS released volume list from the experiment Part I, and the q10_sc2 and q50_sc5 scanned profiles [39].These two scans represent the initial and final state of the 0.05 m 3 /s/m storm condition experiment.The functions were calculated for each STP as: This equation was obtained by substituting Equation ( 8) into Equation ( 4).This equation allowed us to generate C E curves for each SPT for an arbitrary τ c range between 1 and 300 N/m 2 in steps of 1 N/m 2 .In each step, the actual release list of the volumes corresponding with the 0.050 m 3 /s/m storm condition was used as input for the RCD emulators to obtain a cumulative value of T for evaluating Equation (9).The initial and final profiles of the 0.050 m 3 /s/m storm condition were chosen for building the C E functions, as for storm conditions with a lower q m no significant scouring was obtained during the experiment.The scour depths (ε) per location due to the q m = 0.050 m 3 /s/m experiment were obtained in each STP by subtracting both profile scans (Figure 7 lower plot).
Ω ( , ) ≈ ( ( , ) − ) These emulators are valid for estimating the values of any given wave volume ( ) value between 0.150 and 3.2 m 3 /m and for any given value between 1 and 300 N/m 2 .The value reflects the grass quality and hence when erosion occurs.At the same time, the resistance is also affected by the coefficient of erodibility (see Equation ( 4)) which represents the resistance of the cover to be eroded which is referred in this manuscript as the grass quality.

Functions
The coefficient determines the resistance of the grass to be eroded as a function of the in the selected erosion model as shown in Equation ( 2) [19].In the present methodology, it is assumed that is a function of instead of a constant value which ensures that the same grass quality is reflected in both the estimation of and the erosion depth estimation so that Equation (2) still remains valid for any .These functions were built using the pre-trained emulators, the WOS released volume list from the experiment Part I, and the q10_sc2 and q50_sc5 scanned profiles [39].These two scans represent the initial and final state of the 0.05 m 3 /s/m storm condition experiment.The functions were calculated for each STP as: This equation was obtained by substituting Equation ( 8) into Equation ( 4).This equation allowed us to generate curves for each SPT for an arbitrary τ range between 1 and 300 N/m 2 in steps of 1 N/m 2 .In each step, the actual release list of the volumes corresponding with the 0.050 m 3 /s/m storm condition was used as input for the RCD emulators to obtain a cumulative value of for evaluating Equation (9).The initial and final profiles of the 0.050 m 3 /s/m storm condition were chosen for building the functions, as for storm conditions with a lower no significant scouring was obtained during the experiment.The scour depths ( ) per location due to the = 0.050 m 3 /s/m experiment were obtained in each STP by subtracting both profile scans (Figure 7 lower plot).They are equivalent to the solution of the negative integral present in Equation ( 9) for each location.For STPs 18, 19, 27, 28, and 29 where accretion was observed, the results were discarded and replaced by the closest one where erosion was found as the erosion model assumes that all material is fully washed to the downstream part.
Different grass cover relative densities (ρ cover ) per unit of width are required for each C E curve depending on whether the final profile is located in the grass layer or inside the deeper clayey soil layer.As no field data was collected regarding this parameter, density values for each type of cover (grass or soil), reported in Ref. [48] were used as reference.They correspond to a saturated sample extracted from a Danish dike which was composed of 0.17 m of soil and 0.03 m of grass with a total average density of 1870 kg/m 3 .For this study, we assumed a reference density per unit of width value of soil (ρ cover soil ) of 2000 kg/m 3 .Based on these values, it is estimated that the saturated density of grass solely (ρ cover grass ) is 1100 kg/m 3 .The resultant C E curves obtained by the use of Equation ( 9) are presented in the results Section 5.2.It is acknowledged by the authors that these values may differ from the actual ones of Millingen.However, they are of good order of magnitude and even result in very similar C E values with respect to the ones reported in Ref. [19].

Probabilistic Safety Assessment
In general, failure mechanisms can be expressed as a limit state equation: where the variable S represents the load (solicitation) and R represents the resistance of the structure.The Z term represents the marginal resistance which defines the state of the system as safe when positive and in failure when equal to zero or negative.Wave overtopping failure is a major threat to flood defence safety, because when the grass cover is completely lost, a rapid scouring process of the soil core may initiate a dike breach.Hence for the present study, failure is defined as the complete loss of the available grass cover in any location (STP).The grass cover layer thickness (δ) is defined as the top 10 cm (Figure 5), based on Ref. [32,38].
The dike grass cover limit state function is obtained by rearranging Equation ( 4) and making it equal to the marginal strength term for each overtopped wave and then summing over all overtopping waves (Equation ( 11)).Here, δ represents the resistance term (R) and the second term describes the erosion that represents the solicitation (S).
The failure state is reached if the grass cover is completely removed (Z ≤ 0).The emulators were used to compute the excess shear stress per overtopping wave.All overtopping waves during a storm condition together yield the accumulated scour depth, which is evaluated for every STP along the dike profiles for the cases with (Z R ) and without (Z G ) a road on top.The δ values in each location and their corresponding emulators were replaced in Equation (11) such that the limit state equation of the grass cover for both cases with and without a road in each location could be evaluated.The resultant limit state equations per STP are defined as the accumulated scour depth per storm condition (Z R and Z G ) as shown in Equations ( 12) and ( 13): The available grass cover per location (δ R n and δ G n ) is defined as a deterministic variable and was calculated as the difference in elevation between the surface profile and the soil core line (dashed lower lines in Figure 5).The cover relative densities of grass (ρ cover ) were assumed to be constant all along the dike profile.
For every storm condition, possible storm realizations that consist of N wo sets of values of V, τ c were randomly sampled via crude Monte Carlo.These values are used to evaluate Equations ( 12) and ( 13) which result in one single value of Z R and Z G per realization for each STP.Statistical distributions of Z R and Z G per location were built from 1000 realizations of overtopping waves.For each tested storm condition (see Table 3, q m ), the random sampling was done following a pre-fitted Weibull distribution.The failure probabilities per storm condition are then calculated from the distributions as For each realization, a final eroded profile was also estimated by subtracting the accumulated scour depth (second term of Equations ( 12) and ( 13)) from the initial elevation profile (green lines of Figure 5).Note that each estimated value is conditioned to the event of the simultaneous occurrence of the boundary conditions of wind speed and water level.
For the present study, the same boundary conditions were used for each storm condition but the crest level was modified so that different q m values could be obtained.Individual wave volumes, V, were sampled from the two parameter (a and b) Weibull distribution [49].For the present study, it is intended to represent storm/wind surge overtopping events of duration D [h] composed of individual single wind generated waves V i that result in an overtopped volume distribution [50].The cumulative exceedance function of this distribution is expressed as: For the scale (a) and shape (b), the method presented in [50] improved the estimation of b for wave overtopping by fitting an empirical equation to experimental results (Equation ( 16)).
where q m represents the average overtopping discharge [m 3 /s/m] of the storm, % ov represents the overtopping percentage [-], R c represents the dike free crest height [m], T m represents the mean wave period [s], and H m0 represents the incident energy based significant wave height [m].Former Dutch guidelines [11] define the allowable q m between 0.0001 and 0.01 m 3 /s/m depending on the type of cover, the outer slope revetment and the importance of the structure in the system [11].This range plus four larger storm conditions were tested (Table 3).All the storm conditions tested in this study (Table 3) corresponded to a wave (H s = 1 m and T p = 4 s), wind speed, and water level event with exceedance probability P(H s , T p , q m ) = 1 × 10 −4 , which is the design event (in the former legislation) for this particular dike location.All estimated failure probabilities obtained from the emulators are conditioned to this event (P f R Z R ≤ 0 H s , T p , q m and ( P f G (Z G ≤ 0 H s , T p , q m ) ).For calculating the failure probability per location on each model we rely on the use of the conditional probability so that P f Z ≤ 0 ∩ H s , T p , q m = P f Z ≤ 0 H s , T p , q m × P(H s , T p , q m ).
For the τ c random variable, there is no available literature to our knowledge that studied or suggested the stochastic nature (distribution or moments) of this variable for grass.However, from the values reported in Ref. [39], an equivalence curve between critical erosion velocities U c and τ c was built based on the table presented in Ref. [36], as shown in Figure 8.This equivalence curve allowed the transformation to define the C E as a stochastic random variable.
For the random variable, there is no available literature to our knowledge that studied or suggested the stochastic nature (distribution or moments) of this variable for grass.However, from the values reported in Ref. [39], an equivalence curve between critical erosion velocities and was built based on the table presented in Ref. [36], as shown in Figure 8.This equivalence curve allowed the transformation to define the as a stochastic random variable.
Figure 8. Equivalence curve between critical velocity, , and critical shear stress, , built from values reported in [14].This curve is used to include as a stochastic random variable by transforming the stochastic from Table 4 to values, which is required as input for using Equation (9).
The stochastic distributions for the of grass and soil, were taken from Ref. [21].This work concluded that can be represented by a log-normal distribution for the grass cover and by a generalized extreme value (GEV) distribution for bare soil locations.The parameters for the GEV distribution could not be estimated from the experimental data collected during the Millingen experiment.Hence, a log-normal distribution was also assumed for the bare soil spots, as presented in Table 4.  [36], it was possible to define mean values of per grass and soil qualities.The coefficients of variation (CoV) were calculated from the statistical distributions presented by Trung [21].The grass quality correction factor (QCF) is obtained by estimating the value which shifts the "average" grass quality function towards the recommended values in Ref. [19].[14].This curve is used to include C E as a stochastic random variable by transforming the stochastic U c from Table 4 to τ c values, which is required as input for C E using Equation (9).
The stochastic distributions for the U c of grass and soil, were taken from Ref. [21].This work concluded that U c can be represented by a log-normal distribution for the grass cover and by a generalized extreme value (GEV) distribution for bare soil locations.The parameters for the GEV distribution could not be estimated from the experimental data collected during the Millingen experiment.Hence, a log-normal distribution was also assumed for the bare soil spots, as presented in Table 4. Based on the indicative average values of τ c presented by Ref. [36], it was possible to define mean values of U c per grass and soil qualities.The coefficients of variation (CoV) were calculated from the U c statistical distributions presented by Trung [21].The grass quality correction factor (QCF) is obtained by estimating the value which shifts the "average" C E grass quality function towards the recommended values in Ref. [19].

CFD Calibration and Validation
Calibration and validation of the RCD model were carried out as in Ref. [26], using the observations of the Millingen experiment part II.In Ref. [26], it is concluded that although the depth simulations show some deviations, the simulated velocities along the crest and slope are within 20% of the observations.A detailed analysis of the effects and estimated errors derived from the mesh refinement, the CFD RANS parametrization, and the temporal effects in the shear stress time series error can be found in Ref. [26].The flow depths and velocities for the GCD case are in good agreement with the ones fitted for the Vechtdijk experiment (Appendixs A and B).Also, in this case, the simulated velocities are more accurate than the depth simulations.In the modelling approach, erosion is computed based on excess shear stress.Appendix C shows that for most wave volumes, both the GCD and RCD models simulate the duration of the overtopping waves quite well.Especially, the magnitude and timing of the velocity and depth at the wave fronts are simulated well.

Effects of Turbulence on the Excess Shear Stress T
Both the RCD and GCD model results show highly turbulent flow over the dike profile with Reynolds between 10,000 and 200,000 for the 0.15 and 3.2 m 3 /m volumes, respectively, with supercritical flow (Froude between 10 and 2 for 0.15 and 3.2 m 3 /m volumes, respectively).Low Froude numbers are observed for the RCD model in STP 14 Fr ≈ 1 at the adverse bottom slope.For abrupt changes in the bottom geometry, the use of a single depth-averaged turbulence factor along the entire profile may either underestimate or overestimate the locally exerted bottom shear stress as the turbulence intensity changes significantly along the dike profile.This can be observed in Figure 9A, where the average kinematic energy in time ( k(t)) for waves of 0.4 m 3 /m and 3.2 m 3 /m in both GCD and RCD models is shown.These plots show that the highest values of turbulence occur at the locations where abrupt bottom changes (irregularities) are observed.

CFD Calibration and Validation
Calibration and validation of the RCD model were carried out as in Ref. [26], using the observations of the Millingen experiment part II.In Ref. [26], it is concluded that although the depth simulations show some deviations, the simulated velocities along the crest and slope are within 20% of the observations.A detailed analysis of the effects and estimated errors derived from the mesh refinement, the CFD RANS parametrization, and the temporal effects in the shear stress time series error can be found in Ref. [26].The flow depths and velocities for the GCD case are in good agreement with the ones fitted for the Vechtdijk experiment (Appendices A and B).Also, in this case, the simulated velocities are more accurate than the depth simulations.In the modelling approach, erosion is computed based on excess shear stress.Appendix C shows that for most wave volumes, both the GCD and RCD models simulate the duration of the overtopping waves quite well.Especially, the magnitude and timing of the velocity and depth at the wave fronts are simulated well.

Effects of Turbulence on the Excess Shear Stress
Both the RCD and GCD model results show highly turbulent flow over the dike profile with Reynolds values between 10,000 and 200,000 for the 0.15 and 3.2 m 3 /m volumes, respectively, with supercritical flow (Froude between 10 and 2 for 0.15 and 3.2 m 3 /m volumes, respectively).Low Froude numbers are observed for the RCD model in STP 14 Fr ≈ 1 at the adverse bottom slope.For abrupt changes in the bottom geometry, the use of a single depth-averaged turbulence factor along the entire profile may either underestimate or overestimate the locally exerted bottom shear stress as the turbulence intensity changes significantly along the dike profile.This can be observed in Figure 9A, where the average kinematic energy in time ( ( )) for waves of 0.4 m 3 /m and 3.2 m 3 /m in both GCD and RCD models is shown.These plots show that the highest values of turbulence occur at the locations where abrupt bottom changes (irregularities) are observed.Nonetheless, the amount of erosion does not solely depend on the amount of turbulence but also on the momentum of the fluid.This implies that high erosion potential may also occur at low turbulence locations due to the acceleration of the flow.Hence, the average of the product between Nonetheless, the amount of erosion does not solely depend on the amount of turbulence but also on the momentum of the fluid.This implies that high erosion potential may also occur at low turbulence locations due to the acceleration of the flow.Hence, the average of the product between mean velocity (U(t)) and kinematic energy (k(t)) averaged in time allows one to determine the highest potential erosion locations, as shown in Figure 9B.Note that these plots do not directly imply which locations are more prone to erosion as the cover resistance is also a determining factor.
The bottom shear stress τ(t) intrinsically contains the effects of turbulence.Hence, they are influencing the probabilistic estimation as they also influence T. To observe this effect, the release volume list used as input for the 0.050 m 3 /s/m test is used to calculate a T value per volume while assuming τ c = 0. From the produced set of T values, the average shear stress excess (T) value and 95% confidence bounds per location for each of the 33 locations in both cases were calculated and plotted in Figure 10.
mean velocity ( ( )) and kinematic energy ( ( )) averaged in time allows one to determine the highest potential erosion locations, as shown in Figure 9B.Note that these plots do not directly imply which locations are more prone to erosion as the cover resistance is also a determining factor.
The bottom shear stress ( ) intrinsically contains the effects of turbulence.Hence, they are influencing the probabilistic estimation as they also influence .To observe this effect, the release volume list used as input for the 0.050 m 3 /s/m test is used to calculate a value per volume while assuming = 0. From the produced set of values, the average shear stress excess ( ) value and 95% confidence bounds per location for each of the 33 locations in both cases were calculated and plotted in Figure 10.The ratio between the values of the RCD and GCD cases (Figure 10, bottom) shows that is around 2 to 4 times larger in the road transition zones, which correspond to the high turbulent kinetic energy regions on Figure 9. Additionally, the values in the upper part of the landward slope (between STPs 26 and 33) were larger in the RCD with respect to the GCD.This behavior was surprising as the landward transition gap was expected to dissipate energy.However, the turbulence generated in the transition propagates downstream (see ( ) ( ) in Figure 9) and probably induces an increase in the potential erosion of the landward upper slope.
From these results, it is also possible to identify the effects of surface roughness and bottom surface irregularities separately.On one hand, due to form roughness, the values increase significantly after STP 9 where evident abrupt bottom changes are present in the RCD case.On the other hand, the in the centerline location (STP8) of Figure 10 is lower in the RCD case, which can be attributed to the smoother surface.This can be corroborated from Figure 9 at the same location where the bottom slope is almost the same in both models.At this location, the turbulent kinetic energy is also very similar and yet the erosive potential is much greater in the crest for the GCD.
Except for STP8, all locations present equal or larger values in the RCD with respect to the GCD.From a probabilistic point of view this means that on average, the probability of scouring at these locations will be higher if grass quality and cover thickness remain constant along both profiles.
Interestingly, the 95% confidence bounds show a wider spreading with respect to the mean values in all locations with high erosive potential.This is concluded by correlating the spreading of the confidence bounds per location of Figure 10 with respect to the locations of highest values of ( ) ( ) in Figure 9B.Similarly, it can be observed that at locations where bottom slope changes are The ratio between the T values of the RCD and GCD cases (Figure 10, bottom) shows that T is around 2 to 4 times larger in the road transition zones, which correspond to the high turbulent kinetic energy regions on Figure 9. Additionally, the T values in the upper part of the landward slope (between STPs 26 and 33) were larger in the RCD with respect to the GCD.This behavior was surprising as the landward transition gap was expected to dissipate energy.However, the turbulence generated in the transition propagates downstream (see U(t)k(t) in Figure 9) and probably induces an increase in the potential erosion of the landward upper slope.
From these results, it is also possible to identify the effects of surface roughness and bottom surface irregularities separately.On one hand, due to form roughness, the T values increase significantly after STP 9 where evident abrupt bottom changes are present in the RCD case.On the other hand, the T in the centerline location (STP8) of Figure 10 is lower in the RCD case, which can be attributed to the smoother surface.This can be corroborated from Figure 9 at the same location where the bottom slope is almost the same in both models.At this location, the turbulent kinetic energy is also very similar and yet the erosive potential is much greater in the crest for the GCD.
Except for STP8, all locations present equal or larger T values in the RCD with respect to the GCD.From a probabilistic point of view this means that on average, the probability of scouring at these locations will be higher if grass quality and cover thickness remain constant along both profiles.
Interestingly, the 95% confidence bounds show a wider spreading with respect to the mean values in all locations with high erosive potential.This is concluded by correlating the spreading of the confidence bounds per location of Figure 10 with respect to the locations of highest values of U(t)k(t) in Figure 9B.Similarly, it can be observed that at locations where bottom slope changes are minimal (e.g., STP8), less spreading of the confidence bounds is observed which indicates that there is less turbulence.

CE Curves Calculated from Millingen Measurements
The C E curves obtained from the experiment data are presented in Figure 11.The results are presented in Figure 5 and correspond to an estimation of the corresponding curve of grass quality of each of the places where erosion was found along the dike cross sectional profile during the Millingen experiment.This figure only includes the resulting CE curves for STPs located after the asphalt cover (STPs > 9), which were influenced by the presence of the road.The right plot of Figure 11 shows the CE curves obtained for SPTs inside the gap where no grass cover was present and consequently correspond to soil erodability values.
minimal (e.g., STP8), less spreading of the confidence bounds is observed which indicates that there is less turbulence.

CE Curves Calculated from Millingen Measurements
The curves obtained from the experiment data are presented in Figure 11.The results are presented in Figure 5 and correspond to an estimation of the corresponding curve of grass quality of each of the places where erosion was found along the dike cross sectional profile during the Millingen experiment.This figure only includes the resulting CE curves for STPs located after the asphalt cover (STPs > 9), which were influenced by the presence of the road.The right plot of Figure 11 shows the CE curves obtained for SPTs inside the gap where no grass cover was present and consequently correspond to soil erodability values.curves for each STP for grass (left) and bare soil (right).Horizontal lines represent the grass quality reference intervals for very poor, poor, average, and good grass cover quality reported in [19].
The curves from the STPs located closer to the road (STPs 14 to 17 plotted as dense dashed lines in Figure 11) cross the grass quality reference intervals suggested by [19].Thus, we classified them between the very poor and poor grass quality.For STPs located over the vertex plotted as continuous thicker lines (STPs 20 and 26) coincide with average and good grass quality reference intervals.For STPs 25 and 26 which are located in the dent after the vertex, poor grass quality curves were found.The presence of such a dent may be a good indication of deterioration due to traffic which explains such low quality just after spots of good quality.The rest of the profile (mostly along the inner slope in STPs 30 to 33) present good grass quality.These results are in good agreement with the field findings for grass quality reported by [39].For the STPs located initially in the soil zone, good clay quality values were obtained (right plot of Figure 11).For the GCD case, no information about grass quality was available.Therefore, we used the average curve, corresponding to average grass quality.It is acknowledged that a dike without a road may have better quality in reality than what we assumed as we used the transition zone data in the averaging procedure for obtaining the "average curve".

Scour Depth Profiles
The results of the simulations for the storm conditions of 0.0001, 0.001, and 0.010 m 3 /s/m showed almost no scouring with the 1000 samples.Therefore, only the scour profiles for the largest (0.02-0.1 m 3 /s/m) were analyzed.Each profile was produced by plotting the mean value of the final scour depths for each STP after the 1000 realizations generated for each (Table 3).These values are equivalent to the mean value of the probability density function obtained for the second The C E curves from the STPs located closer to the road (STPs 14 to 17 plotted as dense dashed lines in Figure 11) cross the grass quality reference intervals suggested by [19].Thus, we classified them between the very poor and poor grass quality.For STPs located over the vertex plotted as continuous thicker lines (STPs 20 and 26) coincide with average and good grass quality reference intervals.For STPs 25 and 26 which are located in the dent after the vertex, poor grass quality C E curves were found.The presence of such a dent may be a good indication of deterioration due to traffic which explains such low quality just after spots of good quality.The rest of the profile (mostly along the inner slope in STPs 30 to 33) present good grass quality.These results are in good agreement with the field findings for grass quality reported by [39].For the STPs located initially in the soil zone, good clay quality values were obtained (right plot of Figure 11).For the GCD case, no information about grass quality was available.Therefore, we used the average C E curve, corresponding to average grass quality.It is acknowledged that a dike without a road may have better quality in reality than what we assumed as we used the transition zone data in the averaging procedure for obtaining the C E "average curve".

Scour Depth Profiles
The results of the simulations for the storm conditions of q m 0.0001, 0.001, and 0.010 m 3 /s/m showed almost no scouring with the 1000 samples.Therefore, only the scour profiles for the largest q m (0.02-0.1 m 3 /s/m) were analyzed.Each profile was produced by plotting the mean value of the final scour depths for each STP after the 1000 realizations generated for each q m (Table 3).These values are equivalent to the mean value of the probability density function obtained for the second term of Equations ( 12) and ( 13) for the 1000 sampled realizations in each STP.The profiles were only analyzed for the STP locations after the asphalt cover.
For the first part of the analysis, the simulations were done assuming a uniform average grass quality along the whole profile.The obtained scoured profiles (see Figure 12) show that the large failure spot in the lower part of the RCD may be explained from the additional energy available due to the presence of a smoother upstream surface (asphalt cover) and the momentum gain when flowing downslope.Moreover, an eventual increase of the turbulent flow due to the slope bottom irregularity around x = 3.8 m near STP 28 may also contribute.Another significant scour zone is observed on the transition between bare soil (STP12) and the dike crest vertex (STP23).These deep erosion spots coincide with the locations where the initial profile has an adverse slope.These zones are identified as the initiators of turbulence in the flow as it can be observed in Figure 9.
J. Mar.Sci.Eng.2018, 6, x FOR PEER REVIEW 18 of 28 term of Equations ( 12) and ( 13) for the 1000 sampled realizations in each STP.The profiles were only analyzed for the STP locations after the asphalt cover.
For the first part of the analysis, the simulations were done assuming a uniform average grass quality along the whole profile.The obtained scoured profiles (see Figure 12) show that the large failure spot in the lower part of the RCD may be explained from the additional energy available due to the presence of a smoother upstream surface (asphalt cover) and the momentum gain when flowing downslope.Moreover, an eventual increase of the turbulent flow due to the slope bottom irregularity around x = 3.8 m near STP 28 may also contribute.Another significant scour zone is observed on the transition between bare soil (STP12) and the dike crest vertex (STP23).These deep erosion spots coincide with the locations where the initial profile has an adverse slope.These zones are identified as the initiators of turbulence in the flow as it can be observed in Figure 9. From these results, the most interesting finding is that STP23 represents a weak spot in the profile despite the fact that the available grass cover is even thicker than most of the other locations.This can be explained based on Figures 9 and 10, where STP23 corresponds to one of the most turbulent locations.Again, this supports the importance of including a highly detailed geometry and its derived turbulence effects in modelling of dike cover erosion.
For a more realistic representation of the grass cover erosion, the grass qualities of each point were replaced by one of the three estimated qualities presented in Table 4, that seemed closer to the local condition per STP found in Section 5.2.The resultant scoured profiles (Figure 13) become safer in the lower part of the slope (STP 31) for both RCD and GCD as better grass quality than average was obtained in that zone.
However, two new close to failure locations were identified in STPs 25 and 26 which were not observed when assuming all grass quality as average.For the 0.1 m 3 /s/m RCD storm profile, it can be observed how STP 25 is very close to failure as the scoured depth is very close to the available grass cover thickness (lower left plot in Figure 13).This location corresponds to the most prone to fail spot in the GCD given the poor grass quality present in that particular location.From these results, the most interesting finding is that STP23 represents a weak spot in the profile despite the fact that the available grass cover is even thicker than most of the other locations.This can be explained based on Figures 9 and 10, where STP23 corresponds to one of the most turbulent locations.Again, this supports the importance of including a highly detailed geometry and its derived turbulence effects in modelling of dike cover erosion.
For a more realistic representation of the grass cover erosion, the grass qualities of each point were replaced by one of the three estimated qualities presented in Table 4, that seemed closer to the local condition per STP found in Section 5.2.The resultant scoured profiles (Figure 13) become safer in the lower part of the slope (STP 31) for both RCD and GCD as better grass quality than average was obtained in that zone.
However, two new close to failure locations were identified in STPs 25 and 26 which were not observed when assuming all grass quality as average.For the 0.1 m 3 /s/m RCD storm profile, it can be observed how STP 25 is very close to failure as the scoured depth is very close to the available grass cover thickness (lower left plot in Figure 13).This location corresponds to the most prone to fail spot in the GCD given the poor grass quality present in that particular location.

Probability of Failure
The failure probabilities for both dike cases are presented in Figure 14.For this study, any failure probability larger than 9.9 × 10 −5 is assumed as an already failed spot.This value was assumed since STPs 10, 11, 12, and 13 (located in bare soil) presented values above this threshold and were already in failure state (no grass cover) before routing each different storm condition and where not taken into account for the dike assessment.For the present study, it is assumed that the failure probability of the configuration to be studied is equal to the largest value found along the profile.This allows us not only to define the failure probability but also the spatial location of the failure which can later be related to the analysis of the origin of the failure.

Probability of Failure
The failure probabilities for both dike cases are presented in Figure 14.For this study, any failure probability larger than 9.9 × 10 −5 is assumed as an already failed spot.This value was assumed since STPs 10, 11, 12, and 13 (located in bare soil) presented values above this threshold and were already in failure state (no grass cover) before routing each different storm condition and where not taken into account for the dike assessment.

Probability of Failure
The failure probabilities for both dike cases are presented in Figure 14.For this study, any failure probability larger than 9.9 × 10 −5 is assumed as an already failed spot.This value was assumed since STPs 10, 11, 12, and 13 (located in bare soil) presented values above this threshold and were already in failure state (no grass cover) before routing each different storm condition and where not taken into account for the dike assessment.For the present study, it is assumed that the failure probability of the configuration to be studied is equal to the largest value found along the profile.This allows us not only to define the failure probability but also the spatial location of the failure which can later be related to the analysis of the origin of the failure.For the present study, it is assumed that the failure probability of the configuration to be studied is equal to the largest value found along the profile.This allows us not only to define the failure probability but also the spatial location of the failure which can later be related to the analysis of the origin of the failure.
Figure 14 shows that the dike without a road presents significant failure probabilities in the slope on STP31 with Pfs of 5 × 10 −5 and 8.2 × 10 −5 for the q m 0.075 and 0.1 m 3 /s/m.The results for this last configuration show that the probabilities of failure along the landward slope are almost identical for both RCD and GCD cases for the largest tested q m .This means that when using the failure criteria of assigning the maximum failure probability along the profile for a single storm it could be concluded that both dikes with and without road have almost the same failure probability.Yet, for lower q m (<0.075 m 3 /s/m), the probabilities of failure of the landward slope are always higher in the dike with a road.It is also observed on STP25, that a sudden potential failure spot is created for a q m larger than 0.050 m 3 /s/m.It is described as sudden because the failure probability was less than 1 ×10 −4 for lower storms with poor grass quality.This can be explained due to the highly turbulent flows expected in this location (see Figure 9).Also, the assessment was performed with the spatially varying grass quality based on observed grass quality in the Millingen experiment to represent a more realistic estimation.Here, a sudden failure spot is observed in STP25 for a 0.1 m 3 /s/m storm (see Figure 15).These results show that for the dike with a road, the maximum probability of failure is equal to 5 × 10 −5 for the largest storm in STP25.Based on all previous analyses, STP25 presented a high T value combined with poor grass quality and medium to high grass cover thickness (δ = 7.8 cms).This means that for this particular place, turbulence plays a more important role than cover thickness as STP25 represents the highest failure potential spot on the dike while also presenting the highest value of average kinetic energy (see Figure 9).From this last figure, it is expected that this dike will present its first failures occurring nearby the already scoured zones (STP14 to STP18) with Pf range between (1 × 10 −5 and 1 × 10 −4 ) for the case where the road is present.For the no road configuration, it is expected to have a Pf lower than 1 × 10 −6 .
J. Mar.Sci.Eng.2018, 6, x FOR PEER REVIEW 20 of 28 Figure 14 shows that the dike without a road presents significant failure probabilities in the slope on STP31 with Pfs of 5 × 10 −5 and 8.2 × 10 −5 for the 0.075 and 0.1 m 3 /s/m.The results for this last configuration show that the probabilities of failure along the landward slope are almost identical for both RCD and GCD cases for the largest tested .This means that when using the failure criteria of assigning the maximum failure probability along the profile for a single storm it could be concluded that both dikes with and without road have almost the same failure probability.Yet, for lower (<0.075m 3 /s/m), the probabilities of failure of the landward slope are always higher in the dike with a road.It is also observed on STP25, that a sudden potential failure spot is created for a larger than 0.050 m 3 /s/m.It is described as sudden because the failure probability was less than 1 ×10 −4 for lower storms with poor grass quality.This can be explained due to the highly turbulent flows expected in this location (see Figure 9).Also, the assessment was performed with the spatially varying grass quality based on observed grass quality in the Millingen experiment to represent a more realistic estimation.Here, a sudden failure spot is observed in STP25 for a 0.1 m 3 /s/m storm (see Figure 15).These results show that for the dike with a road, the maximum probability of failure is equal to 5 × 10 −5 for the largest storm in STP25.Based on all previous analyses, STP25 presented a high value combined with poor grass quality and medium to high grass cover thickness ( = 7.8 cms).This means that for this particular place, turbulence plays a more important role than cover thickness as STP25 represents the highest failure potential spot on the dike while also presenting the highest value of average kinetic energy (see Figure 9).From this last figure, it is expected that this dike will present its first failures occurring nearby the already scoured zones (STP14 to STP18) with Pf range between (1 × 10 −5 and 1 × 10 −4 ) for the case where the road is present.For the no road configuration, it is expected to have a Pf lower than 1 × 10 −6 .

Uncertainties in the Modelling Process
In this paper, a methodology is presented to predict the spatially varying probability of failure along the grass cover using emulation of a highly detailed CFD model.The Pf values obtained for each location are conditional to the upstream and the downstream locations.This means that the turbulence effects are captured in the emulators as they are trained based on the time series from the CFD model which includes them by solving for the kinematic energy variable in time and space.One of the main assumptions in this study is that the CFD models can be validated by checking the order of magnitude and trends of water depths and maximum velocities in two points along the profile (See appendix A, B, and C).However, the main output of the CFD models used in this study is the

Uncertainties in the Modelling Process
In this paper, a methodology is presented to predict the spatially varying probability of failure along the grass cover using emulation of a highly detailed CFD model.The P f values obtained for each location are conditional to the upstream and the downstream locations.This means that the turbulence effects are captured in the emulators as they are trained based on the time series from the CFD model which includes them by solving for the kinematic energy variable in time and space.One of the main assumptions in this study is that the CFD models can be validated by checking the order of magnitude and trends of water depths and maximum velocities in two points along the profile (See Appendixs A-C).However, the main output of the CFD models used in this study is the produced bed shear stress which cannot be validated at this stage, because no data was collected regarding this quantity during the experiments.Yet, we are confident that the model is giving consistent results in terms of order of magnitude.The uncertainties of the CFD model and assumptions in the erosion model determine the accuracy of the emulators and therefore the probabilities.Therefore, in order to increase the accuracy of the final failure probabilities, the first step is to improve and better validate the CFD models.In general, the presented methodology can be applied to other CFD models as well.
Additionally, the scouring process is highly time-dependent as the dike profile changes in time (not only due to erosion, but also due to accretion).Consequently, hydrodynamics change in time as well.In that sense, the emulators built for this study become more uncertain for larger storms as the chances of significant change in the dike profile become larger.This means that the updating of the bottom profile should be done after each routed wave and not after the complete storm.However, for the present methodology, the resultant profiles are a function of the emulators, which use a set of 1000 realizations for each storm, generated independently from the set of volumes used in the Millingen experiment and yet they result in very similar scour patterns compared to the one observed during the Millingen experiment.This indicates that even though the CFD model and its emulators were not perfect, and the profiles were not updated after each wave overtopping, the bed shear stress stochastic nature is fairly well represented by the emulators, plus it allows one to make scenario assessments without rerunning the original models.
Regarding the erosion model adopted for this study, it is also acknowledged that it fully attributes the complete erosion process to the bottom shear stress, whereas, in reality, normal stresses also contribute.This choice was made to actually reduce the uncertainty as using more complex erosion models requires a larger amount of data with their associated uncertainties potentially making the model even more uncertain.Still, the results obtained for the erodibility functions (C E ) also show good agreement with the suggested reference values deducted from field measurements and literature.This may be a good indicator that the simpler adopted erosion model is sufficiently robust for probabilistic modelling.

Sensitivity to the Grass Quality Estimation
Grass quality is highly variable in space and difficult to determine [36] in field experiments.Therefore, we inferred the grass quality for the test location from the observed erosion depths and emulators and we subsequently compared different configurations of grass quality.The inferred grass quality was better along the slope than on the crest where thickness was smaller, but the results showed that grass quality is not necessarily correlated to the available grass cover thickness which is another factor that determines the resistance to erosion.An example of this can be found for STP21 and STP26 where ε is larger than 0.10 m for both locations and still the estimated grass qualities were good and poor.These examples show that grass quality also has a significant effect on the occurrence of reduced failure probabilities.For example, the failure probability in the RCD with more realistic grass quality distribution was 5 × 10 −5 while presenting a poor local grass quality in STP 25.The failure probability in this same location for the same dike case is close to zero by just changing the grass quality to average (see Figure 14).This shows that the safety assessment cannot be done solely from the hydrodynamic point of view as the erosion pattern differs from the expected one using highest T values only.Poor grass quality spots have a high probability of failure despite the magnitude of q m .Additionally, the RCD more realistic model also presents a failure probability of 6 × 10 −5 in the landward vertex (STP 25, Figure 15) during the most severe storm, whereas the realistic GCD has almost a negligible probability of failure in this same location.The variability on the grass cover may be attributed to an external factor such as traffic which, for the case of a dike with a road on top, is inevitable.
Accurate observations of grass quality were therefore essential in our study case.However, we used a relatively simple approach to infer the grass quality from the erosion scanned profiles.For further applications of the methodology, detailed grass quality data is required, which can be derived from dike inspections or, possibly, remote sensing methods.
It is acknowledged that grass quality will also vary in the longitudinal direction.The present study assumes that the dominant direction of the generated shear stresses occurs in the cross-sectional direction (from crest to landside) and, therefore, the most probable way to determine failure is in that direction.However, with the constant modification of the cross section during the erosion process it can happen that grass quality and shear stresses importance vary their dominant direction.These effects cannot be considered with the present method and may be worth further studying.

Applications of the Probabilistic Method
The results show that the probability of failure if a road is present is higher along the simulated dike section, with respect to the case without a road, for all storms and grass qualities tested.However, for the GCD case with variable grass quality, all estimated failure probabilities (Figure 15) were less than 1 × 10 −6 .This result does not imply that the dike is 100% safe at these spots, because the small failure probabilities are predicted less accurately.The reason for the latter is that the 1000 samples of the Monte Carlo assessment were probably not sufficient to explore the tail of the distribution in depth as their associated P f can be smaller than the ones obtained for a Monte Carlo standard error >3%.The methodology presented in this paper enables the use of such a high number of computations within a feasible time period.However, a more extensive probabilistic analysis, including the uncertainties in the modelling process and increasing the number of samples (lower standard error in the Monte Carlo) and using smart sampling techniques to cover the tail will greatly increase the accuracy of the predicted probability of failure.
The method can be extended by including the wave parameters and water levels during the random sampling by making them stochastic variables fitted to data series of storm wave conditions and water levels [51].This allows one to extend the method to a fully probabilistic analysis by including the wind speeds water levels uncertainties explicitly.The resulting overtopping volumes are used as input for the emulators presented in this study, allowing one to derive fully probabilistic failure values along the dike profiles.Nevertheless, it remains a challenge to fully analyze the implications of such an extensive analysis and this is recommended for further study.
One of the main assumptions in the methodology was the definition of failure of the dike cover if the erosion depth exceeded the grass cover thickness of 10 cm.This assumption is highly arbitrary and debatable.However, no definition of failure for wave overtopping is available, other than a dike breach, which is a highly complex process caused by many interacting failure mechanisms and therefore too difficult to model.The adopted 10 cm threshold has obviously large consequences for the resulting probabilities.In order to achieve realistic values for the failure probability due to wave overtopping, we recommend more research to achieve a more advanced definition of the threshold of failure.
To summarize, the precise quantification of the low failure probabilities is more uncertain, but it is expected that they will still be in the obtained order of magnitude, which is sufficient to draw a good conclusion from the study in terms of trends and locations.The spatially distributed failure probabilities become useful, not only to determine the most prone to failure location, but also to stress their analysis during the yearly visual inspections for Dutch dikes.The quantification of the local probabilities of failure is input for any probabilistic dike safety assessment or probabilistic design.These results become even more relevant in cases of dike concepts such as the unbreachable dike and the multi-functional flood defence [51] which are designed to cope with climate change [5] and to withstand large amounts of overflow and storms similar or larger than the largest ones tested in this study.

Conclusions and Recommendations
This study shows a methodology to assess the spatially distributed probability of failure of a dike, with and without a road on the crest, due to wave overtopping erosion.In our case, the presence of a road on the dike increased the erosion of the grass cover and increased the probability of failure locally and for the entire dike profile.Probabilities were derived using a novel emulation method.Based on the obtained results of the study, specific conclusions for the two cases can be summarized as: 1.
The dike with a road showed higher probabilities (5 × 10 −5 > P f > 1 × 10 −4 ) of failure with respect to a dike without a road (P f <1 × 10 −6 ) if realistic grass quality distribution was assumed.
The probability of failure in the zones where higher initial deterioration was observed (e.g., locations immediately next to the road and over the vertex) is significantly higher with respect to the remaining part of the dike.Yet, these locations also correspond to spots where higher turbulence was observed with respect to the remaining part of the locations along the dike profile.However, if good quality grass was present all along the dike, both dike cases will present very low failure probabilities along the slope and the vertex (P f < 5 × 10 −5 ).

2.
Local erosion depth is highly influenced by the momentum of the water volume, which is reduced (energy dissipation) by a rough surface such as grass.If a road is present, both the smoother asphalt surface and the resultant bottom irregularities (possibly derived by traffic and material change) increase the scour potential for failure along the crest and the slope with respect to the dike without a road as they influence the generated and dissipated turbulence (kinetic energy) and its variability (turbulence intensities).
Finally, the coupled emulator-erosion model was able to yield realistic probabilities and trends, given all uncertainties in the modelling process.Yet, further improvements are required, especially for the CFD model validation, grass quality estimation, and erosion modelling.The present method can be extended by including model uncertainty coefficients as stochastic variables, to account for the errors originated from the imperfection of both the erosion and the hydrodynamic models.This is especially important to improve the erosion estimation of wave volumes located in the tail of the distribution.
The velocity measurements from the Millingen experiment and the results from the RCD and GCD models are in fairly good agreement with each other for waves with volumes which are less than 3.5 m 3 /m.The flow depths for both GCD and RCD models differ significantly from the recorded values during the experiment.Again, these differences could be attributed to the presence of the geotextile.The accuracy of the measuring devices should also be included in this analysis but unfortunately, no data was available in the literature.Yet, it is concluded that maximum velocities and maximum flow depths of the GCD model are in good agreement with the measurements taken at the Vechtdijk presented in [15].Additionally, the expected behavior presented by the GCD model (i.e., lower velocities and larger flow depths) with respect to the RCD model gives sufficient confidence for using the model in the present study.

Appendix C. CFD Validation for Overtopping Duration
The simulated overtopping durations presented in Figure A3 are in good agreement with the ones measured by the surf boards on the crest and slope but not with the ones computed from the paddle wheel velocity meters.Apparently, there is a discrepancy between both measuring devices.No thorough analysis of the sensitivity of the CFD model was carried out in this study, but based on validation results, a 20% error is estimated.In a further study, a thorough sensitivity analysis and extended validation of the CFD model is required.Yet, as we concluded in the analysis of the validation of depth and velocity on the crest and slope, the order of magnitude and trends are reasonably well predicted given all uncertainties in the model and measurement data.This suggests justified validation of the model for both crest and slope in terms of overtopping durations.
geotextile.The accuracy of the measuring devices should also be included in this analysis but unfortunately, no data was available in the literature.Yet, it is concluded that maximum velocities and maximum flow depths of the GCD model are in good agreement with the measurements taken at the Vechtdijk presented in [15].Additionally, the expected behavior presented by the GCD model (i.e., lower velocities and larger flow depths) with respect to the RCD model gives sufficient confidence for using the model in the present study.

Appendix C. CFD Validation for Overtopping Duration
The simulated overtopping durations presented in figure A3 are in good agreement with the ones measured by the surf boards on the crest and slope but not with the ones computed from the paddle wheel velocity meters.Apparently, there is a discrepancy between both measuring devices.No thorough analysis of the sensitivity of the CFD model was carried out in this study, but based on validation results, a 20% error is estimated.In a further study, a thorough sensitivity analysis and extended validation of the CFD model is required.Yet, as we concluded in the analysis of the validation of depth and velocity on the crest and slope, the order of magnitude and trends are reasonably well predicted given all uncertainties in the model and measurement data.This suggests justified validation of the model for both crest and slope in terms of overtopping durations.

Figure 1 .
Figure 1.Excess shear stress integral over time for different critical thresholds.

Figure 1 .
Figure 1.Excess shear stress integral over time for different critical thresholds.

J 28 Figure 2 .
Figure 2. Millingen wave overtopping simulator (WOS) experiment with a road on top of the crest.Volumes released by the WOS are obtained from random sampling list from a fitted Weibull distribution.

Figure 2 .
Figure 2. Millingen wave overtopping simulator (WOS) experiment with a road on top of the crest.Volumes released by the WOS are obtained from random sampling list from a fitted Weibull distribution.

Figure 3 .
Figure 3. Flow chart of the steps of the proposed methodology.

Figure 3 .
Figure 3. Flow chart of the steps of the proposed methodology.

4. 1 .
CFD ModelsDetailed two dimensional (2D) CFD Reynolds averaged Navier-Stokes equations (RANS) with k-ε turbulence closure solution models are considered better approximations of flow models where transient bed shear stresses are required with respect to uniform flow-based methods.They simulate turbulent flows for which rapid dissipation of kinematic energy is converted into internal energy in the form of turbulence.The inclusion of turbulence derived effects in the bed shear stresses becomes more important when abrupt geometric bottom changes and local surface roughness variations are present (e.g., RCD bottom profile, see Figure4).

Figure 4 .
Figure 4. CFD models boundary conditions of (A) Road over crest dike model and (B) Grass over crest dike model (schematic meshing not in actual size on both models).

Figure 4 .
Figure 4. CFD models boundary conditions of (A) Road over crest dike model and (B) Grass over crest dike model (schematic meshing not in actual size on both models).

Figure 5 .
Figure 5. Cross section of (A) Road Covered Dike RCD; (B) Grass Covered Dike GCD model domains with study points (STPs) and estimated grass cover layer thickness (0.1 m).

Figure 5 .
Figure 5. Cross section of (A) Road Covered Dike RCD; (B) Grass Covered Dike GCD model domains with study points (STPs) and estimated grass cover layer thickness (0.1 m).

Figure 7 .
Figure 7. (A) Cross section of the GCD showing two scanned profiles of the Millingen dike before (q10_sc2) and after (q50_sc5) the 0.05 m 3 /s/m experiment; (B) Available grass cover thickness ( ) calculated as the difference initial state q10_sc2 and estimated soil core, scour depths (ε) calculated as the difference between the q10_sc2 and q50_sc5 scanned profiles.

Figure 7 .
Figure 7. (A) Cross section of the GCD showing two scanned profiles of the Millingen dike before (q10_sc2) and after (q50_sc5) the 0.05 m 3 /s/m experiment; (B) Available grass cover thickness (δ) calculated as the difference initial state q10_sc2 and estimated soil core, scour depths (ε) calculated as the difference between the q10_sc2 and q50_sc5 scanned profiles.

Figure 8 .
Figure 8. Equivalence curve between critical velocity, U c , and critical shear stress, τ c , built from values reported in[14].This curve is used to include C E as a stochastic random variable by transforming the stochastic U c from Table4to τ c values, which is required as input for C E using Equation(9).

Figure 9 .
Figure 9. (A) Average kinetic energy (k(t)) of single wave volumes of 0.4 m 3 /m and 3.2 m 3 /m and (B) average turbulent erosive potential (U(t)k(t)) for single wave volumes of 0.4 m 3 /m and 3.2 m 3 /m.

Figure 10 .
Figure 10.Mean values of T (red line, right axis) and 95% confidence intervals (dashed lines) for RCD (upper plot) and GCD (middle plot).Lower plot shows ratios of RCD and GCD mean T values at each location.

Figure 10 .
Figure 10.Mean values of T (red line, right axis) and 95% confidence intervals (dashed lines) for RCD (upper plot) and GCD (middle plot).Lower plot shows ratios of RCD and GCD mean T values at each location.

Figure 11 .
Figure 11.curves for each STP for grass (left) and bare soil (right).Horizontal lines represent the grass quality reference intervals for very poor, poor, average, and good grass cover quality reported in[19].

Figure 11 .
Figure11.C E curves for each STP for grass (left) and bare soil (right).Horizontal lines represent the grass quality reference intervals for very poor, poor, average, and good grass cover quality reported in[19].

Figure 12 .
Figure 12. Results of scoured profiles for different conditions on average grass quality for (A) RCD and (B) GCD.Lower plots show their initial grass cover thickness (δ) and the eroded depths (ε) after each storm condition for (C) RCD and (D) GCD.

Figure 12 .
Figure 12. Results of scoured profiles for different q m conditions on average grass quality for (A) RCD and (B) GCD.Lower plots show their initial grass cover thickness (δ) and the eroded depths (ε) after each storm condition for (C) RCD and (D) GCD.

Figure 13 .
Figure 13.Results of scoured profiles for different conditions on with the estimated grass quality for Millingen at each STP for (A) RCD and (B) GCD.Lower plots show their initial grass cover thickness (δ) and the eroded depths (ε) after each storm condition for (C) RCD and (D) GCD.

Figure 14 .
Figure 14.Failure probabilities per storm at each STP for different grass qualities for the RCD (upper three bar plots) and for the GCD (lower three bar plots).Bars in read coincide with locations in failure state before the storm condition was routed.

Figure 13 .
Figure 13.Results of scoured profiles for different q m conditions on with the estimated grass quality for Millingen at each STP for (A) RCD and (B) GCD.Lower plots show their initial grass cover thickness (δ) and the eroded depths (ε) after each storm condition for (C) RCD and (D) GCD.

Figure 13 .
Figure 13.Results of scoured profiles for different conditions on with the estimated grass quality for Millingen at each STP for (A) RCD and (B) GCD.Lower plots show their initial grass cover thickness (δ) and the eroded depths (ε) after each storm condition for (C) RCD and (D) GCD.

Figure 14 .
Figure 14.Failure probabilities per storm at each STP for different grass qualities for the RCD (upper three bar plots) and for the GCD (lower three bar plots).Bars in read coincide with locations in failure state before the storm condition was routed.

Figure 14 .
Figure 14.Failure probabilities per storm at each STP for different grass qualities for the RCD (upper three bar plots) and for the GCD (lower three bar plots).Bars in read coincide with locations in failure state before the storm condition was routed.

Figure 15 .
Figure 15.Failure probabilities per storm at each STP for the Millingen estimated grass quality per location for the RCD (Left side bar plot) and for the GCD (Left side bar plot).Bars in read coincide with locations in failure state before the storm condition was routed.

Figure 15 .
Figure 15.Failure probabilities per storm at each STP for the Millingen estimated grass quality per location for the RCD (Left side bar plot) and for the GCD (Left side bar plot).Bars in read coincide with locations in failure state before the storm condition was routed.

Figure A3 .
Figure A3.Depth and velocity time series for the location on the crest (upper 6 plots) and the side slope (lower 6 plots) of different wave volumes from the RCD and GCD model and the measured values for the same waves during the Millingen experiment.

Figure A3 .
Figure A3.Depth and velocity time series for the location on the crest (upper 6 plots) and the side slope (lower 6 plots) of different wave volumes from the RCD and GCD model and the measured values for the same waves during the Millingen experiment.

Table 1 .
Experimental storm conditions for part I of Millingen experiment.

Experiment Part I Overtopping Discharge Interval 1 Interval Scanning Instant Profile Scan Label
The intervals shows the actual test durations and scanning instant shows the exact moment in which the laser scanning takes place.The 2 profiles are used for model the CFD model construction and validation.

Table 3 .
Characteristic values of dike configurations and 7 tested storm boundary conditions for riverine dikes for the tested range of overtopping discharges.

Table 4 .
Stochastic random variables of cover quality used as input for CE calculation.CoV is the coefficient of variation and QCF is grass quality correction factor for poor, average, and good.

Table 4 .
Stochastic random variables of cover quality U c used as input for C E calculation.CoV is the coefficient of variation and QCF is grass quality correction factor for poor, average, and good.