A numerical study of tsunami wave run-up and impact on coastal cliffs using a CIP-based model

There is a general lack of the understanding of tsunami wave interacting with complex geographies, especially the process of inundation. Numerical simulations are performed to understand the effects of several factors on tsunami wave impact and run-up in the presence of submarine gentle slopes and coastal cliffs, using an in-house code, named a 10 Constrained Interpolation Profile (CIP)-based model in Zhejiang University (CIP-ZJU). The model employs a high-order finite difference method, the CIP method as the flow solver, utilizes a VOF-type method, the Tangent of hyperbola for interface capturing/Slope weighting (THINC/SW) scheme to capture the free surface, and treats the solid boundary by an immersed boundary method. A series of incident waves are arranged to interact with varying coastal geographies. Numerical results are compared with experimental data and good agreement is obtained. The influences of submarine gentle slope, 15 coastal cliff and incident wave height are discussed. It is found that the rule of tsunami amplification factor varying with incident wave is affected by angle of cliff slope, and there is a critical angle about 45. The run-up on a toe-erosion cliff is smaller than that on a normal cliff. The run-up is also related to the length of submarine gentle slope with a critical about 2.292m in the present study. The impact pressure on the cliff is extremely large and concentrated, and the backflow effect is nonnegligible. Results of our work are in high precision and helpful in inversing tsunami source and forecasting disaster. 20


Introduction
Tsunami one of the coastal hazards in the world, can be caused by earthquake, volcanic eruption and submarine landslide.
The 2004 tsunami in Southeast Asia is one of the most destructive tsunami events in human history.There were over one hundred thousand victims in 11 countries during the tsunami period (Liu et al., 2005).In the recent Great East Japan Earthquake and Tsunami in 2001, over 24 thousand people were killed or missing and 300 thousand building damages occurred (Mimura et al., 2011).More serious nuclear disaster at the Fukushima Nuclear Power Plants No.1 was caused by the powerful run-up and destructive force of the tsunami wave.All these events and lessons from the previous tsunami wave disasters indicate that the actual tsunami wave run-up and the corresponding destructive forces were underestimated (Dao et al., 2013).
Nat. Hazards Earth Syst.Sci. Discuss., doi:10.5194/nhess-2017-37, 2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 2 February 2017 c Author(s) 2017.CC-BY 3.0 License.existence of cliff can not only influence the impact and run-up of a tsunami wave but also the erosion-deposition.Different layers provide variations in resistance to erosion (Stephenson et al., 2011).Particularly, some coastal cliffs consisting of soft rocks are eroded at the toe (Yasuhara et al., 2002), which makes it easier to be destroyed.It is indispensable to understand the function of coastal cliffs and submerged gentle slopes during the tsunami wave approaches near shore.In this study, tsunami wave run-up and impact on coastal cliffs is simulated using an in-house code, named the CIP-ZJU model.
Considerable attention is paid to the influence of different coastal topographies on the tsunami inundation in near shore areas.
Steep cliffs on the beach and submerged gentle slopes are considered.The submerged gentle slope such as the continental shelf, affects the waveform evolution and wave celerity before tsunami waves reach the shoreline.Tsunami amplification factor, relative wave height, run-up on the cliff and impact pressure will be analyzed in this work.
In this paper, Section 2 describes the governing equations and the numerical methods; Section 3 presents the initial condition and model validation; Dimensionless analysis is then used to examine the effect of front slope length, depth ratio and cliff angles on the run-up, and impact pressure; finally, the conclusion is made in Section 4.

Governing equations
Our model is established in a two-dimension Cartesian coordinate system, based on viscous fluid theory with incompressible hypothesis.The governing equations are continuity equation and Navier-Stokes equations written as follows where u, t, , p,  and fare the velocity, time, fluid density, hydrodynamic pressure, dynamic viscosity and momentum forcing components, respectively.
Multiphase flow theory is employed to solve the problem of solid-liquid-gas interaction.A volume function  m is defined to describe the percentage of each phase in a mesh 0 where m = 1, 2, 3, indicating liquid, gas, solid respectively, and Physical property,such as the density and viscosity in a mesh can be calculated by: Nat. Hazards Earth Syst.Sci. Discuss., doi:10.5194/nhess-2017-37, 2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 2 February 2017 c Author(s) 2017.CC-BY 3.0 License.

Numerical methods
A fractional step approach is applied to solve the time integration of the governing equations ( 1) and ( 2).The first step is to calculate the advection term, neglecting the diffusion term and pressure term, as equation ( 5) shows. (5) A CIP (Constrained Interpolation Profile) method is employed to solve equation ( 5).The CIP method was first introduced by Takewaki et al. (1985) as an efficient method to solve the hyperbolic partial differential equation.The basic principle of CIP is an interpolation in a grid by a cubic polynomial with the value and differential on the grid node.The advantage of CIP is that it can provide a third order interpolation function in a single grid, which makes it a compact high order scheme.
The second step is to solve the diffusion term by a central difference scheme where u * is the solution of equation ( 5), and u ** is the solution to calculate in this step.
The final step is the coupling of the pressure and velocity by considering equation ( 1): Equation ( 7) is solved by a successive over relaxation (SOR) method.More details can be found in our previous works (Zhao et al., 2016a and2016b).
The free surface is captured by a tangent of hyperbola for interface capturing with slope weighting (THINC/SW) scheme, which is based on the principle of VOF method.The THINC/SW was put forward by Xiao et al. (2011), using a hyperbolic tangent function with adjustable parameter to interpolate.The solid boundary is treated by an immersed boundary method (IBM) (Peskin et al., 1973).
3 Numerical results

Numerical wave tank
In this section, 2D numerical wave tanks including incident waves and geographies are introduced.Simulation cases are divided into two categories according to different parameters and purposes.
The first simulations are completed in 1# tank, as shown in Fig. 1.This wave tank is 10.0 m in length and 1.0 m in height.
Four slopes compose the topography-profile, representing continental slope, continental shelf, beach and cliff, respectively.
The still water depth in front of the topography-profile is fixed at 0.35 m, so that the still-water shoreline is located at the starting point of the beach.This point is regarded as the original point of this tank to determine other positions mentioned in  In this work, considerable attention will be paid to 2# tank.It is necessary to number the simulated cases in 2# tank to avoid confusion, as shown in Table 1.

Model validation
To verify the accuracy of our model, numerical result from one of the cases in 1# tank is compared with available experimental data (Sim et al. 2015).The incident wave height and the cliff slope angle of this case are H =0.055m and  4 = 79 o , respectively.A variable grid is used for the computation, in which the grid points are concentrated near the free surface and the topography.Three non-uniform grids are used to perform a grid refinement test.The grid number and the minimum grid size are shown in Table 2. .Fig. 3 concerns the predicted time series of water elevations at different locations S1~3 and the physical measurements (Sim et al., 2015) are also presented for comparison.Fig. 3 (a) illustrates the comparison results at S1.It can be observed that wave has not reached the topography, so that the waveform has not transformed grossly and is similar to the original waveform.Good general agreement is found for all computations.The relative wave height at S1 is 1.0, which reveals the accuracy of the target incident wave. it is hard to determine the true trace of water surface in such a complex condition.The result of fine-grid is between the result of video data and sensor data, the result of middle-mesh is similar to the video data, and the result of coarse-mesh is 3% smaller than video data.In general, our model shows a good performance in this verifiable example, even when the flow regime is extremely volatile.More verification of our model can been found in Zhao et al. (2014).However, for the coarsegrid has a little underestimation and the fine-grid has low time efficiency, the middle-grid will be adopted to complete the remaining case studies.

The tsunami amplification factor in 1# tank
Fig. 4 describes the results of the tsunami wave amplification factor in 1# tank.The vertical coordinates are H m /H r , in which H m means the local wave height and H r means the reference wave height at a reference location x r .The reference location in present study is x r = -0.87m (same as Sim et al., 2015), and the reference wave height H r is provided by present numerical result.From the results in Fig. 4, it is observed that there exists a critical cliff slope about     When the angle of cliff slope is smaller than the critical angle, the tsunami wave amplification factor increases with the increase of the cliff slope angle.When the slope is steeper than the critical slope, the effect of the cliff slope becomes insignificant.This result is similar to Sim et al. (2015).As for Fig. 4  contrast to the gentle cliff.As the cliff slope becomes steep, the differences between different incident waves become small firstly.The critical cliff slope is about  4 45 o .When the angle of cliff slope is greater than 45 o , the differences start to get bigger again.A possible reason for this overturned phenomenon is that the velocity of water particle in the high wave is faster, which allows the high wave easier to run up on the cliff.So that when the cliff slope is gentle, the water of high wave rushes along the cliff and reaches a rearward area, instead of accumulating in the front of cliff as water of small wave cases does.Then, as the cliff slopes get steep, the so-called rearward area becomes hard to reach for the high wave.This change makes the high wave to accumulate water in the area of these two gauge stations.Moreover, the tsunami amplification factor at these two stations keeps increasing with the increase of the cliff slope angle for a given original incident wave no matter the angle is greater or smaller than 45 o , which is different to the results of other four stations mentioned before.Hence, the presence of a cliff does amplify the water elevations on the beach.The influence is particularly evident for the high wave.In present study, the largest tsunami amplification factor is 2.86, as shown in Fig. 4 (f).The predicted run-up is normalized to the incident wave height, H.It can be seen in Fig. 7 that there exists a critical length of submarine gentle slope about L=2.292m.When L < 2.292 m, with a given cliff, the relative wave run-up increases as the length of gentle slope increases.When L goes over the critical length L=2.292m, the wave run-up for the case H = 0.04 m decreases with the increase of the gentle slope length, but results of the case H = 0.06 m still increase with the increase of the gentle slope length.Moreover, the main notable result is that the normal cliff gives an increase result and the toe-eroded cliff gives an opposite result for cases H = 0.05 m.It reveals that the critical value relates to both the incident wave height and the inclination of a cliff.The increase of the incident wave height and the erosion at the cliff toe exaggerates the critical value of submarine gentle slope length.As a whole, the run-up on a normal cliff is higher than that on a toeeroded cliff.The maximum relative run-up on a normal cliff reaches up to 4.2.The possible reason is that the inclination of normal cliff is accordant with the direction of the incident wave while the toe-eroded cliff is contrary.d), (f), (h), (j) and (l).The recorded pressure is normalized to the hydrostatic pressure due to incident wave amplitude, gH.The peak pressure can be divided into two categories: the early peak and the later peak, which are caused by the wave impact and backflow, respectively.It is obvious that the high incident wave height is favorable for the appearance of early peak.Once the early peak appears, it is very larger comparing to the later peak.In the cases of small incident wave height, only later peak is observed at both measurement points.In the action period, the main scope of relative pressure ranges from 1.5 to 2.5, taking no account of the pressure peak.The maximum pressures at five pressure sensors of all cases are also shown in Fig. 9.According to the analysis of Fig. 8, there is no peak when the maximum relative pressure is smaller than 2.5, here our attention will be paid on maximum relative pressure greater than 2.5.The inclination of a cliff affects the appearance of the pressure peak.Under the condition 5 of a toe-eroded cliff, the generation of pressure peak is frequent.As for the normal cliff, the pressure peak is rare, but it can be terribly large once it appears.In addition, the pressure peak is found to be related to the length of submarine gentle slope L. When L is small, it is easier for a small wave to generate an extreme pressure, which may be caused by backflow; and when L becomes large, a high wave has a trend to generate an extreme pressure, which is probably caused by the direct wave impact.10 Combining with the results in Fig. 8, the extreme pressure caused by the direct impact is very large and extraordinarily concentrated, which may destroy the structures and human beings near the coast.While the backflow also produces extreme pressure with a widely affecting area.The worst thing is that it may carry plenty of floats from damaged construction and vegetation, causing secondary damage.

Conclusions
In this study, tsunami wave impact and run-up in the presence of submarine gentle slopes and a coast cliff are investigated numerically using a CIP-based model.Numerical results are firstly compared with available experimental data and the good agreement revealed the ability of our model to solve the complex flow field, such as wave breaking, water-air mixing and violent impact.The results can be summarized as follows.
(1) The angle of cliff slope has a critical value about 45 o , different characteristics of tsunami amplification factor with different incident wave heights has been found when the angle is greater or smaller than 45 o .
(2) The length of submarine gentle slope influences the tsunami wave arrival time and the run-up, and has a critical value about L = 2.292 m in this study.
(3) When wave transforms in front of cliff, the cases with small incident wave height has a large amplification factor, which means a devastating tsunami may be caused by a moderate source.
(4) It is easier for tsunami waves to run up on a normal cliff than on a toe-eroded one.
(5) There are two opportunities for the appearance of pressure peak during the process of tsunami wave run-up and impact.One is the direct impacting pressure when tsunami waves first hit the coastal cliff, and the other is caused by the backflow from the cliff after run-up with a widely affecting area.
The present study gives time history of tsunami evolution from open sea to coastal area, which is rare in field study.
Several topographies and different incident waves has been considered.Comparing with the SWE result, which may underrate and need amendment, present results can simulate the tsunami in near shore areas more accurately.The present model is helpful for tsunami forecast, dangerous prediction and post-disaster analysis.Furthermore, combining with geology knowledge, the earthquake source magnitude and generation location can be determined.
Fig. 3 (b)  shows the results at S2 x = -0.87 m.This gauge point is in the area of the submarine gentle slope, and shoaling happens when wave propagates here.It can be seen from Fig.3(b) that the wave front face becomes steep and the back face becomes gentle, which means wave asymmetry appears.Results of three grids are in good agreement with experimental data.Fig.3 (c) is the most significant among these three graphs, for the gauge station of this graph is located in the front of the cliff where the flow structure is extremely complex.Wave rushes from the coastline in a shape of water jet.Then, it impacts and runs up on the cliff and falls back to the beach.Large quantity of air is entrained in water when backflow interacts with the incident flow.The velocities of water particles fluctuate violently due to the water-cliff interaction and the drastic water-air mixing.Intense spray of water makes it hard to measure the water elevation with a wave gauge.Hence,Sim et al. (2015) employed three HD Pro c910 web cameras to observe wave transformation Nat.Hazards Earth Syst.Sci.Discuss., doi:10.5194/nhess-2017-37,2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 2 February 2017 c Author(s) 2017.CC-BY 3.0 License.besides the Ultralab sensors.Data of video recordings are also presented in Fig. 3 (c) marked by × from Sim et al. (2015).It can be seen from Fig. 3 (c) that the crest value of video data is 18% smaller than the value of sensor data, which reveals that
(e) and (f), the wave gauges are close to the cliff, when the cliff slope is gentle, close to 22 o , the tsunami amplification factor increases with the decrease of the incident wave height.It is noteworthy that under the condition of steep cliff, the tsunami amplification factor increases with the increase of the incident wave height, in Earth Syst.Sci.Discuss., doi:10.5194/nhess-2017-37,2017 Manuscript under review for journal Nat.Hazards Earth Syst.Sci.Published: 2 February 2017 c Author(s) 2017.CC-BY 3.0 License.
Fig.9Maximum impact pressure at five measuing points.

Table 1
Summary of basic parameters calculated in 2# tank

Table 2
Parameter of three sets of grids (Unit: m)