### 1. Introduction

*b*,

*s*). Here

*b*and

*s*are vertical locations from a reference line. In river flow or debris flow, people often use the Saint-Venant equations with the variable of depth of water or debris instead of

*b*and

*s*. When using the shallow water equations with the coordinate of (

*b*,

*s*), it is possible to simulate flow in river and coastal zones simultaneously. We use the finite volume method to discretize the governing equations to deal with numerical discontinuities. In section 2, we develop a model of subaerial and submarine landslides based on the shallow water equations with the coordinate of (

*b*,

*s*) and the finite volume method is used to discretize the governing equations. In section 3, we validate the developed model by comparing numerical results with analytical solutions and experimental data for dambreak water flow and debris flow. In section 4, we apply the developed model to subaerial and submarine landslides and find reasonable results. In section 5, we summarize the present investigation and suggest future study.

### 2. Development of Landslide Model

*b*,

*s*) and the finite volume method is used to discretize the governing equations.

### 2.1 Derivation of governing equations

*s*and

*b*are the elevations of the surface and bottom of the debris from a reference line,

*h*(=

*s*−

*b*) is the debris depth, and (

*u, v, w*) are the velocities in the (

*x, y, z*) directions. The continuity equation for the incompressible fluid flow is given by

*x, y*) domain as

*x*) domain as

_{xx}, τ

_{yx}, τ

_{zx}are the shear stresses. We only consider the

*x*−

*z*plain, so τ

_{xx}, τ

_{yx}are equal zero. In Eq. (6), the pressure

*p*is assumed to be hydrostatic and thus

*p*=

*ρg*(

*s*−

*z*). We add

*u*(

*∂u*/

*∂x*+

*∂v*/

*∂y*+

*∂w*/

*∂z*) to Eq. (6) in order to get the following equation

*∂b*/

*∂x*(=

*S*

_{0}) is the bottom slope. We integrate Eq. (7) from the bottom to the surface, use the Leibnitz’s rule and the boundary conditions given by Eqs. (2) and (3). After neglecting the shear stress at the free surface, the resulting momentum equation is given as

*x*) domain, Eq. (8) is reduced to

*S*may consist of two components, i.e., the bottom friction and the internal soil friction given by

_{f}*n*and

*μ*are the Manning’s roughness and the effective friction coefficients, respectively. Manning’s

*n*values range from 0.032 to 0.2 s/m

^{1/3}while the

*μ*values for debris flows range from 0.06 to 0.175 (Paik, 2015).

*s*−

*b*), the resulting momentum equation is given by

### 2.2 Numerical scheme

*U*is the vector of conserved variables. (Eq. (14) is solved using the splitting method. At the first step, the homogenous pure advection problems are solved as

*U*are obtained by adding the source terms as

*S*is evaluated. To avoid numerical instability due to the friction term, an implicit scheme is employed following Paik (2015).

*c*(=

_{k}*S*Δ

_{k}*t*/Δ

*x*) is the Courant number related to the wave speed

*S*. The flux jump across wave

_{k}*k*at the cell interface is given by

*r*

^{(k)}is the ratio of the upwind change to the local change.

*B*(

*r*

^{(k)}) is the limiter function (Toro, 2001). In this study, we use the van Leer limiter function for

*B*(

*r*) as

### 2.3 CFL criterion

*c*is the celerity, uf is the total friction velocity. In the Voellmy resistance relation, the total friction velocity can be expressed as

### 3. Model Validation for Water Flow and Debris Flow

### 3.1 Validation for water flow after dam break

*t*= 0.2 s using the conservative forms and alternative conservative forms, respectively. In the figures, analytical solutions (Mangeney et al., 2000) as well as initial conditions are shown for comparison. The conservative forms of the equations yield the same solutions as the analytical solutions. However, the alternative conservative forms of the equations yield different solutions. At a point of dry bed, (

*s*−

*b*) is equal to zero. In deriving alternative conservative forms of the momentum Eq. (13), we divide some terms by zero-valued (

*s*−

*b*) which is mathematically meaningless. That is why numerical solutions of alternative form of the equations are different from the analytical solution in the dry/wet discontinuity points. In further investigation, we use the conservative forms of the equations only.

*t*= 0.12 s, 0.20 s, 0.36 s, 0.50 s. On the whole, the numerical results show good agreements with the analytical solutions and the experimental data even though wave breaking occurred in reality. Especially, the numerical solutions are closer to the measured data at

*t*= 0.36 s and

*t*= 0.50 s when the water surface elevations are milder than at previous times.

^{−1/3}s, 0.0125m

^{−1/3}s, 0.01875 m

^{−1/3}s which are different with a ratio of 1.5 each other. The upstream boundary is closed while the downstream boundary is open. Liang and Marche (2009) measured time series of water depths at 7 gauges which are located at 2 m, 4 m, 8 m, 10 m, 11 m, 13 m, and 20 m downstream of the dam. For numerical experiment, the total simulation time is 90 s and the grid size is 5 cm.

*t*= 5 s. Because of an interaction between the incoming waves and the reflecting waves from the side of the hump, a shock wave formed and propagated back to the upstream boundary. On the other side of the hump, a rarefaction wave occurred and moved to the downstream boundary, which caused increase of water depth at the right side of the hump (gauge 7). Around

*t*= 25 s, the shock wave reached and reflected from the wall at the upstream boundary. After this reflection, the second wave propagated to the downstream. These complex processes continued until the momentum energy was damped due to bottom friction and water mass passed through the downstream boundary.

*t*= 90 s. At gauges 3 to 6, numerical results around the first peak of water depths are higher than the measured data. At gauge 7 which is at the other side of the hump, numerical results are always higher than the measured data. These discrepancies might happen because the numerical solutions of the shallow water equations were made based on the assumption of neglecting dispersive terms which would spread out several components with different speeds. If we simulate these flow phenomena using the Boussinesq equations which include dispersive terms, we would get solutions closer to the experimental data. When the Manning’s roughness coefficient n increases, the flow speed decreases at the first wave and the flow thickness increases at all the 7 gauges. Especially, at gauges 3 to 6, due to the adverse slope of the hump, with the increase of the Manning’s n, the flow thickness decreases at the first peak. The Manning’s n equal to 0.0125 gives the numerical solutions closest to the measured data.

### 3.2 Validation for debris flow

^{3}consisting of sand-gravelmud (SGM) material and is suddenly opened by a vertical headgate. More detailed information about the experimental facilities was reported in Iverson (2003) and Iverson et al. (2010).

*ρ*= 2100 kg/m

^{3}, the Manning’s roughness coefficient of n = 0.034 m

^{−1/3}s, and the effective friction coefficient of

*μ*= 0.12. The grid size is 0.1 m and the CFL number is 0.9 which guarantees numerical stability. During simulation we use a constant value of 2 m channel width because the numerical model is horizontally one dimensional. In the real experiment, however, the width changed from 2 m at the main channel to 4 m at the runout surface. This width change would cause numerical solutions different from the experimental data on the area of the runout surface.

*x*= 32 m, 66m, 90 m downstream of the headgate. It should be noticed that the location x is not the horizontal distance but the downslope distance on the bottom. The first and second locations are in the main channel while the third location is on the runout surface. Due to turbulent high speed of the debris flow on the steep channel, high fluctuation of the solutions occurred at the first and second locations in the present high-resolution numerical model. To smooth the highly fluctuating solutions, we use the moving average algorithm as given by

*y*is the raw (noisy) data, (

_{k}*y*)

_{k}_{s}is the smoothed data, the odd number 2n + 1 is named as the filter width or span. The greater the filter width is, the smoother the data are. In order to smooth these time series data of the flow thickness, we choose the value of span 2

*n*+ 1 = 101 which is equivalent to value of n = 50 and repeat smoothing process again. The results after smoothing are shown in Figs. 6(a) and 6(b).

*x*= 32 m and

*x*= 66 m in Fig. 6(a) and Fig. 6(b), respectively, show better agreements than those at

*x*= 90 m. These differences are of the curvature connected between bed slope of 31° and the flat bottom of 2.4° is not represented well in this simulation. Moreover, in the present 1D model, the width of the channel of 2 m is taken into account for the whole domain while that of runout surface of 4 m is set at this point of the experiment as mentioned previously. It results in the overestimation of approximate 50% in term of the flow thickness on the runout surface.

### 4. Model Application to Subaerial and Submarine Landslides

*g*

_{0}= 9.81 m/s

^{2}whereas the submarine landslide is affected by the buoyancy gravity acceleration

*g*

_{1}=

*g*

_{0}(

*ρ*−

_{s}*ρ*/

_{w}*ρ*where

_{s}*ρ*and

_{s}*ρ*are densities of soil and water, respectively. Rzadkiewicz et al. (1997) conducted hydraulic experiment of a submarine landslide and the induced tsunami. The volume of soil in a triangular section with 0.65 m × 0.65 m suddenly slid down on a slope of 45° which was connected to a horizontal bottom and the tsunami occurs due to the landslide. They also simulated the phenomenon using 2-dimensional (

_{w}*x*,

*z*) Navier-Stokes equations. In this study, we conduct numerical experiments of subaerial and submarine landslides using their model scales (see Fig. 7). However, we use different properties of

*ρ*= 2100 kg/m

^{3}, n = 0.034 m

^{−1/3}s,

*μ*= 0.08 following Paik (2015) who simulated landslides using the Saint-Venant equations. The grid size is Δ

*x*= 1 cm and the CFL number is 0.9 which guarantees a numerical stability.

*t*= 0.4 s. For the submarine landslide, due to the buoyancy effect, the acceleration is reduced to

*g*

_{1}= 0.52

*g*

_{0}and the soils move more slowly compared to the subaerial landslide.

*t*= 0.4 s for the water surface at

*s*= 0.8 m. Subaerial landslide occurs above the water surface while submarine landslide occurs below the water surface. For the submarine landslide, due to the buoyancy effect, the acceleration is reduced to

_{w}*g*

_{1}= 0.52

*g*

_{0}and the soils move more slowly compared to the subaerial landslide. Therefore, for this combined subaerial-submarine landslide, the soils move more slowly compared to the subaerial landslide (Fig. 8). However, they move faster compared to the submarine landslide (Fig. 9).

*θ*= 45° to

*θ*= 30°. Fig. 11 shows the topography for simulating subaerial landslide with bottom slope of θ = 30°.

*t*= 0.4 s. After comparison with the solution with θ = 45° shown in Fig. 8, we find that the velocity increases with the increase of the bottom slope. These results are physically reasonable.

^{−1/3}s and n = 0.051 m

^{−1/3}s which are different with a ratio of 1.5. It is clear that the propagation speed of the debris flow depends on the bottom roughness coefficient.

^{−1/3}s and 0.051 m

^{−1/3}s, respectively. The speeds of landslides of the two cases are slightly different at

*t*= 0.2 s, but the two speeds become more clearly different at

*t*= 1.0 s. This implies that larger bottom roughness cause slower speed of the debris flow which was mentioned also by the previous works (Mikoš et al., 2006; Paik, 2015).

### 5. Conclusions

*x*) shallow water equations. We used the (

*b*,

*s*) coordinate system to express the shallow water equations which can be applied in coastal zone as well as river. To discretize the governing equations, we used the well-known high-resolution scheme which combines the TVD limiter function and the second-order Godunov method (WAF method) employing HLL approximate Riemann solver. This scheme can deal with the wet/dry problems occurring in the real landslide phenomena. We applied both the conservative and alternative conservative forms of the shallow water equations to a dam-break water flow and found that the conservative form of the equations yielded accurate solutions in the dry/wet boundaries but the alternative conservative form yielded inaccurate solutions. We further verified the developed numerical model by comparing numerical results with the analytical solutions and the experimental measurements for both the dam-break water flow and the debris flow. We also applied the model to subaerial landslides with different bottom slopes and bottom frictions and found that both steeper bottom slope and smaller bottom friction cause faster movement of soils. In the future, the two-dimensional (

*x*,

*y*) shallow water equations as well as the non-hydrostatic equations (e.g., the Boussinesq equations) will be developed to accurately solve the subaerial and submarine landslide phenomena.