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

There is a general lack of understanding of tsunami wave interaction 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 gentle submarine slopes and coastal cliffs, using an in-house code, a constrained interpolation profile (CIP)-based model. 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 gentle submarine slope, coastal cliff and incident wave height are discussed. It is found that the tsunami amplification factor varying with incident wave is affected by gradient of cliff slope, and the critical value is 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 a gentle submarine slope with a critical value of about 2.292 m in the present model for most cases. The impact pressure on the cliff is extremely large and concentrated, and the backflow effect is non-negligible. Results of our work are highly precise and helpful in inverting tsunami source and forecasting disaster.

Abstract. There is a general lack of understanding of tsunami wave interaction 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 gentle submarine slopes and coastal cliffs, using an in-house code, a constrained interpolation profile (CIP)-based model. 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 gentle submarine slope, coastal cliff and incident wave height are discussed. It is found that the tsunami amplification factor varying with incident wave is affected by gradient of cliff slope, and the critical value is 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 a gentle submarine slope with a critical value of about 2.292 m in the present model for most cases. The impact pressure on the cliff is extremely large and concentrated, and the backflow effect is non-negligible. Results of our work are highly precise and helpful in inverting tsunami source and forecasting disaster.

Introduction
A tsunami is one of the most disastrous coastal hazards in the world and can be caused by earthquakes, volcanic eruptions and submarine landslides. The 2004 tsunami in South-East Asia was one of the most destructive tsunami events in human history. There were over 100 000 victims in 11 countries during the tsunami period (Liu et al., 2005). In the recent Great East Japan Earthquake and resulting tsunami in 2001, over 24 000 people were killed or went missing and 300 000 buildings were damaged (Mimura et al., 2011). A serious nuclear disaster at the Fukushima Daiichi Nuclear Power Plant 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).
Investigation of tsunami wave transformation in the nearshore area is a feasible approach to learning the action mechanism and inverting the tsunami source. Post-disaster studies were mostly done through field observations and numerical simulations. The sediment in the near-shore area is frequently regarded as traces of a tsunami. Monecke et al. (2008) analysed sand sheet deposited by the 2004 tsunami and extended tsunami history 1000 years into the Aceh past, pointing out the recurrence frequency of damage-causing tsunamis. However, due to the strong backflow of tsunamis, sediment can be brought back to sea. Hence, there may be an underestimation of the run-up using sediment information to study palaeotsunamis. Dawson (1994) pointed out that the upper limit of sediment deposition lay well below the upper limit of wave run-up, which was marked by a well-defined zone of stripped vegetation and soil. Goto et al. (2011) found that previous estimates of palaeotsunamis were probably underestimated by considering their newly acquired data on the 2011 Japanese tsunami event.
However, limited by the complexity in field observations, numerical simulation is another effective approach to investigate tsunamis. The computational domain should include a large area in which a tsunami generates, propagates and inundates (Sim and Huang, 2015). Therefore, shallow water equations (SWE) were popularly used due to their high efficiency. However, when there is interaction between wave and complex geography, SWE is unable to capture flow structures in detail and often underestimates the result (Liu et al., 1991). An amendment is needed to improve the result of SWE by using data from physical experiments or field observation.
Since the available time-series data of tsunami waveform and flow field in the near-shore area are scarce, there is a general lack of understanding of tsunami interaction with complex geography. A more accurate method is required to reproduce the process of tsunami evolution in the coastal area. Due to the development of supercomputing technology and a precise numerical algorithm, computational fluid dynamics (CFD) with viscous flow theory and fluid-solid coupling mechanics are capable of dealing with the complex flow problems when geographies exist. The finite difference method is widely used in various CFD models as a flow solver. Hitherto, the accuracy of the finite difference method has been a great challenge. In this paper, we introduce a CFD model based on the constrained interpolation profile (CIP) algorithm. The CIP method was first introduced by Takewaki et al. (1985) as a high-order method to solve the hyperbolic partial differential equation. Tanaka et al. (2000) proposed a new version of the CIP-CSL4 that overcomes the difficulty of conservative property. Hu et al. (2009) simulated strongly non-linear wave-body interactions used a CIP-based Cartesian grid method, and the results were in good agreement with experiment data. Kawasaki and Suzuki (2015) developed a tsunami run-up and inundation model based on the CIP method, and a highly accurate water surface profile was observed by using slip conditions on the wet-dry boundary. Fu et al. (2017) simulated the flow past an in-line forced oscillating square cylinder using a CIP-based model. The CIP method can be applied in CFD and shows good performance in other areas. Sonobe et al. (2016) employed the CIP method to simulate sound propagation involving the Doppler effect.
It is a significant research project to deal with the freesurface problem in CFD. Hirt et al. (1981) put forward a mass conservation method named "volume of fluid" (VOF), which is flexible and efficient for treating complicated free boundary configurations. Based on the principle of VOF, several improved methods were developed: PLIC (Youngs, 1982), THINC (Xiao et al., 2005), WLIC (Yokoi, 2007) and tangent of hyperbola for interface capturing with slope weighting (THINC/SW; Xiao et al., 2011). Yokoi et al. (2013) proposed a numerical framework consisting of the CLSVOF method, multi-moment methods and density-scaled CSF model. The framework can capture free-surface flows with complex interface geometries well. More recently, conventional VOF has been widely used by combining it with various additional schemes. Malgarinos et al. (2015) proposed an interface sharpening scheme on the basis of the standard VOF method, which effectively restrained interface numerical diffusion. Gupta et al. (2016) used a coupled VOF and pseudotransient method to solve free-surface flow problems, and the numerical solution compared well with analytical or experimental data. Quiyoom et al. (2017) simulated the process of gas-induced liquid mixing in a shallow vessel and found that the mixing time predicted by EL + VOF was in good agreement with the measurements.
When coastal geographies are included, special handling of solid boundary is required. Peskin (1973) proposed an immersed boundary method (IBM) to treat the blood flow patterns of the human heart, which was later introduced to simulate the interactions between solid objects and incompressible fluid flows (Ha et al., 2014;Lin et al., 2016).
CFD is more convenient and economical than laboratorial experiments and field observation. The most attractive feature is that it can provide time-series data of waveform and flow field, which are helpful for a better understanding of the inundation mechanism of waves in near-shore areas. Markus et al. (2014) introduced a virtual free-surface (VFS) model, which enabled the simulation of fully submerged structures subjected to pure waves and combined wave-current scenarios. Vicinanza et al. (2015) proposed new equations to predict the magnitude of forces exerted by the wave on its front face. The equations were added to five random-wave CFD models and good agreement was obtained when compared with empirical predictions. Oliveira et al. (2017) utilized PFEM to simulate complex solid-fluid interaction and free surface, so that a piston numerical wavemaker was implemented in a numerical wave flumes. Regular long waves were successfully generated in the numerical wave flumes.
When the efficient numerical algorithm is adopted, CFD numerical simulation can be applied to study the tsunami inundation in coastal areas. However, limited by the computational efficiency and numerical dissipation of finite difference method, the quantity and slenderness ratio of computational grids should be moderate. Because of the high proportion of spatial span in horizontal and vertical coordinates in most cases, reasonable abnormal model scale in horizontal and vertical is necessary, similarly to laboratory experiments.
The purpose of the present work is to understand the characteristics of a tsunami, help invert the generation mechanisms and provide reference for tsunami forecasting and post-disaster treatment. In this study, tsunami wave impact and run-up on coastal cliffs are simulated using an in-house code, a CIP-based model. Considerable attention is paid to the influence of different coastal topographies, steep cliffs on the beach and submerged gentle slopes. Coastal cliffs are one of the most common coastal landforms, representing approximately 75 % of the world's coastline (Rosser et al., 2005), such as the coast of Banda Aceh in Indonesia and the steep slope at San Martin (Baptista et al., 1993). The existence of a cliff can influence not only the impact and run-up of a tsunami wave but also the erosion deposition. Different layers provide variations in resistance to erosion (Stephenson and Naylor, 2011). Particularly, some coastal cliffs consisting of soft rocks are eroded at the toe (Yasuhara et al., 2002), which causes them to be more easily destroyed. This is indispensable to understanding the function of coastal cliffs and submerged gentle slopes when the tsunami wave approaches the shore. The submerged gentle slope, such as the continental shelf, affects the waveform evolution and wave celerity before tsunami waves reach the shoreline. The tsunami amplification factor (Satake, 1994), relative wave height, run-up on the cliff and impact pressure will be analysed in this work.
In this paper, Sect. 2 describes the governing equations and the numerical methods, Sect. 3 provides the initial condition and numerical wavemaker and Sect. 4 presents 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 discussion and conclusion are in Sects. 5 and 6, respectively.
2 Numerical models

Governing equations
Our model is established in a two-dimensional 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: where m = 1, 2, 3, indicating liquid, gas and solid, respectively, and φ 1 + φ 2 + φ 3 = 1.
Physical properties, such as the density and viscosity in a mesh, can be calculated by

The fractional step approach
A fractional step approach is applied to solve the time integration of the governing Eqs. (1) and (2). The first step is to calculate the advection term, neglecting the diffusion term and pressure term, as Eq. (5) shows.
A CIP method is employed to solve Eq. (5). The second step is to solve the diffusion term by a central difference scheme: where u * is the solution of Eq. (5), and u * * is the solution to calculate in this step. The final step is the coupling of the pressure and velocity by considering Eq. (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, b). The free surface is captured by a THINC/SW scheme, which is based on the principles of the VOF method. The solid boundary is treated by an IBM (Peskin, 1973).

CIP method
The basic principle of CIP is that when computing the advection of a variable f , both the transportation equation of f and the transportation equation of its spatial gradient g = ∂f/∂x are used (seen in Fig. 1).
To explain the steps of the CIP method, we take the following 1-D advection equation as an example.
By differentiating Eq. (9) with respect to x, we have the equation of the spatial derivative: where g = ∂f/∂x. For simplicity, we assume a constant advection velocity. Then Eq. (10) has the same form as Eq. (9). For the case of u > 0, we approximate a profile for f n inside the upwind cell [ The spatial gradient of f n can be written as As shown in Fig. 2, the profile at the time step n + 1 is obtained by shifting the profile with −u · t; i.e. the time evolution of the function f and g can be obtained by using the following Lagrangian invariants.
Therefore we call the CIP scheme a semi-Lagrangian method.
The four unknown coefficients in Eq. (11) can be determined by using known quantities f n i , f n i−1 , g n i and g n i−1 : By introducing the spatial gradient of variable f , CIP can provide a third-order interpolation function in a single grid. Compared to traditional upwind schemes, CIP method has not only sub-cell resolution but also compact high-order structure.

THINC/SW scheme
The THINC/SW scheme, first put forward by Xiao et al. (2011), is used for free-surface capturing of incompressible flows. Some test examples have indicated that the scheme has the features we need, such as mass conservation and a lack of oscillation (Ji et al., 2013). The basic idea of THINC/SW is that, for the profile of a volume function ϕ (0 ≤ ϕ ≤ 1), a hyperbolic tangent function with adjustable parameter is used to interpolate inside an upwind cell, which is shown as follows: where α, γ , δ and β are parameters to be specified. α and γ are used to avoid interface smearing and are given by Parameter δ is used to determine the middle point of the hyperbolic tangent function and is solved by Parameter β determines the steepness of the jump in the interpolation function varying from 0 to 1. In traditional THINC, a constant β = 3.5 is usually used that may result in ruffling the interface, which aligns nearly in the direction of the velocity. To avoid this problem, in THINC/SW, β is determined adaptively according to the orientation of the interface. In a 2-D case, β can be written as 3 Simulation set-up

Initial condition
In this section, 2-D 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 Tank 1, as shown in Fig. 3. 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. The still-water depth in front of the topography profile is is used as an analogue of tsunami in the numerical modelling. S1-3 in Fig. 3 are three gauges of water elevation, located at x = −7.67, −0.87 and 0.11 m, respectively. Outcomes of S1-3 are used for the comparison between numerical results and experimental data (Sim, 2017). Besides, six gauges of water elevation are employed to calculate the tsunami amplification factor, fixed at x = 0.0, 0.06, 0.11, 0.13, 0.16 and 0.21 m (not drawn in Fig. 3).
The second simulations utilize Tank 2, similar to Tank 1 except for slight differences, as shown in Fig. 4a. In this tank, the still-water shoreline also lies on the starting point of the beach, which is the original point of this tank.  Fig. 4b. The scales of these two kinds of wave tanks are same to the previous works of Huang et al. (2013), Sim et al. (2015) and Sim (2017).
In this work, considerable attention will be paid to Tank 2. It is necessary to number the simulated cases in Tank 2 to avoid confusion, as shown in Table 1.

Numerical wavemaker
By declaring a velocity of water particle in the left-most grid and assigning it a value from laboratory wave-paddle veloc-ity, a numerical paddle wavemaker is set at the left side of wave tank (Figs. 3 and 4a).
For a solitary wave, the approximate solution of wave profile near the wave paddle can be described as follows (Boussinesq, 1872): where H , h, c and ξ are the amplitude of the solitary wave, still-water depth, wave celerity and wave-paddle trajectory, respectively. The wave-paddle velocity can be calculated as For a long wave, the depth-averaged horizontal velocity of water particle derived from continuity equation is expressed as (Mei, 1983) .
The horizontal water particle velocity adjacent to the paddle is equal to the wave-paddle velocity, which means that when x = ξ in Eq. (21), u 1 = u 2 . Using Eqs. (21), (23) and (24), wave-paddle trajectory can be derived as an implicit expression: The stroke length of a wave paddle can be calculated as In theory, a period of solitary wave is infinite. In the application, it can be approximately defined as follows: The water particle velocity imposed in the left-most grid can be given by In our model, the left-most grid is not moveable as the laboratory wave paddle is, and Eq. (24) should be modified. By some numerical tests, it is finally determined as

Model validation
To verify the accuracy of our model, numerical results from one of the cases in Tank 1 are compared with available experimental data (Sim, 2017). The incident wave height and the cliff slope gradient of this case are H = 0.055 m and θ 4 = 79 • , 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 quantity and the minimum grid size are shown in Table 2. Figure 5 concerns the predicted time series of water elevations at locations S1-3, and the physical measurements (Sim, 2017) are also presented for comparison. Figure 5a illustrates the comparison results at S1. It can be observed that wave has not reached the topography, so that the waveform has not transformed excessively and is similar to the original waveform. Good general agreement is found for all computations. The relative wave height at S1 is 1.0 m, which reveals the accuracy of the target incident wave. Figure 5b shows the results at S2, x = −0.87 m. This gauge point is in the area of the gentle submarine slope, and shoaling happens when waves propagate here. It can be seen from Fig. 5b 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. Figure 5c is the most significant among these three graphs because the gauge station of this graph is located in the front of the cliff, where the flow field is extremely complex. Waves rush from the coastline in the form of a water jet. Then, the water jet impacts the cliff, accompanied by high pressure acting on the toe of cliff. Great acceleration is produced by the impact, making the water run up on the cliff. Under the influence of gravity, water finally falls back and a 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. The intense spray of water makes it hard to measure the water elevation with a wave gauge. Hence, Sim (2017) employed three HD Pro c910 web cameras to observe wave transformation in addition to the Ultralab sensors. Data of video recordings from Sim (2017) are also presented in Fig. 5c, marked by ×. It can be seen from Fig. 5c that the crest value of video data is 18 % smaller than the value of sensor data, which reveals that it is hard to determine the true trace of water surface in such a complex condition. The result of the fine grid is between the result of video data and sensor data, the result of the middle mesh is similar to the video data, and the result of the 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 unstable. More verification of our model can been found in Zhao et al. (2014). However, as the coarse grid has a slight underestimation and the fine grid has low time efficiency, the middle grid will be adopted to complete the remaining case studies. Figure 6 describes the results of the tsunami wave amplification factor in Tank 1. The tsunami amplification factor is defined as a ratio of the local tsunami height to the tsunami height at a reference location. The vertical coordinates are H m /H r , in which H m is the local wave height and H r is the  reference wave height at a reference location x r . The reference location in the present study is x r = −0.87 m (same as Sim, 2017), and the reference wave height H r is provided by our numerical model. From the results in Fig. 6, it is observed that there is a critical cliff slope of about θ 4 = 45 • . When cliff slope is gentler than the critical value, the tsunami wave amplification factor increases with the increase of the cliff slope gradient. When the slope is steeper than the critical slope, the effect of the cliff slope gradient becomes insignificant. This result is similar to Sim (2017). As for the influence of incident wave height, it can be found in Fig. 6e and f, of which the wave gauges are close to the cliff. When the cliff slope is gentle, close to θ 4 = 22 • , the tsunami amplification factor increases with the decrease of the incident wave height. As the cliff slope becomes steeper, the effect of different incident waves, negligible at first, becomes important. The critical cliff slope is about θ 4 = 45 • . It is noteworthy that under the condition of steep cliff, the tsunami amplification factor increases with the increase of the incident wave height, in contrast to the gentle cliff. A possible reason for this contrary phenomenon is that the velocity of water particles in the high wave is higher, which allows the high wave to more easily run up on the cliff. Thus, when the cliff slope is gentle, the water of a high wave rushes along the cliff and reaches a rearward area, but water of a small wave accumulates at the front of cliff. Then, as the cliff slopes becomes steeper, the so-called rearward area becomes hard to reach for the high wave. This change allows 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 incident wave, whether or not the cliff is steeper or gentler than 45 • . Hence, the presence of a cliff does amplify the water elevations on the beach. The influence is particularly evident for the high wave. In the present study, the largest tsunami amplification factor is 2.86, as shown in Fig. 6f. It is similar to the result of Sim (2017), which has a value of 2.8.

Time evolution of relative wave elevation in Tank 2
Figure 7 depicts the time series of relative wave elevation in Tank 2 for cases 1-12. Four gentle submarine slope lengths and three incident wave heights are considered. The predicted results at x = 0 m are shown in Fig. 7a, c and e, and those at x = 0.4 m are shown in Fig. 7b, d and f. It can be seen in Fig. 7a, c and e that there are conspicuous distinctions between the incident and the reflected wave. The relative height of incident wave at x = 0 m increases with the decrease of the initial wave height. The reflected wave fluctuates remarkably because of the complex flow pattern, and the crest reflected wave is higher than the incident wave. As for Fig. 7b, d and f, when the wave gauges are close to the cliff, it is hard to distinguish the incident and reflected wave. The superposition of incident and reflected wave makes the crest much higher than the results of Fig. 7a, c and e. The effect of initial wave height and length of gentle submarine slope is not clear in Fig. 7 and will be further studied here. The time-series results of cases 13-24 are similar to cases 1-12, which are not shown here.

Relative wave height in front of the cliff in Tank 2
Wave heights at five gauges in Tank 2, x = 0, 0.1, 0.2, 0.3 and 0.4 m, are shown in Fig. 8. The predicted wave height is normalized to the incident wave height, H . The trend line is also presented as the black lines in Fig. 8. In Fig. 8a and b, the maximum relative wave heights are greater than 2.5, and the gradients of trend lines are 2.22 and 2.16, respectively. In Fig. 8c and d, the maximum relative wave height is 2.4 and the trend lines have gradients of 1.77 and 1.75, which are not as steep as those in Fig. 8a and b. In Fig. 8e and f, the trend lines are gentle, with gradients of 1.46 and 1.13. In summary, in the case of smaller wave height, the development of wave height along with the decrease of distance to the cliff is more obvious, and finally a larger relative water height appears. As for larger waves, the rate of wave height increase is very small, especially when the cliff is toe-eroded. The possible cause of this interesting phenomenon can be explained as follows. According to the result of Fig. 7, the crest of wave elevation is produced by the mixing of incident and reflected wave. Under the condition of large waves, the reflected wave is very strong, which causes the mixing to occupy a wide area on the beach. As a result, energy distribution of large waves is not as concentrated as that of the small wave. The energy concentration helps the small wave to produce a higher rela- tive wave height near the cliff. It reveals a moderate surface wave magnitude, which may cause enormous destruction in near-shore areas.
4.5 Wave run-up on the cliff in Tank 2 Figure 9 displays the wave run-up on the cliff in Tank 2; different lengths of gentle submarine slope are compared. The predicted run-up is normalized to the incident wave height, H . It can be seen in Fig. 9 that there is a critical length of gentle submarine slope of about L = 2.292 m. When L < 2.292 m, with a given cliff, the relative wave runup increases as the length of gentle slope increases. When L goes over the critical length 2.292 m, the wave run-up fluctuates. When L > 2.292 m, for both normal cliff and toe-eroded cliff, run-up of the case H = 0.04 m decreases, but results of the case H = 0.06 m still increase with the increase of the gentle slope length. Moreover, for cases H = 0.05 m, the normal cliff results in an increase and the toe-eroded cliff gives an opposite result. It reveals that the critical value relates to both the incident wave height and the inclination of a cliff. In contrast, the run-up on a normal cliff is higher than on toeeroded cliff. The maximum relative run-up on a normal cliff reaches up to 4.2. The reason for this is that the inclination of normal cliff is accordant with the direction of the incident wave, while the toe-eroded cliff is contrary. Waves can easily climb up on the normal cliff and regurgitate slowly along the cliff. As for the toe-eroded cliff, a wave reflects on the cliff and some of the water can run up on the cliff; finally, under the action of gravity, water falls back earlier.
4.6 Impact pressure analysis in Tank 2 Figure 10 shows the time histories of impact pressure on the cliff at two pressure sensors, P3 and P4. The detailed locations of pressure sensors are shown in Fig. 4b. Figure 10a, c, e, g, i and k give the predicted results at location P3, and the results at location P4 are shown in Fig. 10b, 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 directly impact and backflow, respectively. It is obvious that the high incident wave height is beneficial to the appearance of early peak. Once the early peak appears, it is much larger than the later peak. In the cases of small incident wave height, only the later peak is observed at both measurement points. In the action period, besides the pressure peak, the main scope of relative pressure ranges from 1.0 to 2.5.
The maximum pressure at five pressure sensors of all cases is also shown in Fig. 11. According to the analysis of Fig. 10, there is no peak when the maximum relative pressure is smaller than 2.5, and so our attention will be paid to maximum relative pressure greater than 2.5. The inclination of a  cliff affects the appearance of the pressure peak. Under the condition 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 significantly large once it appears. In addition, the pressure peak is found to be related to the length of gentle submarine slope L. When L is small, it is easier for a small wave to generate an extreme pressure, which may be caused by backflow; when L is large, a high wave tends to generate an extreme pressure, which is probably caused by the direct wave impact. Figure 12 demonstrates the snapshots of the pressure distribution at different time instances. Two typical cases of a normal cliff θ 4 = 80.02 • (case 12) and a toe-eroded cliff θ 4 = 91.91 • (case 24) are shown. The predicted results for the case 12 are shown in Fig. 12a, c and e, and those for case 24 are shown in Fig. 12b, d and f. Firstly, the wave impacts on the cliff at the very onset and extreme pressure distribution appears in a small area at the toe of the cliff, near the loca-  tion of pressure sensor P4 (Figs. 12a and b). Because of the cliff, the wave front is deviated and is deflected to a water jet along the cliff. As the water runs up the cliff, it is slowed down by gravity. Also it seems that the pressure impacting on the cliff has been decreased a little (Figs. 12c and d). Moreover, significant water splashing of small water droplets can be observed, contributed from our numerical model. Finally, the water falls back from the cliff due to gravity. This causes the formation of a backward plunging backflow hitting the underlying water, entrapping air. At this stage, extreme pressure appears again and the scope is large, mainly on the beach in front of the cliff. Due to the cliff inclination angle effect, the detailed water flow characteristic on the toe-eroded cliff is found to be different from that on the normal cliff. When the wave propagates along the gentle slope and hits the cliff, the water front horizontal velocity is changed upward along the cliff, as shown in Figs. 12a and b. As time goes on, the wave run-up on the toe-eroded cliff deviates from the cliff due to gravity (Fig. 12d), which will cause a violent plunging breaker onto the underlying water (Fig. 12f). While for a nor-mal cliff the wave runs up along the cliff, complex backflow is generated and later wave breaking and air entrainment occur between the backflow and the underlying water. A comparison of the left and right columns of Fig. 12 reveals that although the general key flow features such as water breaking, water splashing and air entrainment are not significantly affected by the inclination angle of the cliff, the detailed flow features are different. The above findings are partly consistent with the results of Huang et al. (2013). Combined with the results in Fig. 10, the extreme pressure caused by the direct impact is very high and extraordinarily concentrated, which may destroy the structures and result in human losses near the coast. The backflow also produces extreme pressure with a widely affected area; this may carry large amounts of floating debris from damaged construction and vegetation, causing secondary damage.

Discussion
The tsunami amplification factor is essentially a kind of relative wave height that normalized to the height at a reference location. The interesting result in the present work is as follows. In Tank 1, we analyse the tsunami amplification factor near the steepest cliff and find that it increases with the initial wave height (as Fig. 6e and f shown). However, in Tank 2, when the gauge is close to the normal cliff, the relative wave height decreases as the initial wave height increases (seen in Fig. 8a, c and e). It seems that the results of Tank 1 and Tank 2 are contradictory. One of the possible explanations is the influence of the beach. The most significant difference between Tank 1 and 2 is the length of beach. The effect of beach can be simply summed up as follows. The longer the beach is, the more energy is lost before wave impact. The beach is also an area for the mixing of incident and reflected wave; for a large wave, which requires a long area for mixing when the beach is not long enough, the drastic mixing will occur under the coastal line. The details of the beach effect, including process of mixing and energy dissipation, are a meaningful research subject for future work.
As for the run-up in Tank 2, a critical length of gentle submarine slope is found for some cases. Before the wave gets across the coastal line, gentle submarine slope facilitates wave deformation and energy focus. A proper slope helps a wave to adequately prepare before it touches the cliff. When the slope is too long, as the wave reaches the shore line it may have broken or be on the verge of breaking, which makes energy dissipate ahead of time. However, the optimum length is affected by several factors such as initial wave height and cliff slope; this is why there is no critical value found in some cases. From the present work, it is reasonable to speculate that higher initial waves require longer, gentle submarine slopes to achieve the critical value. This can be connected to the analysis of Fig. 11 that when L is small, a small wave generates extreme pressure, and when L becomes large, a high wave tends to generate extreme pressure. In contrast, a normal cliff also enlarges the critical value compared to toe-eroded cliff. The complicated relationship between these factors needs an even deeper investigation.
The present work is only a start as the understanding of tsunami inundation needs to be expanded upon and quantified.

Conclusions
In this study, tsunami wave impact and run-up in the presence of submarine gentle slopes and a coastal cliff are investigated numerically using a CIP-based model. Numerical results are initially 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 gradient of cliff slope has a critical value about 45 • ; different characteristics of tsunami amplification factor have been found when the angle is greater or smaller than 45 • .
2. The length of gentle submarine slope influences the tsunami wave run-up and has a critical value of about L = 2.292 m in this study for some cases.
3. When wave transforms near the cliff, the cases with small incident wave height have a larger relative wave height, 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 cliff.
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 affected area.
The present study gives the history of tsunami evolution from open sea to coastal area, which is rare in field study. Several topographies and different incident waves have been considered. Comparing with the SWE result, which may underrate and need amendment, the present results can simulate the tsunami in near-shore areas more accurately. The present model is helpful for tsunami forecast, danger prediction and post-disaster analysis. Furthermore, combined with geology knowledge, the earthquake source magnitude and generation location can be determined.
Data availability. No data sets were used in this article.
Competing interests. The authors declare that they have no conflict of interest.