THREE-DIMENSIONAL MODELING OF LONG-WAVE RUNUP : SIMULATION OF TSUNAMI INUNDATION WITH GPU-SPHYSICS

Tsunamis need to be studied more carefully and quantitatively to fully understand their destructive impact on coastal areas. Numerical modeling provides an accurate and useful method to model tsunami inundations on a coastline. However, models must undergo a detailed verification and validation process to be used as an accurate hazard assessment tool. Using standards and procedures given by NOAA, a new code in hydrodynamic modeling called GPU-SPHysics can be verified and validated for use as a tsunami inundation model. GPU-SPHysics is a meshless, Lagrangian code that utilizes the computing power of the Graphics Processing Unit (GPU) to calculate high resolution hydrodynamic simulations using the equations given by Smooth Particle Hydrodynamics (SPH). GPUSPHysics has proven to be an accurate tool in modeling complex tsunami inundations, such as the inundation on a conical island, when tested against extensive laboratory data.


BACKGROUND Motivation
Since the 2004 Sumatra event, tsunamis have received increased attention at all levels of government.Not surprisingly, assessments of the tsunami hazard for coastal communities have been initiated and tsunami-warning systems have been extended or newly developed.For USA coastlines, Dunbar and Weaver (2007) compiled a report on tsunami hazard assessment.From a scientific point of view, this report highlights two major items.First, that numerical modeling of tsunami generation, propagation, and inundation needs to be included in hazard assessments to predict the environmental and social impact of future events more robustly.Second, that geological information needs to be included for a reliable assessment because large earthquakes may have recurrence periods that are larger than the recorded history in different regions, given the heterogeneous settling history of the North American continent.In most areas, tsunamis leave sedimentary deposits that serve as forensic evidence of events; meaning that such sedimentary layers serve as records of past events.The relevance of tsunami hazard assessments, especially in a probabilistic sense, is described in Synolakis et al. (2007).
However, the sedimentary layers not only record the incidence of a tsunami event, they also contain information on the causative tsunami.As not all tsunami deposits are alike, the complex interactions of the tsunami wave with the movable bed and geometry of the bathymetry cause a generation of specific features of deposits that are a function of the number and characteristics of the incident waves as well as their withdrawal.If deciphered, the information on the tsunami characteristics can be inferred from deposits if --most importantly--the mechanics of the transport, erosion and deposition processes are understood.Our motivation is to tackle the complex physics of the sediment transport processes due to tsunamis and contribute in this way to a more meaningful numerical modeling of past events that will make tsunami hazard assessments more robust and more reliable.

Classical Tsunami Modeling
As every problem in fluid mechanics, the dynamics of tsunamis obeys Newton's Second Law of Motion, and can hence be described with the Navier-Stokes equations.Depending on the processes and the scale these processes operate, the Navier-Stokes and Euler equations can be simplified.For example, for tsunami wave propagation without nonlinear effects, such as dispersion and breaking, the Euler equations can be integrated over depth to arrive at nonlinear and linear versions of the Shallow Water equations.With these equations, non-breaking or weakly breaking waves can be modeled, as demonstrated by Titov andSynolakis (1995, 1998), with the MOST code.The Weather Service of the National Oceanic and Atmospheric Administration (NOAA) uses the MOST code for the early tsunami warning system for Pacific Ocean and US coastlines.However, if processes that operate on a smaller physical scale are considered, such as dispersion or breaking, the Navier-Stokes equations need to be used for simplifications.The Boussinesq equations (Mei, 1983) can be directly derived from the 1 Texas A&M University 2 Johns Hopkins University 3 Istituto Nazionale di Geofisica e Vulanologia 4 Universita di Catania Navier-Stokes equations and allow for dispersion and breaking, but small-scale processes such as plunging cannot be considered.For linear tsunami wave propagation, scales in the order of hundreds of meters are necessary.Wave breaking and dispersion scales between one meter and about ten meters are possible; although, for wave plunging or any processes involving interactions with the movable bed, scales on the order of centimeters to decimeters are involved.
For these small-processes, the three-dimensional Navier-Stokes equations are necessary to be solved.At the scale needed for sediment transport, analytical solutions of the Navier-Stokes equations produce time-finite singularities (e.g., Doering and Gibbon, 1995) which shows evidence that turbulence is an important process on this scale due to the mathematically uncontrollable behavior of this kind of singularity.Given the stochastic characteristics of turbulence, this might mean that a stable and meaningful analytical solution of the Navier-Stokes equations is impossible to achieve.Therefore, numerical solutions need to be computed.In this paper, we validate the model GPU-SPHysics with the conical island experiments, as presented in Briggs et al. (1995), to ensure reliable computations of tsunami inundation on a larger scale.This is a necessary and important step before more complex and non-linear processes are modeled and also validated.For the model validation of GPU-SPHysics, we follow NOAA's standards and procedures for the verification and validation of tsunami models (Synolakis, 2007).

MODELING GPU-SPHysics
SPHysics uses Smooth Particle Hydrodynamics (SPH) to solve equations.SPHysics is freely available and has been used for a variety of scientific problems.SPHysics solves the weakly compressible Navier-Stokes equations on a Lagrangian mesh using kernels as originally proposed by Monaghan (1982Monaghan ( , 1992Monaghan ( , 2005)), Gingold and Monaghan (1982), Benz (1990), and Lui and Liu (2003).Applications of SPHysics for free surface flows, such as waves (Dalrymple and Rogers, 2006;Gomez-Gesteira et al., 2010), range from wave impact on structures (Gomez-Gesteira and Dalrymple, 2003;Crespo et al., 2007;Rogers et al., 2010), the dam-break problem (Crespo et al. 2008), and tsunamis (Crespo et al., 2008;Rogers and Dalrymple, 2008).The SPH method is inherently parallel.A new side branch of the SPHysics development uses the advantages of the parallel properties of the computational problem and the recently emerged computational resources of graphic cards.In concert with these new capabilities, GPU-SPHysics uses C++ and CUDA to access the resources on NVIDIA graphic cards (Herault et al., 2010).The GPU-SPHysics code allows us to run numerical experiments with about four million particles which correspond to four million grid points in a fixed-grid sense.This resolution is obtained using four Tesla C1060 GPUs that provide a total of 960 processors and a speed of four Teraflops, a trillion floating-point calculations performed per second.

Experiments
On December 12, 1992, a 7.5 magnitude earthquake off of Flores Island, Indonesia, caused a tsunami that was responsible for over 750 deaths on the nearby Babi Island (Briggs, et al 1995).The catastrophe caused devastating effects to the island and its inhabitants.One of the most unique aspects of the tsunami at Babi Island was the refraction of the leading wave of the tsunami around the island and its convergence on the backside of the island.The convergence caused an extremely high runup on the back side of the island due to the combination of two waves with equal amplitudes that were in phase.Motivated by the event, the Coastal Engineering Research Center in Vicksburg, Mississippi performed large-scale laboratory experiments in a 30-m wide, 25-m long, and 0.60-m deep wave basin (Briggs, 1995;Synolakis, et al. 2008).The dimensions for the island are 7.2-m toe diameter, 2.2-m crest diameter, 0.625-m height, and a 1:4 slope (Fig. 1a).In the experiments, a piston wavemaker was used to generate solitary waves of different steepness.

Simulations
The geometry of the conical island experiment was entered in GPU-SPHysics.The wavemaker implemented in GPU-SPHysics is the one proposed in Goring (1978), which is slightly different than the wavemaker in the original experiments.Even though Goring's wavemaker produces a slightly larger dispersive tail than the one developed by Synolakis (1986Synolakis ( , 1987)), the characteristics of the main wave remain largely similar.The experiments were done with the maximum number of particles technically possible on the given hardware.
The simulations were done with normalized wave heights of 0.045-m, 0.091-m, and 0.181-m, which are normalized with the water depth of 0.32-m.Wave gauges 6, 9, and 22 were used to compare the water elevation (Fig. 1b).Gauges 9 and 22 are located in vicinity of the still-water shoreline; whereas gauge 6 is located above the toe of the conical island.Furthermore, the runup is also in the comparison.Synolakis et al. (2008) demanded that mass is conserved, which may be a problem for fixed-mesh methods, such as finite differences; however, the SPH method is a Lagrangian method and since no particle can leave the domain during the simulations, mass is conserved.

Data Analysis
For SPH simulation results, water elevations and runup can be directly computed from the location of the particles.For the timeseries at the gauges, a virtual cylinder of radius r c at the locations of the gauges was created.The second step is to find the particle within the radius that has the largest distance from zero height.The third step is draw a virtual sphere around the particle that is farthest away from zero with radius r c .The water elevation is computed by averaging over all particles within the sphere.This is done for all time steps and is used to find the still water surface.
To find the runup on the conical island, virtual circles are drawn from the center of the conical island with increasing radius Δr and with segments of size Δs.If a particle is recorded at a certain position on the slope, the position is flagged and cannot be overwritten by a different position.Circles are created until all segments are filled.The resolution of determining the runup along the conical island is a function of Δr and Δs, which are on the order of two particle diameters.

Results
The comparison between measured and simulated data is excellent.Figure 2 depicts the comparisons at the gauge location 6 (Fig. 2a), 9 (Fig. 2b), and 22 (Fig 2c).The simulated results agree strongly with the measured data for the phase and amplitude of the leading wave.For all waves, the time series at gauge 6, located at the toe of the conical island, exhibits the best match of all the comparisons.At gauge 9, the amplitudes match well.Notably, for H=0.181-m the gauge falls dry between 17-s and 22-s, and the model simulates the dry period of this part of the beach correctly.Also the comparison at the backside of the conical is excellent, and the model is able to reproduce the amplitude and phase accurately.
Comparing the runup around the conical island between experiments and simulations also reveals an excellent comparison.The dots in Fig. 3 indicate the measurements; the solid black line is the simulated runup.The colored intervals around the simulated runup correspond to the 5, 10 and 15% error.For all three wave heights, the comparison on front half of the island is exceptional.For a wave height of H=0.045-m, the data is within the 15 % of the simulated data (Fig. 3a).The runup on the front side is about 0.12-m, and on the backside is 0.08-m.On the backside of the island, the model was able to reproduce the runup distribution and absolute runup excellently.Only a few points are outside of the 15% error range.For H=0.091-m, the model reproduced the runup distribution precisely.On the backside of the island, the runup exceeds, with 0.25-m, the runup on the front side (Fig. 3b) where the maximum runup is 0.22-m.However, for the wave with H=0.181-m, the wave was reported to break, and the maximum runup on the front side is 0.55-m and on the backside is 0.38-m.

CONCLUSIONS AND OUTLOOK
GPU-SPHysics solves the three-dimensional weakly compressible Navier-Stokes equations.Accordingly, GPU-SPHysics will be used to study tsunamis in a three-dimensional and nonlinear framework with a movable bed.However, for this endeavor to occur, validation is needed.The validation presented is part of a series of validation tests for GPU-SPHysics and compares measurements with simulations of the conical island experiments by Briggs et al. (1995).The conical island experiment is part of NOAA's standards and procedures for tsunami inundation models (Synolakis et al, 2007(Synolakis et al, , 2008)).
For the conical island experiments, we compared times of water elevation and the maximum runup around the island at three different gauges located at the front and back of the conical island relative to the wavemaker.Three different wave conditions with wave heights of 0.045-m, 0.091-m, and 0.181-m were used.For the time-series comparison, arrival, amplitude and length of the leading wave matched nearly perfectly for all wave conditions.Also, the simulated runup matches well with the measured runup.For a wave height of 0.091-m, the runup on the backside of the island is larger than the incident wave on the front side.
More validation is needed, especially for processes acting on smaller scales like wave breaking, to show that Smooth Particle Hydrodynamics (SPH) is an efficient and robust alternative to more regular numerical methods, such the finite difference or finite volume approaches.GPU-SPHysics makes use of the parallelism setup on recent graphics cards, which, more recently, allow for an eight times increase in double-precision floating point operations speed and can carry up to 480 CUDA cores each.Because of the parallel characteristics of the SPH method and seemingly continuous developments of graphics cards with more GPUs and more memory, GPU-SPHysics can be an appropriate and efficient tool, in terms of the solved equations, for the modeling of small-scale processes within large-scale features, such as tsunamis.

Figure 1 .
Figure 1.(a) A top view of the conical island experiment showing maximum runup caused by the two refracted waves converging on the back of the island.(b) A schematic sketch illustrating the positions of tidal gauges used to measure runup around the island.

Figure 2 .
Figure 2. A wave elevation η vs. time comparison at each gauge (6, 9 and 22) of the simulated results and the measured data.(a) Shows H=0.045-m at each gauge, (b) shows H=0.091-m at each gauge, and (c) shows H=0.181-m at each gauge.

Figure 3 .
Figure 3.The runup as measured around the island for the simulation and the measured data.The solid black line with 3 categories of color representing increments of 5% error illustrates the simulation runup while the dots represent the measurements.A closer view of the runup around the back side of the island is provided to better show the accuracy of the model.(a) shows H=0.045-m initial wave height, (b) shows H=0.091-m initial waveheight, and (c) shows H=0.181-m initial wave height.