Numerical modelling of the erosion and deposition of sand inside a filter layer

This paper treats the numerical modelling of the behaviour of a sand core covered by rocks and exposed to waves. The associated displacement of the rock is also studied. A design that allows for erosion and deposition of the sand core beneath a rock layer in a coastal structure requires an accurate prediction method to assure that the amount of erosion remains within acceptable limits. This work presents a numerical model that is capable of describing the erosion and deposition patterns inside of an open ﬁ lter of rock on top of sand. The hydraulic loading is that of incident irregular waves and the open ﬁ lters are surface piercing. Due to the few experimental data sets on sediment transport inside of rock layers, a sediment transport formulation has been proposed based on a matching between the numerical model and experimental data on the pro ﬁ le deformation inside an open ﬁ lter. The rock layer on top of a sand core introduces a correction term in the Exner equation (the continuity equation for sediment and change in bed level). The correction term originates from the fact that the sand can only be deposited in the pores of the ﬁ lter material. The numerical model is validated against additional data sets on the erosion and deposition patterns inside of an open ﬁ lter. A few cases are de ﬁ ned to study the e ﬀ ect of the sinking of the ﬁ lter into the erosion hole. The numerical model is also applied to several application cases. The response of the core material (sand) to changes in the wave period and wave height is considered. The e ﬀ ect of di ﬀ erent layouts of the ﬁ lter is studied in order to investigate the e ﬀ ect of di ﬀ erent ﬁ lter pro ﬁ les on the resulting erosion. Finally, it is studied how much the design of a hydraulically closed ﬁ lter can be relaxed to obtain a reduction in the design requirements of the ﬁ lter thickness, while the deformation to the sand core remains acceptably small.


Introduction
Rock filters are commonly used in the design of coastal structures and scour protection to prevent erosion of the underlying material either by internal transport processes at the interface or by suction removal.The traditional design approach is to avoid any type of erosion.One such approach is the geometrically tight filter design, where a number of filter layers are placed on top of each other in such a way that the grains in the finer layer cannot pass through the pores of the coarser, adjacent layer due to geometrical constraints.This type of filter design is expensive and difficult to install, because it typically consists of a large number of layers.
One alternative approach is to use fewer layers, where the material is not geometrically constraint to pass through the pore structure of the coarser layers.The adopted approach to suppress any sediment motion in the sandy bed is to increase the filter layer thickness to such a degree that the hydraulic loading at the rock-sand interface is smaller than the threshold for mobilisation of the sediment [1].This type of filter design is termed hydraulically closed filters.The hydraulically closed filter is easy to install, but the consumption of the filter material is still considerable.The design of hydraulically closed filters relies on experimental data for the initiation of motion of sediment grains at the interface.Data for unidirectional flow [21,36] and for oscillatory flow and real waves [21,5] are available.Furthermore, Bakker et al. [1] approach the stability of the sandy base under an armour and a filter layer (hydraulically closed) from a theoretical point of view and they found that large ratios between base and filter materials can lead to a more stable bed than a small ratio.Klein Breteler et al. [21] specified the incipient motion in terms of the critical hydraulic gradient inside the granular filter, while Dixen et al. [5] defined the initiation of suction of the base material as a function of the critical mobility number, where the applied velocity was measured outside of the armour blocks.
The present work is concerned with the concept of open filters (see e.g. the experimental works by [40]and Fig. 1), where the hydraulic loading at the interface between the filter and the base material (sand) is large enough to mobilise the sediment at the interface between the base and the filter.The sand at the core is not subject to any geometrical constraints, consequently, the sand profile is deformed with time.Van Gent and Wolters [40] presented the results from an experimental campaign, where the main focus was on the response of the sand profile as a function of the material properties of the filter and the hydraulic loading.Their investigation did not include variations in water level and the variation in the wave steepness was limited.
The present work presents a newly developed numerical model, which is calibrated with and validated against the experimental data from Van Gent and Wolters [40].This model is used to study the effects like variations in the layout of the filter and changes in water level.Van Gent et al. [41] observed that the hydraulic gradient along the rocksand interface depends on the wave steepness, why it is expected that the wave steepness (wave period) must affect the erosion and deposition patterns as well; this dependency is analysed in the present work.
One of the model requirements is a numerical model that is capable of simulating wave breaking and porous media flow.The first VOFmodel that solved the Navier-Stokes equations for the free-surface flow in-and outside of permeable coastal structures was presented by Van Gent et al. [38] and they showed that this approach can provide valuable insight into the physical processes.Several other VOF models have since been developed for the interaction between waves and permeable structures [20,23,25,32,11].Jacobsen et al. [17] demonstrated that the version of the Navier-Stokes equations by Jensen et al. [20], which includes the interaction with permeable coastal structures, works successfully for breakwaters.The implementation by Jensen et al. [20] will be adopted in the present work.
Morphological modelling with detailed CFD models has also been carried out in both river and coastal problem over the last couple of decades, but the effect of a free surface on the morphodynamics [10,24,28,34,14] is not yet included on a regular basis.Morphological modelling including free surfaces is time consuming, because there is a large difference between the time scale for morphological equilibrium and quasi-steady hydrodynamics.This means that the simulation time typically exceeds the time it takes to reach quasi-steady hydrodynamic conditions and the intra-wave motion must still be resolved: the combination of these requirements results in a huge number of computational time steps.This problem was addressed by e.g.Stahlmann [34], who used a snap-shot of the hydrodynamic loading to calculate the development of a scour hole around a tripod structure over many wave periods.Subsequently, the hydrodynamics was re-calculated over the new profile of the bed and a new snap-shot of the hydrodynamics was used to evaluate the scour development.An alternative is to use the standard approach of morphological acceleration.The latter will be adopted in this work, since the hydrodynamic forcing is that of irregular waves, where the definition of a representative regular forcing (if it exists) is a research topic in itself.
As seen from the above descriptions, neither the modelling of the interaction between permeable structure and waves nor the modelling of sediment transport and the resulting bed deformation due to waves are new research areas.The combination of these two branches, on the other hand, is to the authors' knowledge not attempted before.
The present work is organised as follows: The numerical model is presented in Section 2 and the calibration of a sediment transport formula is presented in Section 3, where the validation of the numerical model is also presented.In Section 4, the numerical model is applied to different configurations of the open filters, the effect of wave heights and periods are investigated and a less strict formulation of a hydraulically closed design rule is analysed.The work is finalised with a conclusion.

Mathematical model
The mathematical model is described in this section.The numerical framework is OpenFoam [43] version foam-extend-3.1 and the framework provides the means of solving free surface flows with the volume of fluid method (VOF).Based on the discussion in Section 1, the following requirements are specified to the numerical model: (i) Solution to the Navier-Stokes equations in-and outside of a permeable layer, (ii) tracking of the free surface in-and outside of a permeable layer, (iii) generation and absorption of free surface waves and (iv) modelling of sediment transport and the resulting change in bed level inside an open filter.Each of these components is addressed below.

Hydrodynamic model
Jensen et al. [20] presented a form of the Navier-Stokes equations that accounts for the presence of permeable, coastal structures.This model has been successfully used to describe the interaction between waves and permeable coastal structures such as breakwaters [17] and to validate the wave-induced pressures inside an open filter [41].Jensen et al. [20] described the Navier-Stokes equations in terms of the filter velocity, and the Navier-Stokes equations took the following form: Here, C m is the added mass coefficient, t is time, ρ is the density of the fluid, u is the filter velocity vector, n p is the porosity of the permeable structure, p* is an excess pressure, g is the vector due to the acceleration of gravity, x is the Cartesian coordinate vector, Γ u is the diffusivity of the velocity field and F p is the resistance force due to the permeable structure.The system of equations is closed with the incompressible form of the continuity equation: The tracking of the free surface is performed with the advection algorithm termed MULES, which is the standard method available in OpenFoam.The advection equation takes the following form: Here, F is the indicator function of the VOF-field and u r is a relative velocity introduced to keep a sharp interface, see Rusche [31] and Berberovic et al. [2] for details.Note the factor n 1/ p that was introduced by Jensen et al. [20] to ensure the conservation of mass, when the fluid passes through a permeable structure.This correction term is required, because the water/air can only occupy the pore volume in the permeable structure.
The indicator function is also applied to evaluate the spatial variation of the density and the viscosity: Here, the sub-indices refer to the fluid properties for F = 0 and F = 1.
In this work, F = 1 means that the computational cell is filled with water and F = 0 means that the computational cell is filled with air.Any cell with a value of F between 0 and 1 will be located at or close to the free surface.The flow resistance term in Eq. ( 1) is described by the Darcy-Forchheimer flow resistance formulation: Here, a and b are closure coefficients that are evaluated based on the parameterisation by Van Gent [39]: Here, ν is the kinematic molecular viscosity, D n f 50 is the nominal diameter of the permeable layer, KC is the Keulegan-Carpenter number and α = 1,000 and β = 1.1 are closure coefficients.It was found by Jacobsen et al. [17] that the standard values of α and β suggested by Van Gent [39] gave good results in comparison with a range of experimental data sets on breakwaters.It was furthermore validated by Van Gent et al. [41] that the α and β coefficients resulted in correct hydraulic gradients inside of an open filter.The KC number is evaluated based on linear wave theory at the toe of the structure, where the maximum orbital velocity is an approximation to the pore velocity in the top part of the filter layer (see [17], for an extended discussion on this choice).
It should be said that there are uncertainties related to both the magnitude of β and the formulation of KC; especially due to the lack of spatial and temporal variation in KC.Both of these uncertainties will affect the results quantitatively, but no data was available for the calibration of the velocities inside of the rock layer, while a validation of the wave-induced pressures inside an open filter is already performed based on the chosen approach for β and KC [41].
The modelling of free surface waves can oftenespecially at laboratory scalebe conducted without the use of a turbulence model, because the bulk of the wave motion is described by potential wave theory.Furthermore, the boundary layers only transition to turbulent wave boundary layers at large Reynolds numbers, see e.g.Jensen et al. [19], and these high Reynolds number are unlikely in small scale laboratory experiments.The presence of a permeable structure changes the magnitude and distribution of the turbulence considerably, because the work conducted on the permeable structure by the waves equals a production of turbulence.It has been seen, however, that there is no need to add a turbulence model to predict the hydrodynamics accurately ( [20,17]).The reason for this is that the resistance coefficients a and b (Eq.( 6)) include the effect of the dissipation of wave energy due to turbulence, because the resistance coefficients are based on experimental data.Therefore, a recalibration of the resistance coefficients α and β is required, if a turbulence closure such as the k − ϵ model or k ω − SST model is added to the hydrodynamic model.Turbulence models are generally important for the modelling of sediment transport, since turbulence influences the bed shear stresses, the shape of the velocity profile and the turbulent mixing of suspended sediment.It was seen in the experimental work by Stevanato et al. [35] that the velocity profile inside of a filter layer becomes practically uniform over the depth, when the thickness exceeds a certain size, consequently, the turbulence has limited effect on the velocity profile inside filters of sufficient thickness.The two remaining propertiesbed shear stress and vertical mixing of suspended matterare not required in this work, because the sediment transport will be limited to bed load transport.A mobility number is used to describe the sediment transport; see Section 2.2 for details.This means that turbulence modelling can be excluded in this work.
It was considered to use the 'free-stream' velocity inside of the rock layer together with a friction factor to evaluate the bed shear stress.In this way, a Shield's parameter could be evaluated without the use of a turbulence model.Primarily, it is an open question whether a welldefined boundary layer even exists inside of a permeable structure, so: Are friction factors even applicable to estimate the bed shear stress inside a permeable layer?However, assume that a boundary layer does exist, the friction factors, known from e.g.oscillatory flows over smooth and rough boundaries, do not apply, since the thickness of the boundary layer is thinner in the presence of a permeable structure; see derivation in Appendix A. This analysis led the authors to completely abandon the concept of Shield's parameters for sediment transport (inside of rock layers) and a mobility number formulation was adopted instead.
The incident waves were generated with the waves2Foam toolbox for OpenFoam [13].This toolbox allows for the generation and absorption of free surface waves in a VOF model.The model relies on the relaxation zone technique, where the weighting between a target incident wave field and the computed field absorbs reflected waves.The weighting ensures that there is no accumulation of water inside the computational domain due to the inward Stokes drift.Unless otherwise specified, the irregular wave field is generated based on a linear superposition of the individual frequencies in the wave spectrum.The frequencies are non-equidistantly distributed with a higher frequency resolution around the peak of the spectrum than at the tails.This gives a better time domain representation of the exceedance statistics of wave height and wave periods for a small number of frequencies.

Sediment transport and erosion model
The sediment transport and erosion model relies on an existing implementation in OpenFoam that previously has been applied successfully to a number of cases such a cross-shore morphodynamics and sediment transport [14,16] and scour beneath pipelines [9].Details on the numerical framework are given in Jacobsen et al. [15] and the handling of mass conservation of the sand in the absence of an open filter is described in Jacobsen [18].
Only a limited amount of experimental data was found on the topic of sediment transport within granular filters.The available experimental data sets were reported in Van der Meulen [37] and the data was later re-analysed by Klein Breteler et al. [21], their Chapter 6.The two main conclusions from Van der Meulen [37] and Klein Breteler et al. [21] were: (i) the initiation of motion of the sand placed underneath the granular filter is a function of the ratio d D / n f 15 , where d is the median grain diameter of the sand and D n f 15 is the 15% quantile of the distribution curve of the filter material (see also theoretical discussion in [1]).(ii): Klein Breteler et al. [21] found that the best fit to the measured sediment transport rates was obtained with the following relationship: Here, q b is the rate of sediment transport, C is a calibration constant, v f is the filter velocity inside the granular material and v fcr is the critical filter velocity.Eq. ( 7) provides information on the functional form of the sediment transport, though C varied by a factor 4 in Klein Breteler et al. [21], so no general transport formulation was suggested.The mobility number, Ψ, is used in this work to describe the sediment transport and it is defined as: Here, ρ s is the density of the sand, g g = 2 is the acceleration due to gravity, u || is the filter velocity vector parallel to the sediment bed, s is the density of the sand relative to that of water and d is the sand size.A combination of Eqs.(7) and (8) gives the functional shape of the sediment transport that will be adopted in this work: Eq. ( 9) describes the sediment transport along the direction of the velocity vector adjacent to the bed.The units of C 1 is m 3 /m/s.The form of Eq. ( 9) has a strong resemblance to other types of bed load sediment transport formulations, e.g.Fredsøe and Deigaard [8] and Ribberink [29].In this work, a bed load type formulation is used, thus any contributions from suspended sediment transport inside of the filter are lumped into calibration coefficients (see below) and phase lags between the hydrodynamics and the sediment transport are neglected.
The critical mobility number also depends on the slope of the bed [3] as known from other areas of sediment transport, see e.g.Fredsøe [7] and Kovacs and Parker [22].The effect of the bed slope inside a filter is described in Bezuijen et al. [3] to behave as Ψ cr0 is the critical mobility number on a horizontal bed, Ψ cr is the slope corrected critical mobility number, ϕ =35 r ∘ is the angle of repose and ϕ is the slope of the bed relative to the direction of the flow.
when the flow is downhill and Ψ < Ψ cr cr 0 , when the flow is uphill.The parameters C 1 and Ψ cr0 will be calibrated in Section 3. The value of ϕ r is kept fixed in this work.
The Exner equation relates the gradient in the sediment transport field to the rate of change in bed level.The Exner equation takes the following form in the presence of a granular filter on top of the transported core material: Here, h is the vertical position of the bed level and n s is the porosity of the sand.Without the correction term n 1/ * p , Eq. ( 11) is the classical form of the Exner equation.The additional correction accounts for the presence of the granular filter, since the sediment can only be deposited in or eroded from the pores of the filter.
Note that h Δ in Eq. ( 12) is evaluated by numerical integration of Eq. ( 11), which again depends on n * p .Therefore, Eqs. ( 11) and ( 12) must be solved iteratively to obey to mass conservation of the sand.The change in porosity from one layer to the next is described as a discontinuous jump, while the porosity changes gradually in the true physical system.
The excess steepness of the bed has to be handled with some kind of sand sliding.This can either be done with an increase in the downslope bed load sediment transport [30], solving for the actual sand sliding with a time dependent sand slide equation [4] or with the help of a geometrical approach ( [26,27]).The latter has been chosen in this work and an implementation for unstructured meshed in three dimensions is applied; see Jacobsen [18] for details.The method by Jacobsen [18] simplifies to that of Niemann et al. [27] in two dimensional simulations.The sand sliding routine accounts for the effective porosity to avoid mass conservation errors in the volume of the sand bed.
The resulting rate of bed level change affects the solution to the hydrodynamics through a displacement of all vertices in the computational mesh.The mesh motion is calculated with a simple linear interpolation between the location of the bed and the top of the computational mesh.

Sinking of the filter
It was outlined above that the presence of the filter is important for the rate of change in the bed level.Therefore, the temporal behaviour of the filter is equally important to address.It was found early in the model development that the sinking of the filter into the erosion hole(s) is important for the bulk behaviour of the developing profile.If the sinking was not included, there was no flow resistance in the momentum equation below the original position of the filter.It had the consequence that the flow speed increased in computational cells without flow resistance and the larger flow velocities transported sand to the sides.The profile of the sand core developed into a wiggled profile, where the wiggles had a short wave length, and the profile did not have any of the characteristics observed in the experiments.Allowing the filter to sink into the erosion hole created a more realistically profile.The outer edge of the filter was kept fixed in the numerical model, because the tracking of the outer edge of the filter would have increased the complexity of the model.
Due to the mesh motion (see Section 2.2) the computational cells inside of the filter layer had to be re-found every time step.The amount of flow resistance at the top of the filter was changed gradually each time step to avoid shocks in the numerical model, which would otherwise have appeared, if a cell had no flow resistance in one time step and a full implementation of the flow resistance in the following time step.The gradual change in the flow resistance is based on the fraction of the computational cell that is inside of the filter.The interface of the filter is thought of as a discontinuous change in flow resistance.The fraction of the cell inside of the filter is found through a cutting routine, which is based on the geometrical cutting routine presented in Jacobsen et al. [13].The weighted flow resistance in a computational cell shared by a rock layer and water is a simple volume averaging between no resistance and the resistance coefficients for the rock layer; this is identical to the method described in Wellens et al. [42].The weighted flow resistance between two rock layers (e.g.filter and armour layers) is evaluated based on an equivalence principal, where it is demanded that the simulated pressure drop over the computational cell drives the same average flow velocity as would have been predicted, if the two layers were modelled separately.The latter approach allows for a better representation of complex layouts of filters and armour layers on simple computational meshes.

Boundary conditions
The bottom boundary (the sea bed) is described with a slip velocity boundary condition, where the slip velocity boundary condition is corrected for the fact that the boundary is moving.The shoreward and vertical boundary on top of the filter is treated as an open boundary, where water and air can exit the computational domain, while only air can enter.The upper (atmosphere) boundary condition imposes an atmospheric pressure corrected for the square of the velocity.The air and water can freely exit, while only air can flow in.

Calibration and validation
This section focuses on the calibration and validation of the sediment transport formulation.Validation of the hydrodynamic model for open filters is already presented in the work by Van Gent et al. [41].The validation was based on measured 2% exceedance values of the hydraulic gradient, i 2% .
All the simulations in this section are performed on a computational grid that was generated with the meshing tool snappyHexMesh, which is distributed with OpenFoam.The discretisation (except close to the slope) was kept uniform at 0.015 m times 0.015 m.This resulted in a total of 73,000 computational cells for the 1:4 slope and 127,000 computational cells for the 1:7 slope (see details below).

Calibration
Two experimental test series from Van Gent and Wolters [40] and one additional (unpublished) test series were used for the calibration of the sediment transport formulae.These test series are referred to as Test cases S04, T09 and T001; see Table 1 for 2 .H m0 is the spectral wave height and T p is the peak wave period.Data set T001 consisted of one set of wave parameters and the profile of the sand core was measured after 3, 6, 9 and 18 h of wave exposure.
The task to measure the developed profile was time consuming, because the entire filter layer had to be carefully removed down to the top of the sand profile.5 or 6 individual profiles were measured over the width of the flume so lateral variability could be quantified.
The calibration procedure for the sediment transport formulation in Eq. ( 9) was performed in a two-stepped approach: 1. Evaluate the mobility number with time for all test cases along the initial slope with the hydrodynamic part of the numerical model.Use this information to approximate the coefficients C 1 and Ψ cr0 for each case.This approach decouples the hydrodynamics from the developing profile, thus the method is merely to be trusted as an initial estimate.This approach is termed "the linear approximation".2. Apply the complete and nonlinear model with the feedback between hydrodynamics and the developing profile.The initial guess for C 1 and Ψ cr0 was based on the linear approximation.Four additional simulations were performed for each validation case with a change of ± 10% in C 1 and ± 30% in Ψ cr0 .
The experimental and numerical profiles of the sand core were used for the calibration.The calibration was based on integrated properties of the profile, where the possible properties are the area of the erosion hole and the area of accretion.A visual comparison between the profiles from the experimental and numerical tests is also presented.The erosion, A e s , , and accretion, A acc , areas are defined as follows (notation adopted from [40]): equals the porosity of the filter layer.This suggests that the compaction of the initial profile of the sand core is larger than the (initial) compaction of the deposited material.Since gradual compaction of the deposited material is not included in the numerical model, the area of the accretion cannot be used as a calibration parameter.The area of the erosion hole, A e s , , is assumed to be unaffected by the gradual compaction and A e s , is therefore used as the target value in the calibration and validation processes.Loss of sediment from the profile could also explain the behaviour of A A / e s acc , with time, but no significant loss of sediment was observed in the flume (offshore deposition) and photos taken during the tests documented that the water column was clear.Following the linear calibration approach, C 1 was estimated to be 10 −5 m 3 /m/s for all cases, while Ψ cr0 depended on the grading of the filter material, σ ; see Table 1.Ψ =0.115 cr0 for the cases T09 and T001 and Ψ =0.090 cr0 for test case S04 based on the linear calibration.The nonlinear development of the sand profile was evaluated with a morphological acceleration factor of 20, thus the full duration of 3 h of wave loading in the experiment corresponded to 540 s in simulated time.It was checked that the resulting dimensions of the erosion hole differed less than 5% between 7 realisations of the incident irregular wave field.Here, a realisation of an irregular wave field refers to a specific set of random phases.It was concluded that one realisation of an irregular wave field was sufficient.
The full set of test cases used for the final nonlinear calibration is tabulated in Table 2.The corresponding hydrodynamic loading is given in Table 1.

Table 1
The properties of the 3 test cases used to calibrate the sediment transport formulation.(*) This case was repeated 6 times and the development in the profile of the sand core was measured after 3, 6, 9 and 18 h of wave loading.Table 2 The sets of calibration coefficients for the sediment transport formulation.The coefficients for T09 and T001 are identical.Each cell contains the following information: The variation of A e s , with time for the 3 calibration cases is depicted in Fig. 3. From the Test cases S04 and T09 it is clearly seen that the increase in wave loading increases the rate of erosion; this is expected.For Test case T001 there is a constant forcing and an almost linear increase in A e s , is observed.Two ways are used to evaluate the area of the erosion hole from the experimental data:

Evaluate A e s
, based on the mean profile of the sand core, where the mean profile is based on the 5-6 profiles measured across the width of the flume.

Evaluate A e s
, as the mean of the area of the erosion hole for each of the 5-6 profiles.
There are only small differences for most of the adopted data sets (including the validation in Section 3.2), but for case T09 there is a considerable difference of 20%.Consequently, the mean of A e s , from the individual profiles are used as the most robust estimate of the mean eroded area.
The calibration Set A is the best candidate for T09 and T001, while both Set B and Set E qualify as a calibration set for S04.Recall that the filter properties were the same for T09 and T001 (wide grading of the filter material).The profile development of the sand core for the reduced number of calibration sets are shown in Fig. 4. The comparison of the profile from experimental and numerical data is good for T09 and T001.From an inspection of the developed profiles, calibration Set B is chosen for S04 (Fig. 4A).For the profile of 1:7, the location of the erosion hole is well captured, whereas the model has a consistent bias towards a more seaward located erosion hole in the numerical results in comparison with the measurements for a slope of 1:4.Again, compaction of the deposited material is less than the eroded (see Fig. 2), hence the accreted area is larger for the experimental data in comparison to the numerical results.The main reason for the discrepancy is that a gradual compaction of the deposited soil is not included in the numerical model.
The transport formula in Eq. ( 9) contains two free parameters: C 1 and Ψ cr0 .The parameters ϕ r and the exponent 1.5 in Eq. ( 9) should also have been included in the analysis, but the available data did not justify an expansion of the set of free parameters.

Validation
The validation of the sediment transport formulation is presented in this section.A total of six additional experiments from Van Gent and Wolters [40] were selected for the validation.Three cases had an initial slope of 1:7 and three cases an initial slope of 1:4.One case had both a filter and an armour layer.The details on the material properties for the filter and armour layers are presented in Table 3.
Each validation case consisted of a number of wave series.The wave properties of these wave series are summarised in Table 4.
The calibrated transport coefficients from Section 3.1 are used in the validation cases.The exception is case S10 that had a different size of the filter material: the median filter material was D =0.038     et al. [21].It was estimated that the critical mobility number would be approximately a factor of 2.0 smaller, so a value of 0.055 was used.The values 0.045 and 0.065 for Ψ cr0 were also tested, thus case S10 can partly be considered as a calibration case.
The resulting profiles for Test cases T06, S01 and S10 are depicted in Fig. 5.All 3 values of the critical mobility number for validation Test case S10 are included in the figure.The results show that the calibrated sediment transport formulation for S01 and S10 yields a good comparison between the experimental and numerical results, though the sediment transport is less diffusive in the numerical model than in the experiments.The results for T06 show an under-prediction of the size of the erosion hole, while the location of the erosion hole is well captured.The results for T06 also exhibit less horizontal diffusion of the profile development.This is likely due to the lack of e.g.suspended sediment transport.
The integrated properties such as A acc and A e s , are depicted in Fig. 6 along with the maximum erosion depth and the maximum accretion height.The panels in Fig. 6 depict a one-to-one comparison with the experimental data.The data from the calibration is also included in Fig. 6.It is seen that the erosion area is well predicted (except for the validation cases for a slope of 1:7).The accretion area is generally under-predicted by the model and this is attributed to the gradual compaction of the deposited material in the physical experiments.
The depth of the erosion hole, z s , and the height of the accretion, z acc , exhibit more scatter (Fig. 6), but this is likely caused by the fact that they are extreme values, whereas the areas are integrated quantities.While z s is either well or under-predicted, the value of z acc is either well or over-predicted.The over-prediction of z acc represents the peaked behaviour of the accretion area, which suggests that a diffusive sediment transport mechanism is missing.This diffusive sediment transport mechanism could for instance be phase lags between hydrodynamics and suspended sediment transport, where any phase lags are neglected in this work.Such phase lags are known to have an effect on the morphological development of for instance river dunes [6] or the net suspended sediment transport in the surf zone [16].

The effect of a sinking filter
The results presented in Sections 3.1 and 3.2 depend on the assumption that the outer profile of the filter layer is fixed in time.
This assumption leads to a gradually increasing thickness of the filter layer in the numerical model with an increase in the erosion depth.In this section, the last wave series in Test case T06 (see Table 4) has been analysed for the effect of the deformation of the outer filter.This test case will be termed T06-4 in the following.
The sinking of the filter was measured between each of the wave series in the experimental campaign by Van Gent and Wolters [40].Consequently, information of the amount of sinking before and after T06-4 was available.In the numerical model, the shape of the filter is  deformed to mimic the shape both before and after T06-4.The two profiles of the filter are representatives for the two extremes of the filter geometry during T06-4.The effect of this modified filter geometry on the resulting erosion and deposition is evaluated.The profile of the sand core and the three filter configurations at the beginning of T06-4 are depicted in Fig. 7A.
The resulting profiles after 3 h of wave exposure are depicted in Fig. 7, where it is seen that the straight filter caused continued erosion, while the deformed filters gave rise to a completely different response.It is seen that the deformed filtercounter intuitivelyresulted in a smaller effective erosion, where the initial erosion hole was partly backfilled and the erosion took place more shoreward.The shoreward movement of the location of erosion makes sense, since the settling of the filter is shifted shoreward relative to the sand profile.The smaller erosion depth, however, requires a careful investigation, which will be based on an analysis of the sediment transport and the mean velocity parallel to the sand-rock interface.
First of all, the 2% exceedance value of the hydraulic gradients parallel to the sand profile was investigated, but the shape and maxima were largely unaffected by the change in the geometry of the filters (results not shown).Therefore, the mean bed load rate and the mean velocity field were evaluated and these quantities are depicted in Fig. 8.While the original (straight) layout of the filter layer resulted in a single peak in the mean sediment transport, the deformation to the filter layer gave rise to a much broader 'peak' (see Fig. 8B).This means that the gradient in the mean sediment transport rate is lowered and the smaller gradient corresponds to a decrease in the rate of bed level change.This behaviour was linked to the mean velocity distribution adjacent to the bed inside of the filter (Fig. 8A) that changed considerably with the deformation of the profile of the filter.The decrease in the offshore directed mean flow is thought to be the cause for the change in the erosion and deposition patterns.Consequently, the sinking of the filter acts as a stabilising effect on the rate of erosion or at least stabilises the maximum erosion depth.The feed-back between the filter profile and the resulting erosion of the sand core is investigated in Section 4.2 in conjunction with a possible optimisation of the configuration of the filter.
A derived effect of this analysis is that the sediment transport formulation as calibrated in Section 3.1 does not account for this settling, thus the sediment transport formulation is only applicable for a limited deformation of the filter.Application of the model for cases with large deformations to the filter layer requires careful interpretation of the results.The present results furthermore show that the sinking of the filter does not explain the under-prediction of e.g.A e s , or z s , see Fig. 6.

The effect of wave height and wave period
The decay of the hydraulic gradient through a rock layer is known to behave exponentially with the wave number scaling, see e.g.Yamamoto et al. [44] and Hsu and Jeng [12].A consequence hereof is that a dependency of the wave period is expected on the resulting deformation to the sand core: a larger wave period will result in a larger response of the sand core.
That a larger response of the sand core is to be expected under   larger wave periods can also be approached based on the empirical formula for the hydraulic loading at the interface by Van Gent and Wolters [40].Van Gent and Wolters [40] stated that the 2% exceedance value of the hydraulic gradient at the rock-sand interface is given as Here, i 2% is the 2% exceedance value of the hydraulic gradient and d d = eq f is the equivalent filter layer thickness, where d d = eq f holds for single rock layers.The experimental data by Van Gent and Wolters [40] was mostly limited to a wave steepness of 0.040 and Eq. ( 14) was reproduced numerically by Van Gent et al. [41] for the same value of s p .Van Gent et al. [41] also computed i 2% numerically for incident waves with different values of the wave steepness (s p equal to 0.015, 0.020 and 0.030) and they modified Eq. ( 14) to account for the steepness in the following manner: Eqs. ( 14) and ( 15) are identical for a steepness of s =0.040 p .Eq. ( 15) essentially states that for a given wave height and a given thickness of the filter layer, the hydraulic loading at the interface between the rock layer and the sand core increases with a decrease in the wave steepness, i.e. an increase in the wave period.This interdependency can be investigated using the developed numerical model.The effect of the wave steepness on the change in the profile is illustrated with a set of 7 simulations, where the filter properties and the initial profile is the same as in test case T001, see Table 1.The results from T001 will be referred to as the reference solution, where T =1.38 p s and H =0.120 m0 m.The magnitude of T p and H m0 were adjusted, such that values of s p equal to 0.0150, 0.0225 and 0.0300 were obtained.This gave 6 simulations: three obtained by a change to the wave period and three obtained by a change to the wave height.The wave height ranged from 0.045 m to 0.120 m and the wave periods ranged from 1.38 s to 2.26 s.The computational domain is the same as outlined in Section 3.
The resulting profile development is shown in Fig. 9A after a total of 3 h of morphological development.The reference profile after 6 h of development is also depicted in Fig. 9A.The deformation to the sand core was as expected, namely a larger rate of erosion in the case of a larger wave period.It is seen that a larger wave period also affects the width of the erosion hole: the width increases with an increase in the wave period.
The relative size of the erosion hole is depicted in Fig. 9B, where the reference solution for s =0.040 p is used for the normalisation.This shows that an increase in the wave period (smaller steepness) increases the area of the erosion hole up to a factor of 2.5, while a decrease in the wave height gives rise to a smaller erosion hole.No difference between the normalised erosion area after 3 h and 6 h was observed.A few of the tests from Van Gent and Wolters [40] have been re-analysed in this work, since they had identical specifications except for the steepness that was 0.015 and 0.040 respectively.It was seen that a smaller wave steepness gave rise to a 1.50 times larger erosion hole for a thickness of the filter layer of 0.20 m, while there was merely an increase in the area of the erosion hole by a factor of 1.05 for a 0.40 m thick filter layer.Hence, the experimental data do suggest an effect of the wave steepness, although the number of tested conditions and the difference in erosion areas were such that the influence of the wave period could not be incorporated in the developed design guidelines.
Van Gent and Wolters [40] found that the area of the erosion hole scales with i 2% and the slope of the structure.The dependency on i 2% is investigated in this work by keeping the value of i 2% fixed for different values of the steepness.The corresponding values of H m0 and T p are evaluated by combining Eq. ( 15) and the definition of the steepness (s H T = /1.56 2 ); the target solution of i 2% is given from the reference case with T =1.38 p s and H =0.12 m0 m.The final sets of wave parameters are presented in Table 5.
The normalised area of the erosion hole is included in Fig. 9B for the additional simulations presented in Table 5.It is seen that the relative size of the erosion hole is the same for all of these simulations.It appears that the results from the numerical model areat least qualitativelyin line with the experimental observations by Van Gent and Wolters [40].

Effect of the KC-number on the results
The sensitivity to the KC-number on the profile development will be analysed in this section.It was discussed in Section 2.1 that there are uncertainties related to the KC-number in the resistance formulation for permeable structures, since a single, global KC-number is pre- scribed for each permeable layer.The results presented in Fig. 9 were based on a KC-number that varied with the incident wave properties based on the method outlined in Section 2.1.
In order to evaluate the effect of the KC-number, the simulations presented in Fig. 9 were re-executed with a constant value of KC = 7.40.The results are presented in Fig. 10.The reference value of A e s , (used for normalisation) is the same in Fig. 9 and Fig. 10.The most noticeable difference between Fig. 9 and Fig. 10 is that the effect from the wave period decreases, i.e. the normalised size of the erosion hole becomes smaller.This shows that a better approximation of the KC-number throughout the structure is important, when the response of the sand core is evaluated, because the sediment transport is governed by the filter velocities close to the profile and not the hydraulic gradients (the former are much more dependent on the exact values of the resistance coefficients than the latter).A more realistic description of the KC-number, such that it becomes a function of space

Table 5
The wave properties for the additional simulations with a fixed hydraulic gradient.The asterisk marks the reference case.s p [-] i 2% [-]  H m0 [m

The effect of the configuration of the filter
In this section, the effect of the configuration (layout) of the filter on the resulting deformation of the sand core is evaluated.First, some idealised profile configurations of the filter layer are analysed to obtain an overall understanding of the response of the sand core due to changes to the filter.The effect of different water levels on the same filter configuration will also be analysed.Secondly, the response of the sand core beneath a naturally shaped filter layer is evaluated.The different natural shapes were taken from existing laboratory experiments conducted at Deltares.

Deformation of the sand core beneath idealised profiles of an open filter
It was seen in Section 3.2.1 that a local reduction in the slope of the filter layer gave rise to a re-distribution of the erosion.This effect will be analysed in this section for idealised profiles of the filter layer.The reference case is again Test case T001 with a slope of 1:4 of the sand core.Modification to the slope of the filter around the still water level is introduced into the numerical model and the slope above and below the modifications remains 1:4.Three classes of layouts (configurations) are presented in Fig. 11.In Fig. 11A, the modification is such that the thickness of the filter is d =0.20 f m at z = 0 m and five local slopes (1:2, 1:3, 1:4, 1:5 and 1:7) are introduced around the intersection of the still water level.This approach gives rise to a change in the thickness of the rock layer on the lower and upper parts of the profile.In Fig. 11B, the thickness of the rock layer is d =0.2 f m from the toe up to z = 0 m and only the profile above z = 0.0 m has a varying filter layer thickness.This class of configurations was included to equalise the dissipation over the lower part of the rock layer.Finally, the configurations in Fig. 11C are similar to Fig. 11B, but the change in slope begins at z = −0.1 m instead of 0.0 m.Only the slopes 1:4, 1:5 and 1:7 were analysed for this case.Different water levels were simulated and the main observations are described per class of configurations.All simulations had a morphological duration of 3 h and the morphological acceleration factor was 20.The computational domain is the same as outlined in Section 3. 11A.Three water levels were tested: z = −0.05m, z = 0.0 m and z = 0.05 m.The resulting deformations to the sand core are depicted in Fig. 12 and the results from the straight profile (slope 1:4) are also included.It is seen that the amount of erosion decreases with an increasing thickness of the filter on the lower part of the slope.This is ascribed to an increase in the dissipation of the wave energy with an increase in the layer thickness on the lower part of the slope.For a change in water level, it is observed that the horizontal movement of the location of maximum erosion is the largest for the 1:7-modification and the smallest for the 1:2modification.Therefore, it could be worth to use a geotextile to protect the most exposed part and use the layout of the filter to geometrically constrain the location of maximum erosion.Edge effects on the geotextile can be important.11B.Three water levels were tested for this configuration.These were z = 0.00 m, z = 0.05 m and   10.As in Fig. 9, however, the KC-number is the same for all cases and equals that of the reference case (KC = 7.40).z = 0.10 m.It was unnecessary to model water levels below z = 0.00 m, because the profile developments would largely be identical.The results are depicted in Fig. 13.

Configuration B from Figure
Contrary to the profiles in Fig. 11A, the dissipation over the lower part of the profile is the same, so it is easier to evaluate the effect of the various modifications to the change in profile.Again, it is seen that there is a more shallow erosion hole with the 1:7-modification and the height of the accretion is also smaller, thus there is a reduced risk of suction of the deposited material into the main part of the water column.The geometrical control of the location of maximum erosion is also predicted for this layout of the filter.The erosion and deposition is smaller for the 1:2 and 1:3 slope configurations in comparison with Configuration A, which is due to the fact that the thickness of the filter is larger than for the configuration in Fig. 11B.This also explains the larger erosion for the 1:5-and 1:7-modifications for Configuration B, since the thickness of the filter layer is smaller than in Configuration A.
For a water level at z = 0.00 m, there is only a minor difference in the results between the various layouts, since the differences between the profiles are all shoreward of the intersection between the still water level and the filter layer.This shows that the geometrical constraint on the location of the maximum erosion only has an effect as long as the modification to the slope begins below still water level (Fig. 13).11C.The tested water levels for Configuration C were z = −0.10m, z = 0.00 m and z = 0.10 m.The accumulated bed level change is depicted in Fig. 14.The interesting observation related to these simulations are that there is a natural limit to the usefulness of a gentle, intermediate slope (see Fig. 14C), where a large erosion is predicted for the 1:7 slope for a water level of z = 0.10 m.The results underline the drawback of the gentle modification, namely that the horizontal translation of the location of the maximum erosion is controlled by the local slope: The location of maximum erosion moves 1.5 m under a water level change of 0.2 m (1:7modification), while the location of the maximum erosion for the reference profile (1:4 slope) merely moves 0.85 m under the same change in the water level.The magnitude of the translations matches the ratio between change in water level and the slope.4.2.1.4.Summary.It was seen that it is possible to change the magnitude (depth) of the erosion hole by locally introducing a gentle slope.This approach, however, has a naturally limit, since the thickness of the filter becomes too thin for certain water levels (see e.g.Fig. 14C).This leads to increased erosion.Another drawback of the gentle modification is that the location of maximum erosion displaces considerably with a change in the water level.

Configuration C from Figure
A locally introduced steeper slope in the profile has the benefit that the location of maximum erosion is geometrically constrained for different water levels.A design that involves this type of filter layout together with a strategically placed geotextile might reduce the magnitude of the erosion to an acceptable magnitude.
An additional note is that the energy dissipation is of importance for the response of the core material, consequently, the addition of an armour layer (or coarser material) on top of the filter will also change the response of the core material.

Deformation of the sand core beneath a naturally shaped open filter
A set of test cases were performed in the 3D wave basin at Deltares, where the deformation of a cobble beach was studied.The laboratory study was performed in a 1:30 scale and the alongshore transport of cobbles was measured based on the change in the cross-shore profile, see examples of the profiles in Fig. 15.The core was a concrete profile with a constant slope of 1:3.5.The alongshore transport was measured under various conditions of which the 1-in-50-year event was adopted.The prototype conditions were H =3.8 m0 m and T =10.1 p s in a water depth of 12.4 m.A JONSWAP spectrum with a peak enhancement factor of 3.3 was used.The meshing approach is the same as outlined in Section 3, but due to the difference in scale, the resolution was set to 0.2 m times 0.2 m.There were a total of 119,000 computational cells.
In this work, the numerical model is used to investigate what would have happened to the core material, if it had been constructed of sand with a median grain diameter of 0.3 mm.The cobbles had a median (prototype) diameter of 0.15 m and a geometrical spreading of approximately D D / =2.0 . The initial profile and the two deformed profiles shown in Fig. 15 were used and the simulations were performed in prototype scale.It was discussed in Section 3.2 that Ψ cr0 depends on the ratioD d /   The test duration in the laboratory was 6 h, thus a morphological acceleration of 10 was chosen for these simulations.The simulated duration was 2160 s.The reason to choose a morphological acceleration of 10 over the factor of 20, which was used for the majority of the simulations in this work, was to have a fair amount of individual waves.The resulting deformation of the core is depicted in Fig. 15.Not surprisingly, the largest deformations were found for the profile at the up-wave edge of the cobble beach (Fig. 15B), while the smoothed, deformed profile in the middle of the cobble beach hardly exhibited any deformation.The initial profile, however, also exhibited large deformations and these are attributed to the kinks in the filter layout at x = −32 m.
In order to analyse the effect of the kinks in the filter, the distribution of the hydraulic gradient was computed over a fixed core (no bed changes) for the same wave loading.The 10% exceedance value of the hydraulic gradient, i 10% , was evaluated instead of i 2% due to the relatively smaller number of waves.The results are depicted in Fig. 16 and it is clearly seen that the maxima in i 10% for the various layouts of the filter corresponds almost exactly to the change from erosion to deposition as depicted in Fig. 15.The maximum in i 10% for the initial profile is located shoreward of the berm and this location is hypothesised to be related to the presence of the kink in the configuration of the filter.The outer flow must change direction in the kink.This change in direction will cause a local pressure peak, which consequently displaces the maximum in the hydraulic gradient to the sides.

Hydraulically closed filters
The design of an open filter is a balance between an acceptable level of deformation of the sand core and the construction costs.The resulting deformation of the rock layer(s) must also be taken into account.The deformation to the rock layers occurs due to sinking into the erosion hole and damage due to external hydraulic loading.For an increasing thickness of the open filter, the filter turns into a hydraulically closed filter which is defined as a geometrically open filter for which no erosion and deposition can occur due to a limited hydraulic loading at the rock-sand interface.Consequently the deformation due to sinking vanishes for hydraulically closed filters.A hydraulically closed filter is obtained, when i i max < cr (16) at the rock-sand interface.For engineering purposes, this condition is relaxed to the following form: Here, a scaling factor α i is introduced.Besides the practical usage of an exceedance value of the hydraulic gradient in place of i max , α 1< i can provide a less conservative (but still safe) design.It was observed in the data by Van Gent and Wolters [40] their Fig. 12 that there was only a noticeable amount of erosion, when α i exceeded a value of 2-3.This essentially means that the dimensions of a filter layer can be decreased under a certain design condition, if α i can be taken larger than 1.0.
In this section, the analysis of hydraulically closed filters will be based on the ratio between i max 2% and i cr , where both of these quantities are linked to the filter velocities through the resistance coefficients of the rock layer.
The methodology is as follows: • Evaluate the changes in the profile after 6 h of morphological development for a set of wave heights and wave periods on a given profile; details are given below; • Decide on a criterion of a maximum deformation, e.g.D 0.1 n f 50 , that is considered as an acceptable deformation; • Evaluate the exceedance value of i 2% (over the profile) and compare i max 2% to either i cr or i cr0 with a measure for the deformed profile.
Irrespectively of the uncertainty related to the influence of the KC-number on the developing profile (see Section 4.1) it is still possible to perform this analysis, because the near bed velocities are transformed into hydraulic gradients, i.e. transformed into quantities, which are only weakly dependent on the actual values of the resistance properties.
The critical mobility numbers Ψ =0.115  critical filter velocities of u =1.8 cr0 cm/s and u =1.4 cr cm/s, respectively.Insertion of these values into the Darcy-Forchheimer resistance formulation yields the critical hydraulic gradients.The resistance coefficient b depends on the incident waves through the KC-number, thus so does the critical hydraulic gradients for each of the test cases in this analysis.
Based on the results presented in Fig. 9, it was found that the critical value of i was in the interval 0.02-0.09.Consequently, sets of wave periods and wave heights were found inside of this interval with use of Eq. (15).Eq. ( 15) was solved for 8 values of i ∈ [0.02,0.09], 5 wave periods (1.40 s, 1.70 s, 2.00 s, 2.30 s and 2.6 s) and three values for the thickness of the filter layer: 0.20 m, 0.30 m and 0.40 m.Eq. ( 15) had no solution for a few of the above parameter combinations.A total of 115 simulations were performed.The remaining properties were kept fixed: n =0.38 m, the water depth at the toe h = 0.85 m and a slope of 1:4.The computational mesh is the same as outlined in Section 3 with a mesh resolution of 0.015 m times 0.015 m, but to accommodate for larger wave periods, the domain was made longer.This resulted in a computational mesh with a total of 96.600 computational cells.
Each case was simulated for 1080 s with a morphological acceleration of 20, which corresponds to 6 h of morphological development at laboratory scale.The processing of the data was as follows: • Evaluate the accumulated bathymetrical change over 6 h and obtain the maximum erosion depth (z s ) and the maximum accretion height (z acc ); • Evaluate the areas of the erosion and accretion regions; • Evaluate the hydraulic gradients based on the near bed velocities and the resistance coefficients; • Determine i 2% and i max along the entire slope in points, where more than 100 individual waves were encountered.This excludes part of the "swash" region; • Evaluate the critical hydraulic gradients based on u cr0 and u cr and the resistance coefficients a and b.There is a discrepancy between the present work and Van Gent and Wolters [40] and Van Gent et al. [41].The last two references used an averaged value of i 2% over the erosion area ( i E 2% ) (no hydraulic gradients could be measured in the accretion zone), while the maximum value over the entire profile was used in this work: i max 2% .For large values of i max 2% , the ratio 2% is 0.6-0.8.This is utilised below.
The non-dimensional variation in A acc and A e s , are depicted in Fig. 18, where two normalisations were used:   that there is a dependency on the thickness of the filter layer, while the normalisation by d D f n f 50 (panels B and D) makes the data come onto a common curve.It is unknown at this stage, whether the dependency on d f observed with a normalisation with d f 2 is caused by the fact that the present (numerical) data set is much larger than the experimental data set in Van Gent and Wolters [40] or caused by a bias in the numerical model.
In Fig. 18A and Fig. 18C the empirical expression taken from Van Gent and Wolters [40] is also included, where their expressions read: Here α cot =4 is the slope of the initial profile and K =0.m, σ = 6.2 in model scale of 1:5.5).The filter was placed on top of a mobile sand core with a median grain diameter of 0.21 mm.The deformation of the sand core was investigated in the physical experiments and it was found that the sand core was stable throughout the experimental campaign, i.e. there was no observable deformation.
In this section, the observations from the physical experiments will be addressed in two fashions: (i) a simple evaluation of the critical and actual gradients along the sand profile and (ii) numerical simulations of the setup with the numerical model.
The value of a is readily evaluated to 33.4 s -1 with an estimated value of n =0.38 p .The value of b requires some discussion.The KC number is first evaluated based on the linear wave theory at the toe of the structure (see Section 2.1), which yields such a large value of KC that the KC-correction to b vanishes and b = 860 m -1 is obtained.However, the present analysis is relevant at the rock-sand interface at the instance, where the critical velocity is exceeded.Utilising u cr for the evaluation of KC yields KC = 10, which gives b = 1500 m -1 .The two bounds on b gives a critical hydraulic gradient in the interval 0.057- 0.068, where the upper value is more likely, since it is based directly on the critical velocity at the rock-sand interface.
The measured 2% exceedance values of the hydraulic gradients for the experiments were reported in Van Gent and Wolters [40], their Fig.18.Values up to 0.27 were reported.Based on the estimated values of i cr , i i / cr 2% is within the interval 4.0-5.0.The lower value is expected to be the most realistic (largest value of i cr ).Comparing i i / =4 cr 2% with the results in Fig. 17, it explains, why no deformations to the sand core were observed in the laboratory experiment.

Numerical evaluation of the erosion and deposition
patterns.The layout of the physical experiments is presented in Fig. 19.The original (non-deformed) profile of the gravel was chosen, since the gravel was eroded at the lower part of the slope (well below still water level) and was deposited more shoreward.Consequently, the original profile has the thinnest layer of gravel around the location, where the largest amount of erosion of the sand core is expected (though, for these simulations, no erosion is expected at all).The simulations are performed at model scale to match the discussion above.
The most severe wave condition was chosen for this work: H =1.44 m0 m and T =5.76 p s with a water level at 0.92 m.The waves were prescribed with a JONSWAP spectrum with a peak enhancement factor of 2.2.The material properties for the gravel layers and the sand are already presented above.During the physical experiments, several configurations of the sand and gravel were considered.The difference was that the amount of deposition of sand inside of the gravel layer was studied, i.e. the gravel was partly buried in sand.In this work, two of configurations from the physical tests are considered (see Fig. 19): (i) The gravel is placed on top of the sand without any deposition of sand inside the pores of the gravel and (ii) the gravel is placed on top of the sand and the pores in the lower half of the gravel are filled with sand.The latter is simply included in the numerical model by prescribing the location of the sand core to be at the transition from gravel to the sandfilled gravel, i.e. erosion takes place inside the gravel.The duration of the laboratory experiments was 1.28 h and the simulations were executed for 920 s (morphological acceleration of 5) and 2300 s (morphological acceleration of 2).The computational mesh has a resolution of 0.075 m times 0.075 m and a total of 137,500 computational cells.
The profile development for the two configurations (with and without sand-filled gravel as initial condition) is depicted in Fig. 20.The results show that there is a limited profile development (less than 1 cm) and these results are in line with the observations made in the experiments.There is a slightly larger erosion and deposition for the sand-filled gravel, which is caused by the following: (i) erosion takes place in the gravel, i.e. there is a visually larger amount of erosion relative to the deposition and (ii) the sand-rock interface is closer to the top of the filter, so there will be less damping the waves.There is an abrupt peak (though still small) in the accumulated change in the profile at x = 26 m in Fig. 20A, and this effect is attributed to the change in slope of the filter at this location.The results also show that there is no quantitative difference between simulations with morphological accelerations of 2 and 5.The results from this study support that the critical hydraulic gradient can be relaxed by a factor α 1< i .In this case α =4 i .0appears to be a realistic estimate.

Conclusions
A novel computational framework has been presented that allows for the evaluation of the deformation of a sand core inside a permeable coastal structure.The layer of rock on top of the sand core required the implementation of a correction term in the Exner equation; a correction term that originates from the fact that the sand can only be deposited in the pores of the layer of rock.Due to the few experimental data sets on sediment transport inside of rock layers, a sediment transport formulation was developed based on a matching between the numerical model and experimental data on the profile deformation inside open filters.
Irrespectively of some discrepancies between the measured and predicted erosion and deposition patterns, it was found that the model could be used to perform parameter studies on the deformation of the sand core underneath an open filter.These applications span the effect of the incident wave properties on the resulting deformation of the sand core, the effect of changes to the profile of the filter and the effect of changing the water level relative to the filter.For non-uniform filter configurations (varying thickness), it was seen that a certain configuration would be beneficial for some water levels, but that the same configuration led to increased erosion for a different water level.These numerical results showed that there are a tight interplay between the filter and the sand core and the existing design guidelines can only be used in pre-design for complex configurations of the core, filter or armour layers.
It was furthermore described, how the definition of the KC-number in the Darcy-Forchheimer resistance term has a direct influence on the magnitude of the erosion and deposition.The present state-of-the-art approach applies a constant value for the entire coastal structure, and this can be improved by accounting for variations of the KC-number in time and space.
The numerical model was also used to provide insight into the hydraulic gradients at the internal rock-sand interface that led to transport of sand.The computations indicate that if the maximum value of the hydraulic gradients, exceeded by 2% of the waves, remains below a factor 3 times the critical hydraulic gradient, then the amount of transport of sand is negligible.It was found that the relaxation of the criterion for the critical hydraulic gradient is also valid in large scale laboratory experiments.
The developed computational framework has shown to be a valuable tool for the design of filters in permeable coastal structures.
the permeable layer can be written as au bu u fωu + = , where f is a real-valued constant and ω is the cyclic frequency of the oscillatory flow.This finally leads to the following momentum equation: Above the permeable layer, f = 0.
Only the solution adjacent to the wall and inside of the filter is of interest in this analysis.It was found that the homogeneous solution to Eq. ( 22) is described as u g z e = ( ) iωt and the vertical variation of g inside the permeable layer reads: By applying a no-slip boundary condition at the bottom, matching conditions at the interface (the top of the permeable layer for g and g z ∂ /∂ ) and enforcing that u u → 0 as z → ∞, it can be shown that This expression reduces to for f = 0.The real part of k 0 is the inverse of the well-known Stokes length for oscillatory boundary layers.The important part is that k k < 0 , consequently the thickness of the boundary layer decreases with an increasing flow resistance due to the presence of a permeable layer.Therefore, the bed shear stress τ ν u z = ∂ /∂ will be larger in the presence of a permeable structure, if the free stream velocity (measured inside the permeable layer) is identical between two cases with and without a permeable layer.The consequence is that friction factors measured in free oscillatory flows are not applicable inside of a permeable layer with a large (relative to unity) representative resistance coefficient f .

Fig. 1 .
Fig. 1.An example of the erosion and deposition pattern in an open filter under irregular wave loading.H m0 is the spectral wave height and η the instantaneous surface elevation.
details on the wave loading and the material properties of the open filter.The experiments were conducted in the Scheldt Flume at Deltares.The flume is 55 m long, 1.0 m wide and 1.2 m deep.The cases S04 and T09 consisted of a number of wave series with increasing wave height and a constant wave steepness s calibration case, while it was D =0.018 n f 50 m in case S10.Consequently, the critical mobility is smaller following Klein Breteler

Fig. 3 .
Fig.3.The variation in the predicted variation in the area of the erosion hole as a function of time for the 5 calibration sets.The data is also shown.Full circles: A e s , based on the mean profile.Empty circles: A e s , based on the mean of the area of the erosion hole for each profile.'Set' refers to the numerical calibration data and 'Data' refers to the experimental data.

Fig. 4 .
Fig. 4. A comparison between the simulated and experimental profiles.The numerical results are shown for the chosen calibration sets.'Set' refers to the numerical calibration data and 'Data' refers to the experimental data.

Fig. 5 .
Fig. 5.A comparison between measured and predicted bed level change for three of the validation cases.C: There are three values of Ψ cr0 in this panel, namely 0.045 (Set MA), 0.055 (Set MB) and 0.065 (Set MC).'Set'/'Model' refers to the numerical calibration data and 'Data' refers to the experimental data.

Fig. 6 .
Fig. 6.Various properties taken from the 3 calibration cases and the 6 validation cases (3 with a 1:4 slope and 3 with a 1:7 slope).The filled symbols are properties based on the mean of the measured profile.The empty symbols are properties based on the individual profiles.A: The area of the erosion hole.B: The area of the accretion.C: The depth of the erosion hole.D: The height of the accretion.'Data' refers to the experimental data.'Model' refers to the numerical results.

Fig. 7 .
Fig. 7.The effect of the sinking of the outer filter on the profile development.A: The outer rock profile of the various filter configurations (hf ) are depicted along with the resulting profile development of the internal rock-sand interface (h).Dashed black line is the initial profile at the beginning T06-4.B: The accumulated change of the profile during T06-4.'Data' refers to the experimental data.

Fig. 8 .
Fig.8.The effect of the change in profile on two quantities.A: The mean velocity (undertow) inside of the filter layer adjacent to the bed.B: The mean bed load transport.The plot only contains numerical results (see Fig.7).

Fig. 9 .
Fig. 9.The effect of the wave period and the wave height on the erosion and deposition patterns.A: The profiles after 3 h of wave exposure; the reference profile after 6 h is added as the dashed line.The yellow arrow indicates the decreasing response of the sand bed with a decreasing wave period.The red arrow indicates the increasing response with increase wave height.B: The area of the erosion hole relative to the reference solution (s =0.040 p

Fig. 11 .Fig. 12 .
Fig. 11.The classes of configurations for the filter layer.Shaded area is the initial sand core.A: The thickness of the filter layer is constant at z = 0 m.Five different slopes between z = −0.1mand z = 0.1 m are analysed.B: Same as in A, but the modification to the slope is only between z = 0 m and z = 0.1 m.C: This configuration only includes more gentle slopes, and the modification to the slope is between z = −0.1 m and z = 0.1 m.

Fig.
Fig.10.As in Fig.9, however, the KC-number is the same for all cases and equals that of the reference case (KC = 7.40).

n f 15 ,=167 n f 15 .
where d is the sand size.The starting point is consequently Ψ =0.063 cr0 for the narrowly graded filter material studied in Section 3.1 where D d / For the laboratory tests with

Fig. 13 .
Fig. 13.The variation in the accumulated bed level change as a function of water level for profile Configuration B. A: Water level at z = 0.00 m.B: Water level at z = 0.05 m.C: Water level at z = 0.10 m.

Fig. 14 .
Fig. 14.The variation in the accumulated bed level change as a function of water level

Fig. 16 .
Fig.16.The hydraulic loading over the fixed and initial profile of the core.A: The crossshore profiles.The black line is the original profile of the core.B. The distribution of i 10% as a function of the cross-shore position.The term 'Initial' refers to the profile in Fig.15A, 'Edge' refers to the profile in Fig.15Band 'Middle' to the profile in Fig.15C.

Fig. 15 .
Fig. 15.The deformation of the sand core as a function of the layout of the filter after 6 h of morphological time at prototype scale.The shape of the outer filters is taken from laboratory experiments conducted at Deltares.Contrary to the experiments, the filter remained fixed during the wave loading.A: Initial profile.B: Deformed profile at the upwave edge of the cobble beach.C: Deformed profile at the middle of the cobble beach.

d f 2
and d D f n f 50 .The former follows the normalisation chosen by Van Gent and Wolters [40].In Fig. 18A and C the normalisation with d f 2 is plotted and it is observed

Fig. 17
Fig. 17.The variation in z D /

.Fig. 18 .
Fig. 17.The variation in z D / s n f 50 and z D / acc n f 50 as a function of i i max / cr 2% .Negative

1
is a scaling factor introduced in this work; see above.

4. 3 . 1 .
Comparison with large scale physical model results Large scale laboratory experiments were conducted in the Delta Flume of Deltares on a gravel type beach.The focus of the experiments was the erosion and deposition patterns of the gravel material itself (D =0.0145 n f 50

4. 3 . 1 . 1 . 15 ./ n f 15 is approximately 25 .
Simple calculations.The grading of the rock material is identical to that in calibration test T001, where it was found that Ψ For the present case, D d With the help of the diagram in Van der Meulen [37] it is found that Ψ =0.058 cr0 , thus with a slope of the sand core of 1:7.5, this gives Ψ =0.047 cr .The corresponding critical filter velocity is u =0.013 cr m/s.The critical filter velocity is transformed into the critical hydraulic gradients as follows:

Fig. 19 .
Fig.19.The two configurations of the experimental setup that were reproduced with the numerical model.A: Initially, no sand in the filter layer.B: With an initial sand-filled gravel layer on top of the sand core.
Fig. 2. The variation in A e s, and A acc with time for test T001 (experimental data).A: The dimensioned values.B: The ratio between A e s , and A acc.

Table 3
Validation cases: dimensions of the rock layers and their material properties.

Table 4
Validation cases: the wave properties.Each validation case consisted of 3 or 4 wave series.
Colours in Fig.17denote the various wave periods and the symbols denote the thickness of the filter layer.It is seen that practically all of the results follow the same curve with merely a few outliers.If an erosion depth of i can be used.