Saint-Venant and Navier-Stokes Models for Tsunami Simulations

: Into the frame of the French TANDEM project (Tsunamis in the Atlantic and the English ChaNnel: Definition of the Effects through numerical Modelling) Principia has been working on the development and qualification of two in-house CFD software: the 2D EOLE-SV (Saint-Venant) model for simulation of large scale tsunami propagation from the source up to coastal scale and the 3D EOLE-NS (Navier-Stokes) model dedicated to tsunami coastal impact modelling. This paper presents a large range of test cases carried out into the frame of the project and dedicated to the validation of numerical codes in various tsunami wave conditions. The main aspects of phenomena such as wave generation, propagation and coastal impact are investigated on academic situations. A real case simulation is concerned as well, the devastating 2011 Tohoku event which is compared with in-situ data.


Introduction 
This work has been performed in Principia within the framework of the French TANDEM project (Tsunamis in the Atlantic and the English ChaNnel: Definition of Effects through numerical Modelling) which is dedicated to the appraisal of coastal effects due to tsunami waves on the French coastlines. TANDEM aims at drawing lessons from the 2011 catastrophic tsunami in Japan. It allows designing, adapting and checking numerical methods of tsunami hazard assessment to define, as accurately as possible, the tsunami hazard for the French Atlantic coastlines in order to provide guidance for risk assessment on coastal civil and nuclear facilities.
To reach these objectives, qualification of numerical codes has been set up, addressing the various stages of a tsunami event: generation, propagation, run-up and inundation. Tsunami modeling involves a major difficulty due to the various scales from the source to the run-up. For instance, considering an earthquake source, the waves Corresponding author: Richard Marcer, PhD, CFD Technical Manager, richard.marcer@principia.fr have initially a small amplitude (around 1 m) and a very long wave length (around 100 km) making strong model assumptions possible. However, close to the coast, complex physical processes like dispersion, nonlinear breaking waves (see for instance [1]), interactions with coastal structures, require more accurate models. Principia have been working on the validation of two complementary models: a depth-averaged Saint-Venant model developed by Principia and Université de Toulon [2] for large scale tsunami propagation simulation, and a fully 3D Navier-Stokes model developed by Principia especially for complex wave breaking and coastal impact problems [3,4].
An overview of the results obtained with both codes aiming at being applicable to tsunami modelling, is presented. The validation process has been done on several academic test cases having experimental data for comparisons, and a real tsunami event:  Breaking of a solitary wave on a constant slope [5].  Propagation of a solitary wave on a 2D reef [6,7]. down a water volume at the rest) [8,9].  Tsunami generation due to a 3D wedge sliding down a slope [10].  Tsunami impact on an urban area [11].  Tohoku tsunami in 2011.
The numerical results are analyzed in order to bring out the ability of both models to simulate the main physical processes of a tsunami life, from large scale propagation using performing remeshing technic, up to coastal processes including wave breaking, wave/structure interactions and different kinds of sources (such as landslide, falling block).

EOLE-NS (Navier-Stokes)
The EOLE-NS code developed by Principia since 1990 is a multi-phase URANS model which solves the two-phase Navier-Stokes equations on structured curvilinear multi-blocks meshes, possibly moving or deforming. It is based on a pseudo-compressibility technique using a dual time stepping method and a second order finite volume scheme for spatial discretization [3]. The interface motion between the different phases is simulated from an implicit VOF model which avoids any CFL constraint and therefore allows globally large time steps. The interface displacement, which corresponds to the transport of the VOF function, may be ensured by a classical Eulerian equation or by an improved Eulerian-Lagrangian method developed by Principia, especially for complex wave breaking problems [3,4]. A mechanical solver with six degrees of freedom is integrated in the code allowing the modelling of fluid/structure coupling. The multi-phase NS model is particularly well adapted for complex flow simulations.

EOLE-SV (Saint-Venant)
The 2D EOLE-SV code is developed jointly at Université de Toulon and Principia. It contains, among others, the Saint-Venant model based on the very long wave assumption leading to non-dispersive depth averaged equations. The solver is based on a finite volume method and unstructured meshes. It uses the robust Godunov solver with a second order MUSCL scheme in space and a Runge Kutta integration in time.
An adaptive mesh refinement (AMR) method based on blocks cutting is implemented (Fig. 1). It enables mesh refinement only where significant phenomena append, allowing optimization of CPU time computation. Details of the methodology are given in Ref. [12].
The Saint-Venant model (often called NLSW for Non-Linear Shallow Water) is well-adapted to simulate large scale tsunami propagation when horizontal processes are preponderant. But when the tsunami approaches the coast, very complex physics phenomena may occur such as wave breaking and interaction with coastal structures, or strong bathymetry gradients. Then 3D modelling is required to get more accurate results. At the edge of complex 3D modelling and the simple Saint-Venant model, Boussinesq-type models [1] are often a good compromise between accuracy and computational efficiency.

Breaking of a Tsunami Wave on a Constant Slope
The experiment, schematically described in Fig. 2,   has been carried out at the California Institute of Technology [5]. A wave maker allows the generation of a solitary wave propagating at first on a constant channel depth area and then on a constant slope beach. The common parameters of this test case are given in Fig. 2. The available experimental data concern the interface evolution of the plunging wave breaking at different times.
The test case is simulated with the NS code. The wave is initialized with a 1st order Boussinesq solitary wave solution with  the wave amplitude equal to 0.45h 0 . Fig. 3 gives comparisons with experiments of the wave interface at two different instants of the breaking process. Results are on the whole satisfactory as the NS model allows capturing with a good phasing and topology the complex physics of the plunging process characterized by a rapid and large deformation of the interface. It can be noticed that a refinement of the mesh would likely improve the representation of the plunging wave crest sharpening. Fig. 4 displays an example of the NS model capability with a plunging wave simulation of splash-up and bubble entrainment processes.

Solitary Wave Reflecting on a 2D Vertical Reef
This study is based on a set of experiments carried out at the O.H. Hinsdale Wave Research Laboratory [6,7]. They involve the propagation, run-up, overtopping and reflection of high amplitude solitary waves on two-dimensional reefs. Their purposes are on one hand to investigate the processes relative to breaking, bore formation, dispersion and transition from sub-to super-critical flows while, on the other hand, providing data for the validation of near-shore wave models in fringing reef.
The geometry of the experiment is shown in Fig The initial depth at still water (h 0 = 2.5 m) gives a partially submerged reef crest. Results are available for the NS model. Comparisons with the experimental

Saint-Venant and Navier-Stokes Models for Tsunami Simulations 46
data of water level distribution along the flume at dimensionless times are given in Fig. 6. For clarity sake, only snapshots representative of propagation, breaking and overtopping phases have been selected.
The first graph, at t' = 66.53, shows the propagation and shoaling of the solitary wave. The NS model provides a peakier wave with regards to the experiments but the code is on the whole satisfactory with the data. At t' = 70.68 the waves overtop the reef, begin to break and initialize the bore formation. Then the sharp bore moves upstream (t' = 80.53) and reflects on the right wall. The resulting reflecting wave passes above the reef for a second time, forming

47
a second bore which propagates upstream (t' = 109.53). On the whole, the model provides precise front positions and amplitudes of the wave during the different phases of the wave propagation. Fig. 7 compares the computed and recorded surface elevation time series at three specific wave gauges, two upstream of the reef in the sloping area and one downstream the reef on the constant depth zone. The first peak of the recorded data at the upstream WG2 and WG8 gauges shows the initial propagation of the wave (WG2) and the moment just before breaking (WG8). For both gauges the later peaks highlight the dispersive and reflective waves induced by the reef. At WG13, initial and reflected bores are well represented. The hydraulic jump developed at the fore reef produces a train of waves over the increasing water depth and the resulting undulations are intensified as higher harmonics are released. The NS

Saint-Venant and Navier-Stokes Models for Tsunami
Simulations 48 model provides a good description of the early phase of the process and a satisfying prediction of undular bore even though higher frequency oscillations are missing in the results.
Detailed comparisons of this Navier-Stokes simulation with results obtained from different codes (NS, NLSW and dispersive depth integrated models) are available in Ref. [13].

Russel's Wave Generator
This case concerns the generation of a long wave induced by the vertical fall of a rectangular rigid body and its interaction with the underlying water. It is based on the experiment published in Refs. [8,9]. The purpose here is to check the accuracy of the NS model in the case of strong interactions between a rigid body and the free surface. In terms of engineering relevance, this case involves the main physics of massive cliffs or ice bodies falling into water and is therefore representative of some tsunami generation processes. Fig. 8 represents a sketch of the experiment which was carried out in a 9 m long flume with a water depth D.
At t = 0 s, a 38.2 kg rectangular block placed just above still water level is released. Experiments are repeated for several water depths, D = 0.288 m, D = 0.210 m and D = 0.116 m. In each case, the block vertical position and the free surface elevation are measured as a function of time, for the latter at a wave gage located at 1.2 m far from the leftward extremity of the flume. Values of H and B which are some characteristic lengths associated with the plunging wave when the wave hits the block (Fig. 9), are also estimated from videos post processing.
To solve fluid/solid interactions, the NS model calculates the resulting pressure force on the block and deduces the free motion from Newton's law. Fig. 10 presents a snapshot of the simulation highlighting the formation of the wave due to the block falling down. Table 1 summarizes the wave height monitored at a distance of 1.2 m from the leftward extremity of the  flume (see Fig. 8), for various initial water depths. On the whole numerical estimations show a good accuracy except for the smaller water depth for which the amplitude is slightly underestimated. Table 2 lists values of H and B (defined in Fig. 9) for D = 0.210 m and confirms the ability of the model to correctly reproduce the complex shape of the wave interface as the block is moving down. Fig. 11 displays a comparison of the free fall solid motion, between NS results and experiments. The curve is to be read from right to left with increasing time. Z/D  1 and Z/D  0 correspond respectively to the extreme positions of the block during the process: release and bottom approach. Satisfactory results are obtained for the fluid/solid coupling model with especially an accurate maximum velocity brought out at the same depth than experiments. Increasing velocity just after the release and decreasing velocity as the block approaches the bottom are correctly reproduced as well.

Tsunami Generation due to a 3D Wedge Sliding down a Slope
This test-case performed with the NS code concerns the modelling of a tsunami generated by a submarine landslide which is represented in the experimental device [10] by a wedge sliding down a 1:2 plane beach slope. Dimensions of the test bench are given in Fig. 12. Comparisons are done with measurements from laboratory data published [10] where different initial wedge positions regarding the free surface height were tested: wedge partially emerged and fully immerged. When the wedge is partially emerged, the top is located at  = 0.1 m above the free surface; when it is fully immerged, the top is located at  = -0.025 m below the free surface. Wave elevation is   monitored at different gauges (wave gauges 1 and 2) whereas the wave run-up is extracted from data on three "run-up" gauges which measure the vertical height (positive or negative) reached by the wave around its initial level (= 0). The detailed positions of the gauges are given in Fig. 12.
In the simulation, the wedge is fully immerged at initial time with a depth of 2.5 mm (distance from the free surface to the top of the wedge) and at rest. Then the block is moved downward along the slope according to a recorded law of motion issued from experiments, with a maximum velocity of 0.22 m/s. Fig. 13 shows visualizations of the wave induced by the sinking of the solid. Red and blue colours express respectively the positive and negative elevations of the wave around level 0. Successive wave crests and troughs are visible and the corresponding fluctuations of wave elevations (with maximum amplitudes ≈ [-10 cm, +5 cm]), and run-up are highlighted on the different gauges (Figs. 14 and 15). Wave gauges show that the first trough and crest are correctly reproduced as well as the wavelength of the waves, except a delay on the third crest for the gauge 2. Amplitudes of run-up are on the whole satisfactorily reproduced, especially when considering the small values of the fluctuations to capture in this problem (less than 5 cm).

Seaside Experiment: Impact on a Urban Area
This experimental case has been carried out in the Oregon State University basin [11]. A complex topography was built including a seawall and several buildings inspired by the city Seaside, Oregon. The offshore wave (height 0.2 m and period 10 s) was designed to correspond to the estimated tsunami wave height for "500 years" Cascadia Subduction Zone tsunami. This test case is particularly interesting to bring out the ability of a CFD model to simulate flows on complex topography including many "macro roughnesses". An overview of the experimental set-up and the probes disposition in the city is shown in Fig. 16.
The simulation is carried out with the NS model on a multi-block mesh of 9 million cells illustrated in Fig. 17. The wave generated experimentally by a piston-type wave maker is numerically represented by a 1st order Boussinesq solitary wave. All the boundaries are supposed to be walls. Fig. 18 shows snapshots of the tsunami impact against the buildings of the city. For numerical results red and blue colors express respectively the positive wave propagation velocity (U max ≈ 1 ms/s) and the reflecting wave velocity (U min ≈ -0.5 ms/s) around the rest (in white). Comparisons with measurements are done at two different instants of the flooding process. Very good results can be noticed for the flooding according to the layout of the city, as well as the reflecting wave by the buildings.  On the whole numerical results of elevation are satisfactory both in amplitude and phasing, whatever the position of the probes. Large amplitude peaks are observed at the impact on the shore-line whereas the wave height decreases strongly as it propagates inside the city. Even for the further inland area (B9 probe) for which the run-up amplitude is expected to be the weaker (maximal amplitude of only 5 cm), the numerical model allows capturing the very small wave elevation amplitudes. Fig. 20 displays cross-shore velocity time series on the same positions as previous probes elevations. Velocities are globally well predicted for the peak values and the time profiles. A quick decrease of the flow velocity is brought out near the shore line (B1) whereas inside the city (B6) the current dynamics is longer maintained due to local ducting effects. It can be noticed that the model tends to underestimate the flow dynamics on the furthest inland area (B9) where the wave amplitude is very small (5 cm). This is probably due to a lack of vertical mesh refinement in this very remote area.
Benchmarks and comparisons of this NS simulation with results obtained with different NSLW Boussinesq codes are available in Ref. [14]. In the same test case, different methods for building modelling in the NS code are presented in Ref. [15].

Simulation of the 2011 Tohuku Tsunami Event
The 2011 earthquake of magnitude 9.1 in Tohoku-Oki (Japan) triggered a tsunami among the most devastating. It broke a lot of seawalls and breakwaters of the Japan coastline, which were designed for smaller critical events. The Tohuku-Oki tsunami has been studied by many authors and reproduced with different numerical models such as [16][17][18][19]. Simulation is here carried out with the SV code using the AMR (Adaptive Mesh Refinement) method. The model covers a large area around Japan in order to assess the reliability of the code for long wave's simulation over very large scales.
The 2011 Tohoku-Oki event is exceptionally well documented. Actually, water buoys and seismic sensors close to the epicentre are available along with water level data on buoys in far field and quite close to the coastline. Therefore, it is possible to qualify the ability of the numerical model to account for the large scale propagation from the source as well as typical physical processes occurring in coastal area like shoaling, refraction and diffraction, due to bathymetry/topography gradients.
The tsunami source is taken from Ref. [20] and considers a slip source alone. The slip distribution is divided into 55 sub-faults. Static sea bottom deformation is calculated by a rectangular fault model assuming an elastic half-space [21]. Horizontal displacement effects are also taken into account.
The model covers an area of about 2,500 km × 5,000 km. The bathymetry discretization provides an accuracy in the range of [10 m-600 m], respectively in the coastal and far field areas. The mesh refinement permitted by the AMR method is 10 m in the coastal zone allowing the modelling of local breakwater impacts, and up to 4,800 m in the far field. Bathymetry and mesh discretization are shown in Fig.  21 with "nc" the local mesh size range given by the AMR method and "nb" the bathymetry interpolation accuracy during the mesh adaptation. In fact, a constraint has been introduced in the AMR method to locally adjust the zone refinement. This method enables the concentration of computational resources on areas of interest such as coastal bays. Fig. 22 yields a global view of the tsunami propagation within the first two hours. The first wave hits the Japan coast before 30 min after the earthquake. As the Pacific Ocean depth is large (about 5,000 m) far off the Japan coasts the transoceanic propagation is very quick, about 800 km/h. Two kinds of quantitative wave elevation comparisons have been performed with far field DART buoys and with GPS buoys closer to the coast (Fig. 23). Fig. 24 presents far-field wave elevation results compared to DART buoys data. Numerical results are comparable to in-situ data especially close to the epicentre (buoy 21418) with correct amplitude and phasing. Elevations are slightly underestimated far away from the source (buoys 21413 and 21419), probably due to a lack of mesh refinement in this far field area or dispersive effects. But it must be noticed that these two buoys bring out quite small wave amplitudes.
Numerical results for nearshore GPS buoys are on the whole very satisfactory (Fig. 25), and better than the preliminary results obtained with the same code based on simplified hypotheses, especially about the

Fig. 22 Tsunami propagation-global view.
Note that the free surface color scale is chosen to highlight the large scale propagation but the leading wave is much higher than 1 m when approaching the coast.  source model [19]. Therefore, wave amplitudes and phases are well represented for all buoys, except for the farthest buoy from the source (GPS 806) where a slight underestimation of the wave amplitude is observed. As the amplitude shift appears right at the beginning of the simulation it may highlight a weakness of the tsunami source model used in the simulation.
A zoom of the tsunami impact in the Kamaishi Bay is presented. During the event, the breakwaters located at the entry of the bay partially collapsed (Figs. 26 and  27).
To highlight the impact of the breakwaters regarding the bay protection, two configurations have been studied. The first case considers the breakwaters as keeping their integrity whereas in the second case they are supposed to collapse. In this last case, as it was impossible to find a detailed description of the damaged breakwaters after the event and then the position of the different elements which have toppled and slipped, it is considered that the breakwaters have instantaneously collapsed on their entire length. The geometrical detail of the removed elements considered in the simulation is given in Fig. 28. The bathymetry cut off has been done at the foot of the breakwaters at the level of -23.5 m. Note that even without breakwaters, the shallow depth of the foot will still have a slight effect on the wave propagation (shoaling, diffraction).
The tsunami propagation in Kamaishi Bay during flooding with and without the breakwaters is shown in Fig. 29. In both cases, large amplitudes of the wave are observed (more than 10 m) leading to a significant flooding of the coastal land. However, the breakwaters clearly limit the wave run-up. As the breakwaters have partially been damaged during the wave entry, the reality might be estimated by considering these two scenarios envelope (i.e. considering an averaged solution between intact breakwaters and fully collapsed).

Conclusion
Following the objectives of the TANDEM project Principia have been working on the validation of its in-house 2D SV and 3D NS codes for tsunami simulations. A qualification protocol has been built mainly from validations on well-documented academic and experimental test cases, each one dedicated to analyze specific tsunami physics (propagation, breaking, impact, run-up, …), as well as a real event. This step-by-step validation approach allowed evaluating the performances of the 2D and 3D models regarding the overall complex physics.
The NS model is basically rather adapted for 3D coastal tsunami impact modelling. Into this frame, the same NS code has been used here for the overall academic test cases of tsunami coastal interaction. This code succeeds in representing with a good accuracy most of the complex coastal processes such as wave breaking, hydraulic bore formation on irregular sea bottom and wave/building interactions. Moreover this model allows the reproduction of the wave generation due to landslide or falling rock, from a dynamic fluid/structure coupling. The 2D SV model performs well for large scale tsunami propagation simulation. The study of the Tohuku tsunami shows that this code is also capable to give interesting results at coastal scale, like for instance the influence of breakwaters regarding their protection against the wave and the issued run-up.