Numerical study of dam-break induced tsunami-like bore with a hump of different slopes

Numerical simulation of dam-break wave, as an imitation of tsunami hydraulic bore, with a hump of different slopes is performed in this paper using an in-house code, named a Constrained Interpolation Profile (CIP)-based model. The model is built on a Cartesian grid system with the Navier Stokes equations using a CIP method for the flow solver, and employs an immersed boundary method (IBM) for the treatment of solid body boundary. A more accurate interface capturing scheme, the Tangent of hyperbola for interface capturing/Slope weighting (THINC/SW) scheme, is adopted as the interface capturing method. Then, the CIP-based model is applied to simulate the dam break flow problem in a bumpy channel. Considerable attention is paid to the spilling type reflected bore, the following spilling type wave breaking, free surface profiles and water level variations over time. Computations are compared with available experimental data and other numerical results quantitatively and qualitatively. Further investigation is conducted to analyze the influence of variable slopes on the flow features of the tsunami-like bore.


Introduction
Recently, earthquake and subsequent tsunami have been making scenes around the world. The December 2004 Indian Ocean Tsunami, the March 2010 Chile Tsunami and the March 2011 Japan Tsunami caused severe damage and losses of life to numerous coastal communities underscored the need for a better understanding of tsunami-structure interaction (Linton et al., 2013). Recorded videos from the December 2004 Indian Ocean Tsunami showed that, in most areas, the tsunami transformed into a hydraulic bore which propagated onshore at a considerably high velocity and advanced towards the shoreline and continues on as an overland flow impacting onshore infrastructures (Nouri et al., 2010). Based on the existing analogies between a tsunamiinduced bore and a dam-break wave in the literature (Chanson, 2006), during the following decade, a dam-break hydraulic bore has been generated by lots of scholars to simulate a tsunami-induced bore (Yeh, 2007;Nouri et al., 2010;Al-Faesly et al., 2012;Wei et al., 2014).
Surge waves resulted from dam break have also been responsible for numerous losses of life. Along with more dams' emergence in the rivers and the rapid development of the coastal areas of the world over the past years, it is obvi-ous that even more lives will be lost and more property will be damaged as the next tsunami or dam breaking comes. According to Nordin and Charleson (2009), Malaysian coastal dwellings may be under threat from tsunami in the future. To date, even though numbers of previous studies have been carried out in this area, there are still doubts on the understanding of the tsunami bore interaction with structures because of its dynamic and complex behavior, and in design of the infrastructure accordingly (Manawasekara, 2013). As a typical structure model of seadike, the triangle obstacle is used to study the interaction between tsunamilike bore and the safety protection facilities.
The positive and negative dam-break bores in a horizontal rectangular channel with dry bed were investigated by Lauber and Hager (1998a) physically. An experimental approach was taken by Nouri et al. (2010) where a dambreak induced flow was used as an imitation of tsunami hydraulic bore, generated by the fast opening of a gate, impacted various free standing structures of different shapes located downstream of the gate. Al-Faesly et al. (2012) presented the results of a comprehensive experimental program focused on the impact of tsunami hydrodynamic forces on structural models generated by a dam-break hy-draulic bore. Ozmen-Cagatay et al. (2014) carried out experimental works on dam-break flows with obstacles and wet bed.
Numerical models of dam-break flows have received considerable attention in the past few years. Academic attempts regarding the numerical models of dam-break flows can be classified into two major approaches, i.e., meshbased methods and mesh-free methods. Meshless method is a kind of modeling method without any grid,and only uses a particle view throughout the problem domain. Mesh-free methods have been applied to hydrodynamics, such as the Smoothed Particle Hydrodynamics (SPH) method (Monaghan, 2012;Zheng et al., 2014) or the Moving Particle Semi-implicit (MPS) method (Koshizuka and Oka, 1996;Tang et al., 2016). Dam break bore evolution over dry and wet beds was analyzed with a SPH model by Crespo et al. (2008). Shakibaeinia and Jin (2011) developed a new meshfree particle model based on the weakly compressible MPS formulation to model dam break over a mobile bed. However, due to the oscillations in the pressure field and low computational efficiency, the mesh-free methods are limited.
Recently, with the advantage of mass conservation and easy to implement, one of the mesh-based methods called Volume of Fluid (VOF) (Hirt and Nichols, 1981) which solves the Navier-Stokes equations with free surface modeling draws much attention in the field of the hydrodynamics. Some improved Volume of Fluid (VOF)-type schemes have been proposed, such as THINC (tangent of hyperbola for interface capturing) (Xiao et al., 2005), VOF/WLIC (WLIC: weighed line interface calculation) (Yokoi, 2007) and THINC/SW scheme (Xiao et al., 2011). Biscarini et al. (2010) studied free surface flows induced by a dam break based on the solution of the complete set of Reynolds-Averaged Navier-Stokes (RANS) equations coupled with the VOF method. Three-dimensional numerical simulations of the dam-break flow were carried out using VOF by Liu et al. (2013) to study the spread of dam-break flow on the complicated terrain. Marsooli and Wu (2014) developed a three-dimensional VOF-based model to simulate dam-break flow over uneven beds in irregular domains.
Though dam-break induced tsunami-like bore problems have been studied both experimentally and numerically for many years and some developments have been made in the study, the understanding of dam-break hydraulic bore with infrastructures is still insufficient due to its complexity, such as the spilling type wave breaking, the air-water interaction and multiphase problem. Therefore, further researches are still needed. In order to improve the solution of the multiphase problem and capture the free surface robustly, a scheme of high resolution based on the framework of the finite difference method has been proposed. The CIP (Constrained Interpolation Profile) method was first presented by Takewaki et al. (1985) in order to reduce the numerical dif-fusion of low order finite difference method. The CIP-based Cartesian grid model has been developed for violent flow problems (Hu and Kashiwagi, 2009, Zhao and Hu, 2012, Cao and Zhao, 2015. In this paper, a CIP-based model is used for the study of dam-break waves in a dry channel with a hump of different slopes.
As a typical structure model of seadike, the triangular hump has been widely applied in offshore engineering to protect the people and properties behind it. A natural hump always has a gentle slope, while the artificial hump always has a steep one. As we know, it is complex and difficult to change the slope of the model in the experiment. The purpose of this paper is to study the dam-break induced tsunami-like bore with a hump of varying slopes by a CIPbased model numerically. An improved model governed by the Navier-Stokes equations with free surface boundary conditions is presented, in which a more accurate VOF-type scheme, the THINC/SW scheme is adopted for interface capturing and the immersed boundary method is used for the treatment of the solid-liquid interface. The model will be shown to accurately fit the experimental results. In addition, the model captures most features of a dam break with a hump; in particular, it allows analyzing the spilling type negative bore, spilling type wave breaking and air-water interaction associated with the interaction between the dam break flood water and the hump. Moreover, based on the capability of this model mentioned above, a series of cases are designed to investigate the effect of the slope of the hump on the features of the tsunami-like wave.
The paper is structured as follows. In the next section, the CIP-based model is introduced. In Section 3, the initial condition of the simulation is first described. Secondly, the effect of the mesh density is analyzed. Then, the numerical results and experimental data of the water surface profiles and water elevations at different locations are compared. Lastly, numerical calculations are performed for fourteen kinds of different slope cases (up1-7 and down1-7) to investigate the effect of the slope of the hump on the features of a dam-break flood wave. Finally, the article is closed with some general conclusions.

Governing equations
In Consideration of the two-dimensional viscous incompressible flow, the governing equations are mass conservation equation and the Navier-Stokes momentum equations written as follows: where Cartesian tensor notation is used, t, u j , p and x j rep-resent time, velocities, hydrodynamic pressure and spatial coordinates, respectively. f i is the momentum forcing components and S ij is the viscous term given by where ρ and μ are the density and viscosity respectively, appropriating for the phase which is occupying the particular spatial location at a given instant. No turbulence model is included and surface tension is neglected. The numerical model is built under a fixed Cartesian grid system. The interaction of the dam-break wave and the downstream obstacle is treated as a multi-phase problem that includes structure, water and air. A volume function ϕ m is defined to constitute the dam-break wave interface, where m=1, 2 and 3 indicates water, air and solid, respectively. The total volume function for all computation domains is solved by the following advection equation, where ϕ 13 =ϕ 1 +ϕ 3 and u i is the velocity vector. The density and viscosity of both solid and liquid phase are assumed to be the same to ensure stability. The volume function for air ϕ 2 is determined by ϕ 2 =1.0-ϕ 1 -ϕ 3 . Thus, after all volume functions have been calculated, the physical property λ, such as the density and viscosity in each grid can be calculated by the following formula: The water depths in the following were the superposition of ϕ 3 .

Flow solver
Following Zhao and Hu (2012), the governing equations are discretized using a high-order finite difference method on a Cartesian grid system. The total discretizing processes are divided into three steps, which are advection phase, non-advection phase (i) and non-advection phase (ii). A fractional step scheme is used to solve the Navier-Stokes equations. The intermediate velocity is computed by the CIP method (Takewaki et al., 1985) firstly; then the pressure is obtained by solving the Poisson equation derived by enforcing the continuity constraint and the final velocity is updated by simple algebraic operations. The key elements of the numerical code are given in the following sections for the sake of clarity, and more details can be found in Zhao et al. (2014). ∂u 2.3 THINC/SW scheme for capturing free surface In this study, the THINC/SW scheme (Xiao et al., 2011) is applied for the free-surface treatment. The THINC scheme is a VOF-type method. The difference is that the original VOF method uses a Heaviside step function as the characteristic function, while the THINC scheme uses a smoothed Heaviside function, the hyperbolic tangent function. In the THINC/SW method, an interface orientation dependent weighting between the original THINC scheme and the first order donor cell scheme is suggested by Xiao et al. (2005), which significantly improves the geometrical accuracy of the original THINC scheme.
The improved multi-dimensional THINC/SW scheme is verified by Zalesak's rigid body rotation benchmark problem (Zalesak, 1979) which was widely used as a benchmark test for the VOF methods. Cases with grid of different sizes are carried out for parametric study. A velocity field is given by u=(y-0.5, 0.5-x) with Δt=2π/628. In general, one revolution is completed in 628 time steps. The numerical error is defined and estimated by: where, is the exact solution of ϕ i, j . The numerical errors are summarized in Table 1. It can be seen from the table that a finer grid produces better shape retention, and the numerical error of the THINC/SW scheme is smaller than that of the original VOF scheme and THINC scheme.

Treatment of solid body boundary
The Fractional Area Volume Obstacle Representation (FAVOR) method is applied to the solid body boundary condition. The FAVOR is one of the most efficient methods to treat the immersed solid bodies, which was developed initially by Hirt (1993). In the present numerical solution procedure, a forcing term is added to the momentum equation to impose the velocity distribution inside and on the solid body boundary, and an alternative way is to do the following updating after the calculation with Eq. (7).
where, is the local velocity of the solid body (here, is set to be zero for a stationary hump) and de- notes the flow velocity obtained from the fluid flow solver. Such treatment of boundary condition on embedded boundaries can be considered from the thought of the immersed boundary method. The advantage of this procedure is its simple extension to three dimensions.

Model validation
In this section, a dam break problem with a hump is computed and compared with physical results for the model validation. A set of laboratory experiments have been conducted by Ozmen-Cagatay et al. (2014), which has the dimensions of 8.90 m×0.30 m×0.34 m with a horizontal bottom, as shown in Fig. 1. The channel bottom and walls were made of glass with the thickness of 9 mm. A gate was located at 4.65 m from the channel entrance and separated the upstream part of the channel, representing the reservoir, from the downstream part. The initial water depth of the reservoir, h 0 , was 0.25 m and the downstream was dry. The removal time of the gate determined from the video records was between 0.06 s and 0.08 s which was shorter than 1.25(h 0 /g) 1/2 , satisfying a sudden removal definition (Lauber and Hager, 1998b). A triangular-shaped bottom obstacle with 1.0 m base length and 0.075 m height was placed 6.15 m downstream from the channel entrance (Fig. 1). Dye was added into the reservoir to color the water for better determination of the free surface levels from the recorded images.
The free surface profiles as well as the time evolution of water levels were obtained by a digital imaging system, which contains three high-speed cameras. A procedure was applied by using filter and edge recognition functions to ob-tain water level variations with time (h versus t) from images of the dam-break flow directly. Fig. 2 shows the initial conditions of numerical simulation. The computational time step is 10 -4 s.
In the numerical simulation, the computational domain is subdivided into a mesh of fixed rectangular cells using Cartesian coordinates. A non-uniform grid system is employed in the x direction and z direction, with the minimum uniform grid size of Δx=Δy=0.2/0.1/0.05 cm in the vicinity of the hump. The local mesh refinement region (uniform grid region) begins at x=4.05 m, z=0 m and finishes at x=7.15 m, z=0.20 m, which covers all locations P1-P7 in Section 3.4 and the hump. Difference only occurs in the uniform grid region and non-uniform grid regions do not change. Fig. 3 shows the experimental and numerical comparison of the free surface profiles at t * =15.16 and 17.54 (t * =t(g/h 0 ) 1/2 ) for three different grids. It can be observed that the computation results of Mesh 2 and Mesh 3 agree well with experimental results. Since Mesh 2 costs less computation time, it will be adopted in the following computations.
3.2 Spilling type reflected bore and the following spilling type wave breaking Fig. 4 shows the comparison of the computed and physical evolution of the dam-break configuration. When the gate is removed, a strong dam-break flood wave is generated which propagates over the dry channel (Fig. 4a). When the flood wave arrives at the hump, a part of the wave is reflected and a spilling type reflected bore appears travelling in the upstream direction while the other part travels towards the downstream of the hump. At t=0.09 s, the dambreak flood wave arrives at the top of the hump (Fig. 4b). At t=1.36 s, the dam-break flood wave just flows past the right toe of the hump (Fig. 4c). In this period the predicted water front agrees well with the experiment images qualitatively.
At t=2.8 s, the spilling type reflected bore begins to emerge (Fig. 4d). After t=3.3 s, the spilling type reflected bore begins to move in the upstream direction (Figs. 4e-4h). The spilling type wave breaking is observed in both experimental snapshots and the present computed results after t=3.3 s, while only the plunging type reflected bore and plunging type wave breaking were observed in the computed results by Ozmen-Cagatay et al. (2014). The water-air interaction can also be observed in both experimental snap-  It can be observed that the proposed method seems to present better accordance with experimental snapshots in the spilling type reflected bore, spilling type wave breaking and water-air interaction, compared with the methods mentioned by Ozmen-Cagatay et al. (2014). This might be due to the different types of VOF methods used in different numerical models. However, slight discrepancy exists in the propagation speed between the experiment observations and numerical predictions; the experiment bore has time delay. With the assumption of frictionless solid boundaries, the propagation of the spilling type reflected bore is over-predicted under the frictionless condition in the present model. The time delay caused by the gate motion in the experiment also causes such discrepancy. Fig. 5 shows the predicted and experimental free sur-   CHENG Du et al. China Ocean Eng., 2017, Vol. 31, No. 6, P. 683-692 face profiles of the dam-break flow at different time. The SWE (Shallow Water Equations) and RANS (Reynolds-Averaged Navier-Stokes equations) results (Ozmen-Cagatay et al., 2014) are also presented for comparison. For convenience, the results are reported in non-dimensional forms of the horizontal distance (X=x/h 0 ), water depth (Y=h/h 0 ) and time t * =t(g/h 0 ) 1/2 , where g is the acceleration of gravity. Free surface profiles are accurately predicted by all numerical simulations for t * =15.16 and t * =17.54 (Figs. 5a-5b). From t * =20.67 to t * =35.83, the negative bore (reflected wave) starts to move in the upstream direction and the spilling type wave breaking is observed. With early wave breaking at t * =23.05 (Fig. 5d), excellent agreements with experimental results are observed in the results of present model, while the SWE and RANS results are significantly different from the experimental data. It can be attributed to the ability of present model to simulate the complex flow phenomenon such as the spilling type wave breaking and airwater interaction that have been observed in the experiment. Different turbulence models k-ω, RNG etc, especially LES (Large Eddy Simulation) or a more accurate VOF method may be useful for the RANS model to capture the spilling type wave breaking. Compared with the results of SWE and RANS, the present model results show a minor free-surface fluctuation (1<x/h 0 <5) at t * =26.69 and t * =35.83 (Figs. 5e-5f). The possible reason is that the hydraulic jump on the free surface and the spilling type wave breaking in the experiment. Three numerical methods can accurately obtain the free surface profiles from t * =41.84 to t * =62.77 . During this period, the flow becomes more stable, the magnitude of the hydraulic jump and wave-breaking in front of the negative bore decreases due to the reduction of the incoming flow. Despite the discrepancies in the prediction of the negative bore front, the simple SWE solution can reproduce the wave profiles. The non-conservative form of SWE in FLOW-3D, neglecting the vertical acceleration and assuming hydrostatic pressure may be attributed to the negative bore difference in the SWE solution.

Comparison of free surface profiles
3.4 Water elevation along the channel Fig. 6 shows the six measurement locations of stage hydrographs along the channel (P1-P6). In the experiment (Ozmen-Cagatay et al., 2014), water elevation variations over time were determined directly from the digital video images of the dam-break flood wave propagation by using virtual wave probes at selected sections. Fig. 7 depicts the comparison between the computed and measured time evolutions of water depths at six selected measuring gauges from P1 to P6. The SWE and RANS results (Ozmen-Cagatay et al., 2014) are also presented for comparison. Agreement between the numerical and measured results is generally satisfactory.
Gauge P1 (Fig. 7a) and Gauge P2 (Fig. 7b) are located just upstream and downstream of the dam gate. When the gate is removed, a strong flow is generated and propagates downstream, while abrupt drop and rise in water depth are observed at P1 and P2 in both numerical results and experimental data (t * <5). There is little change in water depth at  P1 and P2 when t * =5-40. An abrupt rise in water depth at t * =40 as the negative bore from the obstacle can be observed in both the numerical results and experimental data. The water depths gradually drop due to the reduction of the incoming flow. From Figs. 7a and 7b, the present model slightly overestimates the water levels at t * =40 in Fig. 7a and t * =35 in Fig. 7b where the negative bore just reaches the water level measurement positions. The reason for this discrepancy is that water levels were determined directly from the digital video images, ignoring the splash of the dambreaking in front of the negative bore in the experiment.
At Gauge P3 (Fig. 7c), a rising of the water depth from t * =2 to t * =29 can be observed, and a sudden rise due to the arrival of the negative bore is seen at t * =29. The results at t * =30 of the present model for P3 reveal similar phenomena to those at P1 and P2. Some water level fluctuations are seen when t * =32-40. This may be due to the appearance of the spilling type wave breaking in front of the negative bore and the splash on the free surface.
At Gauge P4 (Fig. 7d), the present model, SWE and RANS show a sharp and rapid rising in water depth during the passage of the negative bore through the measurement position. A possible reason for this is that the splashes of the wave breaking in front of the negative bore and the hydraulic jump on the free surface are ignored in the experimental measurements. Since the rapid rising of the water level occurs behind the negative bore at Gauge P5, there is no elapsed time for travelling between flood positive flow and negative flow. Therefore, water levels linearly and constantly rise up to their maximum levels. At the measurement locations close to the toe of the hump (P4 and P5), the maximum water level occurs. At Gauge P6 (Fig. 7f) both the computational and experimental results show a general agreement. All the computed water depths on the falling limbs of all graphs (P1-P6) are always in good agreement with the experimental data.

Numerical results for dam breaking wave with varying slopes
The slope of a natural hump is always gentle, while that of the artificial hump is always steep. Based on the capability of this method mentioned above, numerical calculations were performed for fourteen cases with different slopes (up1-7 and down1-7) to investigate the impact of varying slopes on time evolution of water level at P1-P6 (six selected measuring gauges in Fig. 6).
In our numerical simulation, we first studied the effect of the upside slope (varying from gentle to steep) on non-dimensional time (T=t(g/h 0 ) 1/2 ) evolution of water level (Y=h/h 0 ) at P1-P6. The results are shown in Fig. 9. Some interesting results can be found in Fig. 9. First, the time delay of the negative bore is more serious at P1-P5 with the increase of slope when the slope of the hump is gentler than that of up5. Second, the negative bore arrives at P1-P5 much earlier with the increase of slope when the slope of the hump is larger than that of up5. Third, a good agreement is achieved on falling limbs of all graphs (P1-P6) for all slopes. When the reflected wave moving towards the upstream direction reaches the end of the dam, inflow begins to fall due to the finite length of the reservoir.
According to the different forms of the wave front caused by the hump, the slope can be divided into three types: gentle slope (up1-4), turning slope (up5) and steep slope (up6-7). The lower surface of the wave front touches with the hump completely under the gentle slope (Fig. 10a). Under such condition, as the wave front rises at the upside of the hump, a part of the kinetic energy of water front is turned into potential energy and thermal energy. This means that the gentle slope definitely leads to a velocity reduction of the wave front. At the top of the hump, the minimum velocity occurs. As the water with low velocity hinders the following dam breaking wave, spilling type negative bore appears. In this case, the friction and hinder effect of the hump are the main causes of the negative wave. Therefore, the time delay of the negative bore is more severe at P1-P5 with the increase of slope when the slope of the hump is gentler than that of up5.
Under the condition of the turning slope (near up5), big bubbles are closed by the lower surface of the wave front at the downstream of the top of the hump (Fig. 10b). With the emergence of those big bubbles, a more serious hinder effect occurs. Those bubbles become a key reason of hinder effect and will inevitably lead to a more delayed and violent negative bore. Then, the closer the slope is to up5, the more serious the time delay will become and a larger scale of negative bore will emerge.
With the influence of steep upside slope, hydraulic jump is found at the downstream of the hump (Fig. 10c). In such cases, as the bubbles and the low velocity water disappear, the reflecting action of steep upside slope plays a key role in the generation of the negative bore. The steeper the slope, the more obvious the reflecting action becomes. Thus, the negative bore arrives at P1-P5 much earlier with the increase of slope when the slope of the hump is larger than that of up5.
As the upside slope is only one of the key factors affecting the spilling type negative bore (speed of advance and scale), similar consistency is achieved on the falling limbs of all graphs (P1-P6) for all slopes (Fig. 9) when the reflected wave moving towards the upstream direction reaches the end of the dam, and the inflow begins to fall due to the finite length of the reservoir. Fig. 11 shows the computed results for non-dimensional time evolution of water level at P1-P6 (down1-7, the downside slope of the hump changes from gentle to steep with the upside slope is constant). The results show similar consistency with different slopes, and only little difference was found when the negative bore arrives at the measured points (P1-P6). The possible reason for such difference may be the spilling type reflected bore and the following spilling type wave breaking. The results indicate that the downside slope is less important for the formation of the negative bore.

Conclusions
A CIP-based multi-phase fluid model has been applied to simulate the dam-break induced tsunami-like bore with a hump of varying slopes. The model governed by the Navier-Stokes equations is solved by a CIP-based high-order finite difference method based on a fixed Cartesian grid system. A more accurate THINC/SW scheme is adopted for the free surface capturing. Moreover, an immersed boundary method (IBM) is used for the treatment of solid body boundary.
The locations of water front, the complex spilling type negative bore, spilling type wave breaking and experimental profiles are properly reproduced by the model through the comparison of evolution of dam-break configuration between the present model results and experimental data. Results indicate that the CIP-based model has shown a more powerful ability to deal with the violent spilling type reflected bore, spilling type wave breaking and air-water interaction compared with SWE and RANS results.
Based on the validated model, computations were per-  formed for fourteen cases (up1-7 and down1-7) to investigate the effect of varying slope of the hump on the features of a dam-break flow. According to the different forms of the wave front caused by the hump, the slope of the upside hump can be divided into three types: gentle slope (up1-4), turning slope (up5) and steep slope (up6-7). The major factor that leads to the appearance of negative bore is different from the change of the upside slope. The closer the slope is to the turning slope (up5), as the influence of the bubbles, the more serious the time delay will become and a larger scale of negative bore appears. Similar consistency is achieved on the falling limbs of all graphs (P1-P6) for both slopes when the reflected wave moving towards the upstream direction reaches the end of the dam, and the inflow begins to fall due to the finite length of the reservoir. Numerical results indicate that the downside slope is not a key factor that leads to the appearance of negative bore.
Results of the research prove that the present model is powerful to deal with the complex behavior in the interaction between tsunami-like bore and structure. Meanwhile the conclusion provides reliable references for the design and improvement of the offshore seadike which plays a crucial role in protecting the people and their property from tsunami and surge.