APPLICATION OF SMOOTHED PARTICLE HYDRODYNAMICS IN EVALUATING THE PERFORMANCE OF COASTAL RETROFIT STRUCTURES

This study develops an accurate numerical tool for investigating optimal retrofit configurations in order to minimize wave overtopping from a vertical seawall due to extreme climatic events and under changing climate. A weakly compressible smoothed particle hydrodynamics (WCSPH) model is developed to simulate the wave-structure interactions for coastal retrofit structures in front of a vertical seawall. A range of possible physical configurations of coastal retrofits including re-curve wall and submerged breakwater are modelled with the numerical model to understand their performance under different wave and structural conditions. The numerical model is successfully validated against laboratory data collected in 2D wave flume at Warwick Water Laboratory. The findings of numerical modelling are in good agreement with the laboratory data. The results indicate that recurve wall is more effective in mitigating wave overtopping and provides more resilience to coastal flooding in comparison to base-case (plain vertical wall) and submerged breakwater retrofit.


INTRODUCTION
Coastal zones have been progressively developed in recent decades and have very significant socio-economic value to nations.Protecting the coasts from natural hazards and in specific coastal storms and flooding has been always a key area of research in the past decades.Recent climate change studies (e.g., IPCC 2018) show that not only the sea-level will rise in the future, but more frequent extreme climatic events and coastal storms will occur, which could lead into catastrophic flooding and inundation.Hence, coastal defenses will be under increasing pressure to provide protection for coastal communities and mitigate the long-term effects of climate change.Given that in many regions of the world the national infrastructures are aging, retrofitting of existing structures is an engineering approach that can be adopted to enhance coastal resilience (Dong et al., 2018).
The key design parameter for coastal defenses is the mean overtopping discharge or mean overtopping discharge per unit width of the structure q.Overtopping volume is used as one of the critical parameters in assessing the performance of defense structures.Whilst determining overtopping under different wave conditions could be easily achieved in laboratory, it is a challenging task to predict overtopping for real-world structures and assess the performance of existing structures with retrofits.Understanding the performance of defense structures and retrofits under extreme climatic conditions is very challenging and not even feasible in majority of experimental facilities.Hence, there is a need for an accurate wave overtopping prediction tool which could take into account various structural configurations and hydrodynamic conditions.The majority of existing overtopping prediction tools are based on empirical relations from lab-based data which are not very accurate as the overtopping dynamics rarely resemble the well-controlled expressions presented by empirical predictions.Wave activities and their interactions with defense structures create various responses to wave overtopping which are intricate.Understanding the underlying hydrodynamic mechanisms resulted from wave-structure interactions and ability to model these interactions could help to develop accurate simulation tools which are necessary for evaluating the coastal defenses and retrofit structures.This study adopts the capabilities of Lagrangian particle-based smoothed particle hydrodynamics (SPH) to develop a numerical tool for assessing wave overtopping from coastal retrofit structures.The numerical model is validated against the laboratory data collected at Warwick Water Laboratory and reported in Dong et al. (2018).

NUMERICAL MODEL
A two-dimensional Lagrangian particle-based SPH model is developed based on the weakly compressible formulations (WCSPH).The governing equations for free surface flow simulation in two-dimensional Lagrangian formalism are given by Eq.1 and 2.
where u is velocity vector, P is pressure, ρ represent density, t is time, g is gravitational acceleration, ∇ is vector differential operator and μ is the kinematic viscosity.τ (Eq.2) is the Sub-Particle Scale (SPS) stress tensor (Gotoh et al., 2001), which was initially introduced in the context of moving particle semi-implicit method (MPS).This paper adopts Favre-averaging scheme (Dalrymple and Rogers, 2006), a density-weighted time-averaging method, to implement SPS stress tensor (Eq.3).
where CI is a constant = 0.066 (Blinn et al., 2002), kSPS is SPS turbulent kinetic energy, Sij is SPS strain-rate tensor,   is turbulent eddy viscosity,   represents Kronecker delta and Δ denotes for initial particle-particle spacing.
The turbulent eddy viscosity is determined using Eq. 4 based on Smagorinsky model (Smagorinsky, 1963).
The SPH model developed for this study uses general SPH interpolation equation (Yeganeh-Bakhtiari et al., 2017) and adopts B-spline kernel function (Monaghan, 1992) described in Eq. 5.
where W is weighting function,  is the separation distance between the particles, ℎ is the smoothing distance of a particle interacts with its neighboring particles and  = /ℎ.Inside the smoothing distance, an arbitrary particle  interacts only with the particles exist in a radial distance of order 2ℎ.
where a and b denote the reference particle and its neighbors, respectively,  ⃗  =  ⃗  −  ⃗  ,   =   −   , and φ is a small number introduced to keep the denominator non-zero (= is 0.1hab).
where γ is a constant equal to 7,  =  0 2  0  ⁄ and  0 is the reference density,  0 is the speed of sound at the reference density, and c 0 = c(ρ 0 ) = √∂P ∂ρ ⁄ .In order to satisfy Courant-Friedrichs-Lewy condition, very small constant time-step of 4×10 -5 was used for simulations.
Particles are moving within the numerical domain according to XSPH approach (Monaghan, 1994) [Eq.9] to ensure neighboring particles will move with approximately the same velocity: Where the constant ε is 0.5 and ̅  = (  +   ) 2 ⁄ .The last term in Eq. 9 is called XSPH correction, which ensures that the neighboring particles have approximately the same velocity components [Monaghan, 1994].
The numerical model developed in this study use an Euler Predictor-Corrector time marching scheme (Monaghan, 1989).Eq. 10 summarizes time-marching discretization for momentum, density and position, respectively.
The pressure equation is then calculated as   +1/2 = (  +1/2 ) using Eq. 8.The algorithm then corrects values calculated with Eq. 10 by using force at the half time step (∆/2): The values determined with Eq.11 are finally calculated at the end of each time-step using Eq.12.
To minimize computational efforts and errors as well as conserving both linear and angular momentum, a second order time scheme is adopted.
The numerical domain is setup with a piston-type wave-maker and containers behind the seawall to measure the overtopping volume and investigate the overtopping hazard zone within proximity of the seawall (Fig. 1).Crespo et al. (2007) Dynamic boundary condition was adopted for wall boundaries treatment.To ensure consistent density at walls with the neighboring fluid particles, lines of dummy particles were created at walls.No-slip boundary condition was considered at boundaries with velocity of walls set to zero.Pressure for dummy particles was determined using homogeneous Neumann conditions.Ghost particle approach was implemented to overcome problem of particles penetrating through boundaries and cross walls due to repulsive forces.Ghost particles are forced to satisfy continuity, momentum and equation of state, the same as fluid particles; however, they remain fixed in their position after each time-step and XSPH scheme (Eq.9) is not applied to them.For moving boundaries, wave-maker, the ghost particles are moving at the end of each time-step according to externally imposed function set by wave input condition.
The numerical wave-maker (Fig. 1) implemented for this study is based on Liu and Frigaard (2001) method.This study use JONSWAP spectra for irregular wave propagation in the numerical domain.Biesel transfer function is applied to discretize each component of the spectrum.The algorithm converts the time-series of surface elevations into time-series of piston movement with use of Biesel transfer function (Eq.13).
0, = 2sinh 2 (  ) sinh(  ) cosh(  ) +   𝑑 (13) where S0 is piston stroke, H is the wave height, d is water depth and k is the wave number (=2π/L).Eq.13 can then convert into the time-series of piston displacement as: where e(t) is the piston displacement time-series, ω is the angular frequency (=2π/T=2πfi) and N is the number of parts that JONSWAP spectrum is divided to (>50).
The numerical model described above is used to simulate flow hydrodynamics, wave-structure interactions and overtopping from vertical seawall as well as retrofit structures under different wave and structural conditions.The numerical model is calibrated and validated against the laboratory physical modelling tests performed at Warwick Water Laboratory (Dong et al., 2018).The model is developed and validated with the purpose of serving scientist and engineers as an effective and accurate tool to understand the performance of existing coastal defense structures and retrofits under various hydro-climatic conditions.

LABORATORY DATA
The experimental study was undertaken in a wave flume with dimensions of 22.0 (l) × 0.6 (w) × 1.0 (h) m and beach slope of 1:20, located at Warwick Water Laboratory (Fig. 2a).The flume was equipped with a pistontype wave generator with an active absorption system.Experiments were carried out with vertical seawall fixed at a distance of 12.21m from wave paddle and retrofit structures in front of it (Fig. 2b & 2c).For each test case approximately 1000 pseudo-random waves were generated based on the JONSWAP (γ = 1.0) spectrum, at a 1:50 scale.The characteristics of incident waves and the free-surface elevation were determined using six of wave gauges which were placed close to paddle and seawall, respectively (Fig. 2a).The distance between the gauges was determined based on the Least-Square Method described by Mansard and Funke (1980).
The overtopping volume was measured by a system of collection tank and load-cell which was placed behind the vertical seawall (Fig. 2a).The load-cell was setup to measure wave-by-wave overtopping volume.A syphon mechanism was fixed over the container to allow continuous sampling in the duration of the tests.
Table 1 summarizes details of test conditions for the vertical seawall and retrofit structures.A range of dimensionless freeboard Rc/Hm0, were tested under both impulsive and non-impulsive wave attack.Dong et al. (2018) reported that, for each prototype, tests were conducted with the significant wave height ranging from The laboratory data show that individual overtopping volumes follow Weibull distribution, in agreement with the literature (Hughes et al., 2012;Zanuttigh et al., 2013).Figure 3 illustrates individual overtopping volumes for two test conditions; it is evident that good agreement exists between the data and Weibull plots.The Figure also shows that for larger overtopping resulted from more extreme conditions, the Weibull plots is under predicting the individual overtopping volume.Figure 4 shows the mean overtopping discharge from base-case (vertical seawall) and compares the results to existing prediction formulae for impulsive (Eq.15 and 16) and non-impulsive (Eq.17) wave conditions, proposed by EurOtop (2016).
where  c is the crest freeboard of structure, h is the water depth at the toe of structure,  is gravitational acceleration (=9.81 m/s 2 ),  is the mean overtopping discharge per meter structure width and  −1.0 is statistical wave steepness.Dong et al. (2018) experimental data show that for dimensionless freeboard exceeding 2.25, the overtopping discharge reduced to 98% for recurve wall and 88% for reef breakwater in comparison to the base-case.They reported that for lower range of relative freeboard (Rc/Hm0 < 1) the retrofit impacts on wave overtopping reductions from vertical seawall is lower, where recurve wall had 63% and reef breakwater had 59% reduction in overtopping discharge.Figure 6 illustrates comparison of mean overtopping discharge for base-case (seawall), recurve wall and reef breakwater for the non-impulsive wave conditions tested during laboratory investigations.Dong et al. (2018) laboratory data are compared to EurOtop (2016) empirical relations for non-impulsive wave conditions (Fig. 6).
Further analysis of Figure 5 and 6 indicate that recurve wall has the best performance in reducing mean overtopping discharge from vertical seawall (base-case).The experimental data reveals that reef breakwater is less effective in mitigating wave overtopping impacts on the vertical seawall.For the case of reef breakwater and non-impulsive wave conditions, Dong et al. (2018) reported that the mean overtopping discharge has increased for Rc/Hm0 ≈ 2.25 (Fig. 6) in comparison to the base-case observations.Such increases in overtopping discharge can be associated to relative low wave height and water depth which enhance wave tripping onto the foreshore berm which was previously reported by Allsop et al. (2003 and2005) and Dong et al. (2018).

RESULTS & DISCUSSIONS
This section describes the results of numerical modelling investigation by use of the model developed for this study.The vertical seawall and two retrofits described above are modelled by the WCSPH model.The numerical domain was developed according to the dimensions of physical models to preserve dynamic similarity between numerical study and experimental investigations.JONSWAP (γ = 1.0) spectrum, was implemented for the numerical wave-maker.Base-case and retrofits are modelled according to the dimensions of physical modelling tests (Fig. 2b and 2c).Due to computational limitations it was not possible to simulate 1000 waves for the numerical investigation.For each test case a total number of 30 waves were simulated by the WCSPH model described above.Numerical convergence was analyzed by considering 3 different spatial resolutions for particle spacing dx = 0.05, 0.01 and 0.005.The relative error was determined as the maximum difference between the wave heights generated numerically with theoretical relations (Eq.18).
The comparison of the results for three different particle spacing showed that largest error (~9%) occurred for the roughest resolution and best results obtained for the finest resolution (~2%).This study has used finest particle spacing (dx = 0.005) for modelling vertical seawall and retrofit structures.The WCSPH model is validated by comparing the numerical free surface elevation time series, at different locations near the paddle and in proximity of the seawall, to the data recorded from wave gauges during the physical modelling tests.The model validated and calibrated for the vertical seawall and then used for prediction of overtopping from retrofits.Figure 7 shows the free surface elevations for numerical model and compares them to the wave gauge data near the paddle (WG1, WG2 and WG3) and seawall (WG4, WG5, WG6).Numerical parameters were fine-tuned during the validation and calibration process to improve the accuracy of the model.The results of fine-tuned model confirm that the numerical predictions and laboratory measurements are in close agreements.However, some discrepancies exist between numerical and laboratory data which could be attributed to underlying physical mechanisms of turbulent bores and broken waves.
Following model validation for the base-case, two retrofit structures were simulated by WCSPH model.For all simulation scenarios a time-series of overtopping events were determined.The analysis of data showed that numerically simulated overtopping events are captured at similar time as laboratory measurements which indicate the capability of described model in capturing wave processes and wave-structure interactions with good precision.Further analysis were performed to compare the numerical overtopping volume with both experimental data (Dong et al., 2018) and empirical relations (Van der Meer and Bruce, 2014) Figure 8 (a) compares the results of numerical modelling for base-case and recurve wall with laboratory measurements and empirical relations.It is evident that while data are in good agreement with both measurements and empirical relations, the WCSPH model overestimated the overtopping volume.This could be associated with numerical wave breaking and wave-structure interactions processes.The analysis of numerical results shows that for dimensionless freeboard of 2.5 and higher, no overtopping event occurred for the case of plain vertical seawall and recurve wall.Figure 8(b) compares the laboratory measurements and empirical relations for the reef breakwater retrofit in front of the vertical seawall.The results are in good agreement with the laboratory data and empirical relations.However, the WCSPH model overestimated the overtopping volume of reef breakwater retrofit in comparison to physical modelling tests.Comparison of the results obtained from numerical modelling of vertical seawall and the two retrofit structures show that recurve wall has the best performance in attenuating the wave overtopping from the vertical seawall and this indicates that the results of WCSPH model are in overall agreement with the laboratory measurements.

CONCLUSIONS
This paper presents investigations focused on developing a WCSPH numerical tool for understanding the performance of vertical seawall and retrofitting structures under different hydro-climatic conditions in order to quantify the overtopping volume due to the impacts of wave actions.The numerical model developed within this study is validated and calibrated against the physical modelling tests reported by Dong et al. (2018).
This study shows that WCSPH method is an effective tool for simulating free surface flow, especially for nearshore wave-structure interactions, where large deformations exist.Numerical treatment approaches such as artificial viscosity are implemented in this study to overcome problems associated with truncated kernel at boundaries and free surfaces and to minimize the computational error.The accuracy of overtopping simulation and prediction by the WCSPH model has direct relation to free surface profiles, as unrealistic particle flux could increase or decrease the overtopping volume.The free surface profiles from WCSPH model developed for this study have successfully compared to the wave gauge measurement from physical modelling tests for the basecase (vertical seawall).The calibrated and validated model is then used for simulating the overtopping from two retrofits in front of the seawall (recurve wall and reef breakwater).The overtopping volume resulted from WCSPH model is in good agreement with physical modelling tests.The numerical model has slightly overestimated the overtopping volume in comparison to the laboratory data.The analyses of numerical results showed good agreement with empirical predictions of overtopping discharge from EurOtop (2016).Furthermore, the numerical simulations showed that, for the tested conditions and configurations, recurve wall is the most efficient configurations for reducing wave overtopping, which is in agreement with Dong et al. (2018) physical modelling data.The results confirm that WCSPH model developed within this study could be used as an accurate tool for understanding performance of coastal retrofits with complex geometrical structures under different wave conditions and structural configurations.

Figure 1 -
Figure 1 -Schematic of computational domain and initial arrangement of particles of experimental setup for the base casevertical wall (b) experimental setup for reef breakwater retrofit (c) experimental setup for recurve wall retrofit

Figure 3 -
Figure 3 -Comparison between individual overtopping volume and Weibull distribution for two sets of laboratory data

Figure 4 -
Figure 4 -Laboratory measurement of mean overtopping discharge on plain vertical wall

Figure 5 -
Figure 5 -Laboratory measurements of mean overtopping discharge on plain vertical wall and two retrofit cases of recurve wall and reef breakwater for impulsive wave conditions

Figure 6 -
Figure 6 -Laboratory measurements of mean overtopping discharge on plain vertical wall and two retrofit cases of recurve wall and reef breakwater for non-impulsive wave conditions between wave gauge data and SPH model surface elevation profiles for the case of vertical seawall of 0.25m height, dseawall = 0.1m, Tp=1.25s and Hmo=0.09m at (a) wave gauge locations from paddle WG1 = 1.511m,WG2=1.965m, and WG3=2.147m(b) wave gauge locations from seawall WG4 = 1.474m,WG5=1.022m and WG6=0.832m between mean overtopping discharges obtained from SPH model with physical modelling tests and empirical relations for vertical seawall and (a) recurve wall (b) Reef breakwater