Modelling of tsunami-like wave run-up, breaking and impact on a vertical wall by SPH method

Accurate predictions of wave run-up and rundown are important for coastal impact assessment of relatively long waves such as tsunami or storm waves. Wave runup is, however, a complex process involving nonlinear buildup of the wave front, intensive wave breaking and strong turbulent flow, making the numerical approximation challenging. Recent advanced modelling methodologies could help to overcome these numerical challenges. For a demonstration, we study run-up of non-breaking and breaking solitary waves on a vertical wall using two methods, an enhanced smoothed particle hydrodynamics (SPH) method and the traditional non-breaking nonlinear model Tunami-N2. The Tunami-N2 model fails to capture the evolution of steep waves at the proximity of breaking that was observed in the experiments. Whereas the SPH method successfully simulates the wave propagation, breaking, impact on structure and the reform and breaking processes of wave run-down. The study also indicates that inadequate approximation of the wave breaking could lead to significant under-predictions of wave height and impact pressure on structures. The SPH model shows potential applications for accurate impact assessments of wave run-up on to coastal structures.


Introduction
The recent tsunami generated by the Honshu earthquake in 2011 caused tremendous damages to residential and industrial installations at the affected area. Approximately 15 000 people were killed and more than 300 000 buildings were damaged by the earthquake and tsunami (Mimura et al., 2011). The deep-inland and powerful run-up of the tsunami wave damaged the power cooling units of the Fukushima Daiichi nuclear power plant, triggering one of the world's scariest nuclear disasters (US Geological Survey, 2012). All of these events and the lessons learnt from the 2004 Indian Ocean Tsunami disaster have shown that in spite of studies during the past decades, the run-up threats and destructive force of powerful waves have been underestimated.
Various numerical methods have been developed to predict the tsunami wave run-up. Because of hydrodynamic similarities, researchers often use solitary waves as initial conditions to investigate the tsunami run-up characteristics. One of the most popular methods is the boundary integral method (BIM) for potential flow. Although the BIM can get accurate free surface of non-breaking waves (Kim et al., 1983;Maiti and Sen, 1999), the method is very limited when simulating the wave breaking of steep wave fronts (Dommermuth et al., 1988;Grilli et al., 1997). The Eulerian nonlinear shallow water equation (ENSWE) solvers have been used to simulate long-wave propagation and run-up (Li and Raichlen, 2003;Liu et al., 1995;Titov and Synolakis, 1998;. The inherent numerical diffusion and dispersion, otherwise available in Boussinesq approximation, are often used to mimic the physical ones while a turbulence closure model is used to approximate sub-grid flows. The ENSWE method could simulate accurately the wave propagation and adequately the wave run-up. However, the method itself could not resolve the build-up and breaking of steep wave fronts without special treatments of the free surface. The ENSWE solvers coupled with free surface tracking methods, such as in volume of fluid or level set methods, have been successful in capturing the build-up and breaking of steep wave fronts. However, for very steep surface deformation and intensive breaking, very fine meshes are needed to resolve the breaking free surface and water jets created by the breaking wave fronts. Even then, mass conservation is still a challenge for the surface tracking methods, leading to inaccurate modelling of post-breaking processes. Recently, meshless methods have become popular in modelling of similar fluid dynamics problems. The smooth particle hydrodynamics (SPH) method is one of the most robust, efficient, stable and accurate meshless methods. The SPH method was first introduced independently by Gingold and Monaghan (1977) and Lucy (1977) and has later been modified for modelling wave breaking on beaches in two and three dimensions, green waters and wave structure interactions (Colagrossi and Landrini, 2003;Dalrymple and Rogers, 2006;Gomez-Gesteira et al., 2005). Dao (2010) and  extended the method and simulated the propagation, focusing and breaking of a wave group in deep water using an advanced modelling methodology and a very fine resolution. Very promising results that have not been achieved in past numerical studies were obtained, including the wave front evolution, fine flow structures under the breaking wave, water jets and sprays above the wave crest and the evolution of entrapped air bubbles. These features match very well with experimental observations. The original SPH methods, although satisfying the mass conservation, still have zero order in the kernel approximation which sometimes leads to significant reduction of wave height (Liu and Liu, 2006). Several studies (Bonet and Kulasegaram, 2000;Vidal et al., 2007;Oger et al., 2007) have introduced correction terms on the governing equations to improve the kernel approximation to first order. Xu et al. (2010) and Xu (2013) extended and implemented an enhanced correction method proposed by Liu and Liu (2006). The enhanced SPH model was used to simulate sloshing and solitary wave propagations and impacts on to structures. Obtained results showed significant improvement of momentum conservation. The enhanced SPH method is able to simulate wave propagations over a long distance, breaking and impact pressures at a satisfactory level of accuracy.
In this study, the enhanced SPH method is used to simulate solitary wave run-up and run-down over a sloping beach with breaking and impact on a vertical wall. The simulation results of the SPH model are compared with that of the famous tsunami model, Tunami-N2 (Goto et al., 1997), and verified against the experiments conducted at the US Army Corps of Engineers Waterways Experiment Station, Vicksburg, Mississippi . In the experiment, three types of solitary wave run-up, which include non-breaking, closed to breaking, and breaking wave, were considered. A brief introduction of important features of the numerical models is presented, whereas for more details, readers are referred to previous publications on SPH (Dao, 2010;Dao et al., 2011;Xu et al., 2010;Xu, 2013) and the Tunami-N2 model (Goto et al., 1997;Dao and Tkalich, 2007).

Governing equations
The wave propagation and run-up can be governed by the Navier-Stokes equations which take the form of where, ρ and p are the density and pressure of the fluid. In two-dimensional space, u = (u x , u y ) is the fluid velocity and g = (g x , g y )is the gravitational acceleration.

SPH model
SPH method uses a set of discrete particles to approximate the fluid domain. The governing equations are approximated by the interaction of particles within a support domain of a kernel function. This approximation allows the divergences and gradients of fluid properties are transferred to the kernel function. The formulas take the following forms: where i and j are the particle indices, m is the mass of the particles, x is the position of the particles, and W is the kernel function. In this study, the following fifth-order Wendland kernel function is used: where, α D = 7/(4π h 2 ) and q = x − x / h is the normalized distance from particle j to particle i. The parameter h, often called smoothing length, is a measure of the length scale of the support domain of the kernel function. The relationship of the pressure and density of a particle follows the equation of state (Batchelor, 1967): where c 0 is the numerical sound speed, p o is the reference pressure, ρ o is the reference density and γ = 7 for water. This formulation allows particle density to change within a range controlled by the numerical sound speed, thus the approximated fluid is weakly compressible. In this study, singlephase SPH is used. The numerical sound speed for water is chosen as 20 m s −1 . Sensitive studies on setting parameters and simulations of wave generation, propagation, breaking and impact on structures by SPH have been verified in the previous studies by the authors (Dao, 2010;Dao et al., 2011;Xu et al., 2010;Xu, 2013)

Tunami-N2 model
The Tunami-N2 model used in this paper was originally developed in the Disaster Control Research Center (Tohoku University, Japan) through the Tsunami Inundation Modeling Exchange (TIME) program (Goto et al., 1997). Tsunami-N2 solves a set of nonlinear shallow water approximation of the Navier-Stokes equations. The model uses second-order explicit leap-frog finite difference scheme to discretize the set of nonlinear shallow water equations. Horizontal eddy turbulence is neglected and the bottom friction is computed from the Manning formulation. No specific treatment of breaking surface was introduced in the model.
The Tunami-N2 model has been intensively verified using analytical solutions, laboratory experiments and real cases, and subsequently applied to simulate tsunami propagation and run-up in the Pacific, Atlantic and Indian oceans, focusing on particular regional seas and coastal areas (Shuto et al., 1990;Shuto and Goto, 1988;Yalciner et al., 2004;Dao and Tkalich, 2007;Dao et al., 2008).

Experiment setup
The experiment was conducted at the US Army Corps of Engineers Waterways Experiment Station, Vicksburg, Mississippi . The wave flume consisted of a wave maker located at the left side of the flume, a flat bottom followed by three slopes (1 : 53, 1 : 150 and 1 : 13) and a vertical wall located at the right side (see Fig. 1 In the experiment, three different solitary waves were generated by the movement of the piston paddle. The target wave height, h tar /d, for each case was 0.05 (case A), 0.3 (case B) and 0.7 (case C), respectively.

Numerical setups and calibrations
The same wave flume configuration as in the experiment is used in the SPH simulations. A total of about 1.1 million water particles with a uniform diameter of dx = 0.002 m is used. The fine particle resolution is used in order to resolve small wave heights of these cases, especially case A. Secondly, more particles per wave length would improve the numerical decay of wave height. And thirdly, we are not only interested in the wave height of wave propagation, but also interested in the wave impact pressure on the wall. With more particles introduced, the pressure oscillation in the simulation (as a result of slightly-compressible SPH approach) will be reduced. The wave impact pressure on the wall will also be more accurate.
The solitary waves are generated by moving the piston paddle in the model in the same ways as the experiment did. However, because the time resolution of the piston paddle recorded in the experiment is not at a high-enough resolution to derive an accurate paddle velocity, we use wave maker theory (Hughes, 1993;Khayyer et al., 2008) to improve the time history of paddle position and velocity. We observed that the paddle strokes computed by the theory and generated waves at gage 4 are slightly larger than those measured in the experiment. Therefore, in order to generate comparable solitary waves, the paddle strokes computed by the wave maker theory are scaled such that the resulting wave heights at gage 4 match with the experimental data. The paddle velocities are also scaled down with the same ratio (see Fig. 2). Calibrations of solitary waves generated in the SPH model against measurements at gage 4 are shown in Fig. 3. In cases A and B, the measured wave heights match well with the experiment while case C shows a slight under-shoot. After the wave has passed gage 4, the SPH result in case A shows a clear increase in the water level, which is due to the particle disorder at the free surface, a typical characteristic of SPH. The water level increase is of order of a particle size, which is 0.002 m. This particle disorder also happens in cases B and C but it looks obvious in case A because the scale of wave height in case A is only 5 times of the particle size and is much smaller than that in cases B and C.  In the simulations by Tunami-N2, the numerical domain is similar to the real wave flume, except that the left boundary is truncated at the location of gage 4. The time histories of wave height recorded at gage 4 in the experiments are used as prescribed water elevation boundary condition at the left boundary. The horizontal resolution of dx = 0.002 m, which is equal to the particle size in the SPH simulations, is used in all simulations. Manning coefficient is 0.0025. No special treatment for wave breaking is added in Tunami-N2.

Numerical simulations of case A
In case A the solitary wave propagates and runs-up on the slopes without breaking. In both numerical simulations and the experiment, the wave propagates towards the vertical wall and is reflected by the vertical wall without prior breaking (see Fig. 4). The maximum wave heights of about two times of that at gage 4 (or about 0.1 times of still water depth) are observed in the simulation results. The waves are reflected by the vertical wall and no breaking happens during the impact.
Comparisons of wave height time histories at gages 5-10 extracted from numerical simulations and the experiment are shown in Fig. 5. The wave heights obtained by SPH and Tunami-N2 simulations match well with the experiment records at all gages. The reflected wave is also captured in both numerical simulations (the second peaks shown in Fig. 5). The reflected wave obtained by Tunami-N2 is closer to the experiment result than that obtained by SPH. Due to the particle disorder characteristic, water levels obtained by SPH increase approximately a particle size after the wave Nat. Hazards Earth Syst. Sci., 13, 3457-3467 has passed the gages. Overall, Tunami-N2 achieves a slightly better result than SPH in comparison with the experiment. However, one must note that the incident wave at gage 4 used by Tunami-N2 takes the exact value of the experiment, while that of SPH is generated by numerical wave paddle and is slightly different from the experiment record (see Fig. 3).

Numerical simulation of case B
In case B, the incident wave in the experiment is steep and breaking possibly takes place. Snapshots of the solitary wave approaching and impacting on the vertical wall and pressure distribution obtained from the SPH simulation are shown in Fig. 6. Comparisons of time histories of wave height at gages 5-10 obtained from numerical simulations and experiment are shown in Fig. 7. Incident wave from both SPH and Tunami-N2 simulations agree well with experiment results at gage 5 and gage 6. At gage 7, both numerical incident waves are lower than that in the experiment but the SPH result is closer to the experiment. At gage 8 and gage 9, the SPH result agrees very well with the experiment while the Tunami-N2 wave height decreases drastically. We can see in Fig. 7 that, the wave is very steep when it is passing by gage 7. High nonlinearity evolution of the steep wave front is captured in the SPH but not in the Tunami-N2; which could be attributed to the sharp reduction of wave height in the Tunami-N2 simulation. At gage 10, the wave height recorded in the experiment decreases remarkably while that of the SPH simulation does not. This could be because the wave breaking may occur in the experiment somewhere in between gage 9 and gage 10 but not occur in the SPH simulation. As we could see from the simulation results in Fig. 6, when the wave approaches the wall the wave front is steepening but no breaking is observed (Fig. 6a, b). The wave directly impacts on the wall creating a water jet shooting upwards that can reach more than three times of the still water depth (Fig. 6c). The SPH simulation also captures well the reflected wave (the second peak in the time series). At gage 10, the reflected wave recorded in the experiment is significantly lower than that in the SPH simulation. It could be because the incident wave in the experiment already broke before hitting the wall (as discussed above). The experimental and SPH time histories at gage 9 show that the reflected wave could be steepened and break again. The reflected wave in Tunami-N2 is significantly lower than that in the experiment. Overall, the SPH performs much better than Tunami-N2 in case B where the wave is very steep and breaking possibly takes place.

Numerical simulation of case C
In case C, the incident wave is definitely broken before impacting on the wall. The wave breaking is captured in the SPH simulation as shown in Fig. 8. As the wave propagates over the slope the wave steepness increases. The wave front becomes very steep and a jet is generated in front of the wave and is projected forward (Fig. 8a). The jet curls over and plunges on to the water surface in front, generating a process called wave breaking (Fig. 8b-h). Details of the wave breaking process are out of the scope of this paper. Interested  readers could refer to Dao (2010) for the discussions on wave breaking.
Comparisons of time histories of wave height recorded at gages 5-10 obtained from the numerical simulations and experiment are shown in Fig. 9. At gage 5, both numerical results of incident wave agree well with the experiment. The wave steepens when passing by gages 6 and 7. The SPH simulation captures very well the steepening process while the Tunami-N2 does not. The incident wave heights recorded at these gages from the SPH simulation match well the experimental data, while the Tunami-N2 results show a sudden drop in wave height after gage 5 due to lack of the wave-breaking modelling capability. The reduction of wave height could be a result of large numerical dissipation in the simulation when the gradient of the free surface is very large. Wave steepening sometimes causes Tunami-N2 to crash if the time stepping is not small enough. Spurious oscillations may appear at the rear of the wave when it is too steep as seen at gages 5 and 6. Figure 9 indicates that wave breaking happens somewhere in between gage 7 and gage 8 as sharp drops in wave heights are observed in both SPH and experiment results. Very good agreements are also observed between SPH and experiment results at gage 9 and gage 10, including the reflected wave (the second peak in the time series). Tunami-N2 consistently under-predicts both the incident and reflected waves. Moreover, due to the predicted lower total water depth, the wave travels slower and arrives later in the Tunami-N2 simulation as compared to that in the experiment and SPH simulation. Overall, Tunami-N2 performs very badly as compared with SPH in case C where wave breaking takes place early.
Snapshots of breaking wave impact on the vertical wall simulated by SPH are shown in Fig. 10. The solitary wave breaks at the second slope from the left and continues to propagate forward to the wall on the right. The wave front is composed by water sprays splashing forwards instead of a sharp interface (Fig. 10a). The splash of water is approaching the wall (Fig. 10b). As the whole breaking wave is hitting the wall, water splashes upward and the water elevation at the wall increases quickly (Fig. 10c). The main part of the reformed wave impacts on the wall and creates an upward water jet. Water surface at the wall continues to rise and the water jet shoots up. The water jet could reach up to 2.5 times of the still water depth (Fig. 10d) before falling and impacting on the water surface near the corner (Fig. 10e). Even though the incident wave is almost double that of case B, the ratio of maximum run-up height over the still water depth is lower. This could be attributed to the large amount of wave momentum that was dissipated through the breaking before the wave hits the wall. Figures 8 and 10 show the occurrence of some fragmentation of the free surface after the breaking events. This fragmentation of free surface is mainly due to lack of air resistance and incomplete kernel operation at the free surface in the single-phase SPH simulation. The inclusion of air particles in two-phase simulation would significantly improve the result (Dao, 2010).

Impact pressure on wall
Wave impacts on marine and coastal structures have been studied intensively in the past. Experimental studies have shown that the peak pressures due to wave impact could be more than 10 times of those generated by waves of similar amplitude but without impact, and that the magnitude of peak pressure is primarily determined by the stage of wave breaking prior to wave impact Chan, 1994;Chan et al., 1995). A recent study by Lugni et al. (2006) showed that a flip-through wave impact (without prior breaking) on a vertical wall could generate a pressure of 3-6 times higher than a non-impact wave run-up at the impact area. Impact pressures of a breaking wave on to a vertical wall have been successfully modelled by the SPH program developed by the authors. The SPH was able to accurately reproduce the magnitudes of the impact pressures and shapes of pressure evolutions observed in the Lugni et al. (2006) experiments (Xu, 2013).
In the SPH simulations of the solitary waves in this study, the numerical pressure sensor is located on the wall at the initial still water level. The recorded pressure is normalized to the maximum hydrostatic pressure due to incident wave run-up at the wall, ρgh wmax . The results of impact pressures for the three cases are shown in Fig. 11 together with the hydrostatic pressure estimated from incident wave run-up at the wall, ρgh w (Fig. 11a).
When the wave hits the wall and reflects, the pressure measured in case A shows a gradual increase and then decrease. The peak pressure is approximately ρgh wmax . This is, indeed, mostly contributed by the increase of static pressure which is characterized by the increase and decrease of water level at the wall. The normalized impact pressure before and after the wave has hit the wall is non-zero due to the increase of water level resulting from particle disorder in SPH (as discussed above). In case B and case C, impact pressure is characterized by sharp rises of 2 peaks and significant oscillations. Case B shows a typical flip-through wave impact. The second peak impact pressure in this case is approximately 1.5 ρgh wmax and occurs after the peak hydrostatic pressure. The normalized impact pressure is less than 1 at the time of peak hydrostatic pressure and approximately 0.5 during the rising period of hydrostatic pressure; which indicates the estimation of hydrostatic pressure based on rising water level at the wall and gravitational acceleration might be not adequate. That could be explained as the water jet at the wall is significantly accelerating upwards, resulting in lower total downward acceleration and hence lower static pressure. The second peak is explained by the falling to the water jet. Case C characterizes a different impact as the reformed broken wave and water splash impact the wall. There is no strong upward shooting jet as observed in case B. The drop of estimated hydrostatic pressure after the first peak (Fig. 11) is probably due to the fragmentation of water splashing on the wall, making the estimation of water level inaccurate. The peak impact pressure of case C reaches 5.5 ρgh wmax . The magnitude of the peak impact pressures of case C is close to the observation (3-6 times) by Lugni et al. (2006).

Conclusions
In this study, three cases of solitary wave propagation and run-up were simulated by the SPH and the Eulerian nonlinear shallow water (Tunami-N2) numerical models. Comparisons of the simulated wave height time series with the experiments showed that Tunami-N2 is only suitable for nonbreaking waves, whereas the SPH method showed full capability to simulate the whole process of a wave run-up, including complex wave breaking and tremendous impact pressure.  The study also highlighted that neglecting or using an inadequate numerical scheme for the approximation of breaking waves could lead to significant under-predictions of runup height and impact pressure on coastal structures.
The SPH method outperforms traditional gridded methods in simulating breaking surface waves due to several advantages. Firstly, the method discretizes governing equations in Lagrangian form, leading to several benefits: surface tracking is not necessary; the advection term can be treated exactly by tracking the movement of the particles; and numerical diffusion due to grid discretization can be avoided. Secondly, because the computation elements are based on the particles, the air-water multiphase flow could be modelled easily. Furthermore, the SPH method is based on the interaction between particles, hence the continuum, coalescence and fragmentation of the materials could be treated in a natural way.
The SPH method is, however, much more computationally intensive than the ENSWE due to the Lagrangian approach. The computational demand sometimes prohibits the SPH method from real-time applications. However, the method is very useful in post studies of extreme events or in designing structures that are subject to complex wave impacts. It is also useful in the construction of wave run-up databases in which combinations of different beach slopes and incident waves are simulated offline. Such databases will be used for data-mining or data-driven models that could provide quick and adequate impact assessments of tsunami or storm wave run-up on beaches given the incident waves at the far-field obtained from a real-time ENSWE model.