NUMERICAL SIMULATION OF CROSS-SHORE PROFILE EVOLUTION USING OPENFOAM

In the present work, an innovative numerical approach was developed coupling two models in order to simulate the wave propagation over a sloping beach and the sediment transport in surf and swash zones. The first model, synthesized on the basis of OpenFOAM (version 2.4.0) is used to describe the hydrodynamic characteristics of the flow and the wave propagation while the second one is applied for the sediment transport and erosion/deposition prediction using the results of the first model. The method above constitutes an iterative procedure which is tested hereby and seems to yield satisfactory numerical results in comparison with experimental data (Dette 1998).


INTRODUCTION
It is widely known that the waves constitute the main culprit behind seabed change every season.Actually, it is the breaking waves on beach slopes that trigger the movement of the sand as the erosion and deposition create the beach profile every season.Erosion is the process of losing sand which can be devastating for coastal zones and due to its common appearance, many soft or hard engineering solutions are proposed to combat this unwanted process.
Therefore, it is obvious that the prediction of the beach profile is very important in order to prevent the consequences of the erosion in coastal areas and take the necessary measures to protect coastal structures and adjacent properties.Nowadays, we have some tools in our hands that can be used for the prediction of beach formation which are called numerical models.Although it is tough to simulate breaking waves due to their highly turbulent and unstable nature, some advanced numerical models can be created and implemented which yield satisfactory results.
An innovative approach has been tried in the present work as the next step of previous works of the same authors (Karagiannis et al., 2016), where two different numerical models have been used.Firstly, a numerical model, synthesized on the basis of OpenFoam, has been developed and implemented for the wave propagation in order to investigate the hydrodynamic characteristics of the flow.Then, a second model, developed in FORTRAN, has been implemented with the aim of the beach formation prediction, using the hydrodynamic results of the first model.
The aforementioned former numerical model was created with the open source toolbox OpenFoam and the additional toolbox waves2Foam (Jacobsen et al. 2012).The RANS (Reynolds averaged Navier Stokes) equations have been solved simultaneously with the transport equations of the turbulence model k-ω SST, which, after extended investigation, was found the most suitable for free surface cases, and the VOF (Volume of Fluid) method ones (Hirt and Nichols, 1981).
The second numerical model has been developed in FORTRAN, estimating the sheet flow sediment transport rates with the Camenen and Larson (2007) transport rate formula, as well as the bed load and suspended load over ripples (Karambas, 2012).Suspended sediment transport rate is incorporated by solving the depth-integrated transport equation for suspended sediment (Karambas, 2012).After that, the conservation equation of the sediment mass is applied for several time steps and the scour is computed.
The method used in this work is iterative.Specifically, the first model is applied resulting in hydrodynamic characteristics of the flow, which are used from the second model to calculate the sediment transport rates and the erosion/accretion as well.The new bathymetry, which is the result of the second model, is used by the first model and the process is repeated until convergence of the results, namely morphological equilibrium of the seabed.
The numerical results are compared satisfactorily to experimental data (Dette 1998) and specifically to the SAFE Test B2 series experiment.

GOVERNING EQUATIONS OF THE HYDRODYNAMIC MODEL (OPENFOAM)
This section provides the governing equations of the hydrodynamic mathematical model, which is synthesized on the basis of OpenFOAM and solves the Navier -Stokes equations numerically using the Finite Volume Method.The InterFoam solver, which was created for multiphase flows, is available in OpenFoam libraries and used in this work as the most suitable for free surface problems.Free surface is tracked by the VOF method, while the k-ω SST turbulence model is used for the simulation of the turbulence effects.

Continuity and RANS equations
The mathematical model solves the continuity equation, which is as follows : in conjuction with the Reynolds Averaged Navier Stokes (RANS) equations, which are the following ones: where U is the velocity, ρ is the density, g is the gravity acceleration, p is the pressure, μ is the dynamic viscosity, and −  ′  ′ is the Reynolds stress tensor which is equal to the following expression : where μ t is turbulent viscosity coefficient, k is the turbulence kinetic energy and   is the Kronecker delta.The last term in equation ( 2) represents the surface tension effect.

Volume of Fluid Equation
Volume of fluid method (Hirt and Nichols, 1981) is provided by OpenFOAM for the "tracking" of the free surface during free surface simulations.Using this method, every free surface computational cell is divided in two parts, one which represents the air volume and the other is equal to water volume.The calculation of the water-air portion in every cell is possible with the help of the scalar quantity γ, with its value fluctuating between 0 and 1.When the cell is full of water, γ is equal to 1 and when the cell is full of air, γ is equal to 0, while it takes intermediate value (between 0 and 1) when the cell contains both water and air.The γ value is calculated with the following equation for every free surface computational cell: where U r is a relative velocity.Using γ, quantities like density or viscosity can be calculated for every free surface cell as follows: ρ= γρ water + (1-γ)ρ air (5) μ= γμ water + (1-γ)μ air (6)

Turbulence Modelling
The transport equations for the k-ω SST model are as follows (Menter F. R., 1993(Menter F. R., -1994)): where k is the turbulent kinetic energy and ω is the dissipation rate.The rest coefficients are given in literature.

Wave generation and absorption
Cnoidal waves are implemented at the left boundary of the model domain in order to investigate the model's hydrodynamic behaviour (comparing the results with Ting & Kirby's,1994 experimental ones), using the additional toolbox waves2Foam (Jacobsen et al. 2012).The equation implemented, describing the above waves at the inlet, is as follows: where H is the wave height, L is the wavelength, T is the wave period and η 2 is the trough elevation.Further Cn is one of the Jacobi elliptic functions and K(m) is the complete elliptic integral of the first kind.Moreover, sponge layers from waves2Foam libraries are implemented at the left and right end of the computational domain in order to avoid wave reflection which would affect the numerical results.
The wave attenuation at the sponge layers is described by the following equation: where the α R is used in the following equation : where φ may be the velocity or the quantity γ from the VOF equation.

GOVERNING EQUATIONS OF THE MORPHODYNAMIC MODEL (FORTRAN)
The mode of sediment movement on the coast is usually divided into bed load, suspended load and sheet flow transport.Different model concepts are being presently used for the prediction of each one, which range from empirical transport formulas to more sophisticated bottom boundary layer models.
In the present work, the bed load transport (q sb ) is estimated with a quasi-steady, semi-empirical formulation, developed by Camenen, and Larson, (2007): where the subscripts w and n correspond, respectively, to the wave direction and the direction normal to the wave direction, s (= ρ s /ρ) is the relative density between sediment (ρ s ) and water (ρ), g the acceleration due to gravity, d 50 the median grain size, a w , a n and b are empirical coefficients (Camenen and Larson 2007), θ cw,m and θ cw the mean and maximum Shields parameters due to wave-current interaction, θ cn the current-related Shields parameter in the direction normal to the wave direction, and θ cr the critical Shields parameter for the inception of transport.The net Shields parameter θ cw,net in eq. 5 is given by: where θ cw,on and θ cw,off are the mean values of the instantaneous Shields parameter over the two half periods T wc and T wt (T w = T wc + T wt , in which T w is the wave period and α pl,b a coefficient for the phaselag effects (Camenen and Larson 2007).The Shields parameter is defined by: with U cw being the wave and current velocity, f cw the friction coefficient taking into account wave and current interaction and the subscript j should be replaced either by onshore or offshore.
Phase-lag effects in the sheet flow layer were included through the coefficient pl a (Camenen and   Larson, 2007) with: where ν is the kinematic viscosity of water, U w,crsf the critical velocity for the inception of sheet flow, U w is the wave orbital velocity amplitude, W s the sediment fall speed and the subscript j should be replaced either by onshore or offshore.The suspended sediment load (q ss ) may be obtained from (Camenen and Larson 2007): where c R is the reference concentration at the bottom, ε the sediment diffusivity, and U cw,net , the net mean current.
The bed reference concentration is written as follows based on the analysis of a large data set on sediment concentration profiles (Camenen and Larson, 2007): where d * is the dimensionless grain size : The sediment diffusivity was related to the energy dissipation from wave breaking according to Karambas and Koutitas (2002).Phase-lag effect in the suspended concentration due to ripples, is also incorporated according to Camenen and Larson (2007).
The nearshore morphological changes are calculated by solving the conservation of sediment transport equation (Leont'yev, 1996): where z b is the local bottom elevation and q x (=q s,x +q b,x ), q y (=q s,y +q b,y ) are the volumetric sediment transport rates in x and y horizontal directions respectively.

Geometry of the model -Mesh generation
The geometry which is implemented for the hydrodynamic model appraisal is based on the wave flume where the experiments of Ting and Kirby (1994, 1995, 1996) were conducted.There is a 40m long, 0.6m wide and 1m deep wave tank.A plywood false bottom was installed in the wave tank to create a uniform slope of 1:35.The still water depth in the constant-depth region of the wave tank was 0.4m in both the experiments.The 2-D numerical model geometry is depicted below (Figure 1).There are 8 gauges for the spilling breaker, one offshore and 7 on the slope and 7 gauges for the plunging breaker, all of them placed on the slope.The start of the x and z axis is at the tow of the slope, as depicted in Figure 1.The mesh was generated with standard OpenFOAM tools, specifically with blockMesh and snappyHexMesh tools.First of all, the basic structured mesh with same-sized computational cells was created with blockMesh and then snappyHexMesh created the impermeable slope refining the appropriate cells.The refinement process has been done with snappyHexMesh as described in the OpenFoam tutorial and have been studied in some works (Gisen D., 2014).The discretisation used is Δx=2cm and Δz=1cm.Appropriate boundary conditions were implemented.

Numerical experimentscomparison with experimental data (Ting & Kirby 1994)
This section provides the governing equations of the hydrodynamic mathematical model, which is synthesized on the basis of OpenFOAM and solves the Navier -Stokes equations numerically using the Finite Volume Method.The InterFoam solver, which was created for multiphase flows, is available in OpenFoam libraries and used in this work as the most suitable for free surface problems.Free surface is tracked by the VOF method, while the k-ω SST turbulence model is used for the simulation of the turbulence effects.Spilling breaker.Two waves were implemented on the same geometry in order for two different breaking processes to be investigated.The first one, the spilling breaker, is created for wave height H=0.125m and wave period T=2.0s.Moreover, 8 gauges have been placed on the slope at specific positions, so that hydrodynamic characteristics can be quantified.
Results for surface elevation, wave height and trough, wave setup, undertow and turbulent kinetic energy are presented below and compared satisfactorily with experimental data.
Figure 2 shows the distributions of wave height, trough and wave setup/setdown for both numerical and experimental results.The mean water level  is measured from the still water level, whereas the maximum and minimum surface elevations are measured from the mean water level.It can be observed that numerical results are very close to experimental ones, especially for the wave setup/setdown and the wave trough.It must be noticed that the wave of the numerical model breaks earlier and it also cannot reach the maximum height of the wave of the experiment.Figure 4 shows the instantaneous surface elevation for both numerical model and experiment at the same two positions with these in Fig. 3.It can be observed that there is a difference between numerical and experimental results after the wave breaking, as expected.Numerical results for the undertow have been obtained and compared with experimental ones, as can be seen from the Fig. 5.It is about the variation of time-mean horizontal velocity with depth and specifically the variation of the non-dimensional term(  / gh) with depth that is depicted in Fig. 5.The results are referred to the last 2 gauges placed on the slope after the wave breaking point.It can be noticed that the model behaves quite satisfactorily in terms of the undertow investigation at both gauges after the wave breaking point.

Figure 5. Variation of time-mean horizontal velocity with depth (undertow)
It can be observed from the Fig. 5 that the model yields better undertow results as the gauge is closer to the shore after the wave breaking point.
The time-mean turbulent kinetic energy in the surf zone is presented in Fig. 6.Specifically the nondimensional term (  / gh) 1/2 is plotted with depth for two gauges placed on the slope after the wave breaking point.It can be derived from Fig. 6 that the model's results are close to the experimental ones.Plunging breaker.The plunging breaker is created for wave height H=0.128m and wave period T=5sec, while the geometry remains the same.Seven gauges have been placed on the slope, where numerical results can be obtained.Results for surface elevation, wave height and trough, wave setup, undertow and turbulent kinetic energy are presented below and compared satisfactorily with experimental data.
Figure 7 shows the distributions of the wave height, trough and wave setup/setdown for both numerical and experimental results as happened for the spilling breaker.Mean water level  and maximum amplitudes measured from the mean water level are presented in Figure 7, where it can be observed that numerical results fit quite well the experimental ones, when it comes to the wave setup/setdown and the wave trough but as far as the wave height is concerned, the wave of the numerical model breaks earlier and it cannot reach the maximum height of the wave of the experiment.Figure 8 shows the phase average water surface elevation for both numerical model and experiment at two different positions, at x=5m before the wave breaking and x=8.5m after the wave breaking respectively.All distances are measured from the tow of the slope.It can be observed that where the wave breaking has not occurred yet (x=5m), the numerical surface elevation fit the experimental one, whereas there is a little difference between the numerical and experimental surface elevation results after the wave breaking point (x=8.5m)due to the earlier numerical wave breaking.Instantaneous surface elevation for both numerical model and experiment are presented in Figure 9, measured at the same two positions, namely x=5m and x=8.5m respectively.The same conclusion can be drawn, that there is a difference between numerical and experimental results after the wave breaking while they fit each other quite well before the wave breaking.It can be concluded from Fig. 10 that the model's undertow results are better when the gauge is closer to the shore after the wave breaking point.
Figure 11 presents the time-mean turbulent kinetic energy in the surf zone.Specifically the nondimensional term (  / gh) 1/2 is depicted with depth at two gauges, which are on the slope very close to the wave breaking point.Numerical turbulent kinetic energy results are very close to the respective experimental ones, so the model behaves satisfactorily in terms of turbulence prediction.Summary.From turbulence kinetic energy and undertow results can be concluded that the model somehow overestimates a little bit the turbulence before the breaking point, so early breaking occurs and higher turbulent kinetic energy and undertow numerical values are observed before the breaking point.Further calibration of the model may enhance the hydrodynamic results.However, it can be concluded that OpenFOAM is capable of simulating satisfactorily the wave propagation on beach slopes and the numerical results can be further used for sediment transport prediction.

Geometry of the model -Mesh generation
The geometry which is implemented for the numerical simulation of sediment transport is based on the wave flume where the SAFE (Dette et al.,1998) experiments were conducted.There is a 324m long, 5m wide and 7m deep wave tank.The beach was formed from a sand with d 50 =300μm and a fall velocity of w f =0.042m/s.The underwater profile from 4m above the flume floor was shaped to the equilibrium profile h=0.12x 2/3 where h is the water depth and x the distance from shore.The beach above the still water level (4m above the flume floor) has a slope of 1:10 in the present work.The case of storm surges was investigated with random waves, so the water level was raised to 5m above the flume floor.The mesh was generated with standard OpenFOAM tools, specifically with blockMesh and snappyHexMesh tools.First of all, the basic structured mesh with same-sized computational cells was created with blockMesh and then snappyHexMesh created the impermeable slope refining the appropriate cells.The refinement process has been done with snappyHexMesh as described in the OpenFoam tutorial and have been studied in some works (Gisen D., 2014).The discretisation used is Δx=1m and Δz=1cm.Appropriate boundary conditions were implemented.

Seabed profile evolutioncomparison with experimental data (Dette et al. 1998, Test B2)
It has been already mentioned that a wave spectrum was applied with H m0 =1.20m and T m =5.5s as the wave characteristics of the SAFE B2 experiment in order to compare the morphodynamic numerical results with the experimental data.The wave is applied at the left boundary of the numerical wave flume using a suitable boundary condition from waves2Foam libraries.It has been already mentioned that the method used in the present work is iterative.Hence, the hydrodynamic model is applied yielding these results that are used as input in the morphodynamic model.Surface elevation, velocities and turbulent kinetic energy are obtained from the OpenFOAM model and implemented in FORTRAN morphodynamic model.The latter yields the new seabed which is compared with the previous one until convergence.This means that the model is run as many times are needed as the last obtained seabed profile almost coincides with that of the previous run.That means that the seabed has been finally formed under the current wave characteristics.
The hydrodynamic results were shown in detail in the previous section, so the morphodynamic results with the seabed evolution are depicted below.Ten runs were needed until the seabed profile acquires its latest form and comparing this with the final profile from Test B2 experiment, it can be concluded that the numerical results fit the experimental ones very well.Figure 13 depicts the bathymetry after the 1 st run in relation to the initial one.It seems that there are distinguished areas of erosion and accretion.The point where the accretion stops and the erosion starts is at around 245m from the start of the flume in the swash zone, very close to the point where the slope of the beach changes to 1:10.In terms of the water depth, the erosion starts at around 4.5m above the floor, just before the starting point of the swash zone (5m).So, the new seabed is used as input in the OpenFOAM model and a new hydrodynamic state is achieved after the 2 nd run.The hydrodynamic results are used as input in the morphodynamic model again and a new bathymetry is obtained.The above process is repeated continuously until every new run yields identical bathymetry with the previous run.
In the next Figure 14, the result from the 3 rd run is displayed and compared to the initial profile and the bathymetry obtained from the 1 st run.It can be clearly seen that the erosion and accretion are developed but the shape of the seabed remains the same with that one after the 1 st run.It means that the new hydrodynamic state continues to shape the morphology of the seabed in the same way as expected in the breaking and swash zones.However, there is small difference between the bathymetry after the 1 st run and the one after the 3 rd run.The above process is repeated and the shape of the seabed is being developed at a relatively slow pace, namely there is no so much difference between consecutive runs.Hence, not all runs are depicted but some of these in order for the profile evolution to be clearly presented.In the next Figure 15, the bathymetry after 7 runs is shown in comparison with the seabed of the 3 rd step and the initial profile.The areas of accretion and erosion remain the same with the accretion to be more intense after 7 runs.Repeating the same process, the equilibrium profile where the seabed remains the same after each run is reached at the 10 th run.The seabed has acquired its final form and the swash zone is dominated by intense erosion while a bar is formed at the breaking point and accretion appears inside the breaking zone in general.In Figure 16, the profile after the 10 th run is depicted in comparison with the profile after the 7 th run and the initial profile.It is clear that the bathymetry has been developed under the dominated hydrodynamic state which defines the areas of sand movement and the way that movement happens.It is the swash zone that it is being eroded first and offers its sand when the wave attacks after its breaking and takes the sand onshore with its movement.The sediments are deposited in the breaking zone forming a bar as it can be seen in Figure 16.It seems that the method which is hereby presented yields satisfactory results, which is obvious in Figure 17, where the numerical results are compared with experimental data (SAFE Test B2, Dette et al. 1998).It seems that the method, hereby presented, yields satisfactory results, which is obvious in Figure 17, where the numerical results are compared with experimental data (SAFE Test B2, Dette et al. 1998).

CONCLUSIONS
 A holistic numerical approach for the prediction of the seabed formation under propagating waves was presented above.It consists of two models, which are used in a repetitive process with the aim of reaching the equilibrium formation of the seabed after consecutive runs of both models.The first one, synthesized on OpenFoam platform, describes the hydrodynamic state of the flow as it arises from the presence of the propagating waves.The second one, written in FORTRAN, is used to describe the changes in bed morphology, given the hydrodynamic state from the 1 st model. The hydrodynamic OpenFOAM model was evaluated comparing its numerical results to experimental data.Surface elevation, undertow and turbulence kinetic energy obtained numerically fit satisfactorily the experimental data (Ting & Kirby 1994). Earlier numerical breaking affects the results in the case of both spilling and plunging breaker but mostly in plunging breaker because the numerical wave heights remain lower than the experimental ones until the end of the domain.Moreover, there is a little deviation between numerical and experimental wave troughs towards the end of the domain in the plunging breaker case.The model needs further calibration especially in terms of turbulence intensity before breaking in order to achieve the actual breaking point and height. From turbulence kinetic energy and undertow results can be concluded that the model behaves very well but overestimates a little bit the turbulence before the breaking point, so early breaking occurs and high turbulent kinetic energy and undertow values are observed before the breaking point.Further calibration of the model may enhance the results significantly. The morphodynamic model yields satisfactory results and describes quite well the bed formation.
Numerical results fit quite well experimental data. A bar is formed in the breaking zone while erosion appears in the swash zone. The numerical approach presented hereby is a new one, hence comparison with many experiments will be the next step, so that it can be better validated.

Figure 1 .
Figure 1.Sketch of the model's geometry

Figure 2 .
Figure 2. Distribution of wave height, trough and setup/setdownFigure3shows the phase average water surface elevation for both numerical model and experiment at two different positions.The left one is referred to the position x=4.5m from the tow of the slope where the wave breaking has not occurred yet and the numerical surface elevation almost coincides with the experimental one.The right one is referred to the position x=8m after the wave breaking point and it can be observed that there is a little difference between the numerical and experimental surface elevation due to the earlier numerical wave breaking.

Figure 3 .
Figure 3. Phase average surface elevation before and after wave breaking

Figure 4 .
Figure 4. Instantaneous surface elevation before and after wave breaking

Figure 8 .
Figure 8. Phase average surface elevation before and after wave breaking

Figure 9 .Figure 10 .
Figure 9. Instantaneous surface elevation before and after wave breaking Figure 10 depicts the comparison between numerical and experimental results for the undertow,where the model behaves satisfactorily as can be observed.Specifically, what is depicted in these figures, is the variation of time-mean horizontal velocity with depth and specifically the variation of the non-dimensional term(  / gh).The results are referred to the last two gauges placed on the slope after the wave breaking.( / gh) (  / gh)

Figure 11 .
Figure 11.Variation of time-mean turbulent kinetic energy with depth

Figure 13 .
Figure 13.Seabed morphology after the 1 st run

Figure 14 .
Figure 14.Seabed morphology after the 3 rd run

Figure 15 .
Figure 15.Seabed morphology after the 7 th run

Figure 16 .
Figure 16.Seabed morphology after the 10 th run