MPM–FEM hybrid method for granular mass–water interaction problems

The present study proposes an MPM (material point method)–FEM (finite element method) hybrid analysis method for simulating granular mass–water interaction problems, in which the granular mass causes dynamic motion of the surrounding water. While the MPM is applied to the solid (soil) phase whose motion is suitably represented by Lagrangian description, the FEM is applied to the fluid (water) phase that is adapted for Eulerian description. Also, the phase-field approach is employed to capture the free surface. After the accuracy of the proposed method is tested by comparing the results to some analytical solutions of the consolidation theory, several numerical examples are presented to demonstrate its capability in simulating fluid motions induced by granular mass movements.


Introduction
In the natural phenomenon known as submarine or/and subaerial landslide [25,26], sediment on the slope slides down from shallower to deeper regions of the ocean, resulting in a series of hazardous events, such as the severing of submarine cables and submarine pipelines, obstacles obstructing the operation of resource development facilities, and large-scale tsunamis. The damage caused by such a tsunami has been considerable on a global scale, and the loss of life has also been significant. However, the process of such a landslide is not well understood. Since multiple physical phenomena simultaneously interact in the generation of tsunamis caused by a landslide (the collapse of the seabed, the interaction between soil and water), determining the mechanism has proven to be challenging. It is, however, extremely difficult to predict the behavior of the soil and seawater with high accuracy using conventional approaches [13]. Therefore, it is necessary to develop effective numerical simulation methods that can properly evaluate the complex mechanical behaviors involved in the generation of tsunami due to landslides.
Among the computational methods to solve the complex solid-liquid coupled problems, the traditional finite element method (FEM) is not good at solving the complex large deformation in solid phase such as debris flow and landslides due to severe grid distortion. In this context, the grid-free methods (including FEM with the Eulerian frame) have been developed as substitutes for FEM. The element-free Galerkin method [29], smooth particle hydrodynamic (SPH) method [8,16,37], multi-phase flow model [19,30,40], Particle-FEM [18], and material point method (MPM) [1,3,38] have recently become popular in this area of research. Notable among these methods is the MPM [32] for its contribution to remarkable progress in this area, and is the method of choice for many researchers.
Within the MPM framework, two approaches are proposed to study the interaction between solid and liquid phases: (1) the two-phase single-point approach [10,11,15,31,39] and (2) the two-phase double-point approach [1,3,22,23,38]. The governing equations used by the two approaches are not the same. Compared with the first approach, in which the formulation is derived by reference to the solid configuration, the second approach uses two sets of material points to separately represent the solid skeleton and the liquid phase. Thanks to this, the motions of each phase can be expressed separately and applied to solve more diverse problems, such as seepage [3,23], over-topping [22] and transportation [38]. Since we have our sights set on simulating a tsunami generated by a submarine or/and subaerial landslide, which requires the rough movements of water and soil phases to be represented, the second approach is considered more suitable.
However, in most of the coupled solid-liquid MPM algorithms, see, e.g., [3,22], the liquid phase is assumed to be weakly compressible, since the MPM was originally developed mainly with focus on large deformation problems for solids. Many of these existing analysis methods result in pressure oscillation, free surface instability and involve a high computational cost due to the large bulk modulus of water when an explicit time integration is adopted. For single phase problems, there have been some remedies for dealing with incompressibility that yield stable solution schemes for MPM calculations; see, e.g., [17,42] for flow simulations of incompressible fluids. Nevertheless, only few attempts have so far been made at managing with incompressibility for solidliquid MPM. In addition, the dynamic addition or deletion of water particles is more complicated when it comes to problems containing inflow or outflow boundaries [22,43]. In view of such issues, a class of stabilization schemes with a Eulerian description provides a simpler process and therefore has clear advantages in terms of calculation cost and accuracy. Thus, a novel scheme needs to be developed to couple the FEM with MPM to simulate granular mass-water interactions. Here, the "granular mass" is a synonym for "soil structure" in this study.
In this study, we propose a hybrid method of MPM and FEM capable of expressing the complex interaction between granular masses and liquids. The proposed method is formulated based on porous media theory [4,6], in which saturated soil is considered as a continuum composed of the soil skeleton and the liquid phase. The MPM is applied to the governing equation of the solid phase by the Lagrangian description, whereas the stabilized FEM [34] is applied to the governing equation of the liquid phase in the Eulerian frame. In addition, a free surface is captured by using the phase-field method [14]. Several numerical examples are presented to demonstrate the performance of the proposed hybrid MPM-FEM on suppressing the afore-mentioned instabilities in the liquid phase (pressure oscillation, free surface instability) often appearing in the conventional solid-liquid coupled MPM and the capability in dealing with landslide-triggered tsunamis involving complex granular mass-water interaction problems.

Governing equations
In this study, saturated soil is assumed to be a continuum consisting of a solid skeleton and a fluid phase, which are separately handled as shown in Fig. 1. It must be noted that either fully saturated or fully dry soils are considered in the present formulation.

Momentum balance and continuity equations
In the present formulation, the superscripts 's' and 'f' are attached to involved variables and symbols to refer to the solid skeleton and fluid phase, respectively. In line with the methods established in a previous study [22,23,38], the momentum balance equations are solved to obtain the accelerations of the solid skeleton and the fluid phase separately. The momentum exchange term is considered in both of them.
In mixture theory, the momentum balance equation for the solid skeleton is written as where ρ s is the true density of solid grains, n is the porosity, and b α i (α = s, f) is the body force. Here, the partial stress tensor of soilσ s i j is defined as where σ i j is the effective stress of the soil skeleton, p f is the fluid pressure and δ i j is the Kronecker's delta.P i is the drag force between the solid skeleton and the fluid phase, and the Forchheimer-like equation [23] is used to compute it aŝ where μ f is the dynamic viscosity of fluid and k is the intrinsic permeability. The intrinsic permeability of the soil is calculated according to the following Kozeny-Carman formula [9]: which is a function of the porosity and the averaged diameter of the solid grains denoted by D 50 . Using Eqs. (1)-(3), the final version of the momentum balance equation for the solid skeleton is written as The momentum balance equation for the fluid phase is expressed as where ρ f is the true density of fluid, andσ f i j , a f i are the partial stress tensor and acceleration of fluid can be written as, respectively, Here, the a f i is in the Eulerian description so the advection velocityv f j also needs to be considered. Using Eqs. (6)-(8), the final version of the momentum balance equation for the fluid phase becomes On the other hand, solid grains can be regarded incompressible, the mass balance equation for the solid skeleton is given as (10) in which the second term can be neglected. Also, on the assumption of the incompressibility the fluid phase, the mass balance equation for the fluid phase is given as in which the second term can also be neglected. In this study, the continuity equation of the fluid is solved together with the momentum equation of the fluid. Therefore, the material time derivatives on the left-hand side of Eqs. (10) and (11) read Using Eqs. (10)- (13), the continuity equation for the whole mixture becomes

Constitutive model of soil skeleton
We employ a standard elasto-plastic constitutive model based on multiplicative decomposition of the deformation gradient such that where F e ik and F p k j are the elastic and plastic components, respectively. The Hencky hyperelastic model is adopted to represent the elastic response of the solid keleton whose energy function is given as where λ e α (α = 1, 2, 3) are the three principal stretches of the elastic component of the left-Cauchy-Green deformation tensor, λ, μ are Lame constant, and J e is the determinant of F e ik . The plastic deformation of the solid phase is defined by the Drucker-Prager yield criterion whose plastic potential function is given as where s i j is the deviatoric stress defined as s i j = σ i j − 1 3 pδ i j . Here, p is the hydrostatic pressure component and J 2 is the second invariant of the deviatoric stress tensor. Also, c is the cohesion, and η, ξ are material constants defined as η ≡ 3 tan ψ 9 + 12 tan 2 ψ , (18) ξ ≡ 3 where ψ is the internal friction angle. In order to suppress excessive dilatancy, the non-associated flow rule is employed in this study. The corresponding plastic potential to determine the specific flow rule is formed by setting the second and third terms in the right hand side in Eq. (17) to zeros.

Weak forms
We adopt the direct method with the conventional FEM, in which the momentum balance and continuity equations, Eqs. (9) and (14), are solved together. The corresponding weak form is written as where v f i ,v f j , p f andp f i are the velocity, advection velocity, pressure fields and prescribed pore pressure that acts on the part of the fluid phase boundary ∂Ω f , respectively, while ω i and q are the corresponding test functions.
By multiplying Eq. (5) by a test function ω i and then integrating over the domains Ω s and Ω f , we also can derive the weak form of the momentum balance equation for the solid phase as in which the term involving the stress has been integrated by parts by the application of divergence theorem. Here, t i is the prescribed traction that acts on the part of the solid skeleton's boundary ∂Ω s , andp f i is the same as that in Eq. (20).

Hybrid MPM-FEM
In the proposed hybrid MPM-FEM, MPM is applied for deformation analyses of the solid skeleton using both Lagrangian and Eulerian frames, the latter of which is used together to apply FEM for flow simulations of the fluid phase. This section presents both of the discretized equations in some detail. For reference, Fig. 2 shows a schematic diagram illustrating the combination of these separate discretization schemes.

Spatial discretization
An arbitrary variable (x) in the solid skeleton and the fluid phase is approximated as integrated over every FE (finite element) and every solid material point as, respectively, where superscripts 'e' and ' p' represent a FE and material point, respectively. Also, Ω e is the domain of an element and Ω p is the volume of a solid material point. The field variables such as velocity, pressure and their test functions can be approximated using the shape functions N α (x) associated with the Eulerian grid as where N n is the total number of nodes and α is the nodal value at node α. Here and hereafter, variables with superscript 'h' indicates their approximations in this manner. It should be noted that the same shape function and the same Eulerian grid are utilized to discretize both the variables associated with the solid skeleton and fluid phase. In this study, we employ the 2nd-order B-spline basis functions as shape functions, the details of which are presented in Appendix A. The stabilized FEM with the SUPG/PSPG stabilization schemes [7,34] is applied to the governing equation of the fluid phase, Eq. (20), to obtain the following discretized equation: The first four integral terms on the left-hand side and the integral term on the right-hand side are Galerkin items.
On the other hand, the fifth term involves the SUPG and PSPG terms, which are introduced to stabilize the advectioninduced unstable behavior and to suppress pressure oscillation, respectively. Also, the sixth term is the shock-capturing term [2] introduced to avoid the numerical instability of free surfaces. The stabilization parameters, τ e m and τ e cont , are defined as, respectively, where Δt, h e , ν,v f, h i and Re e are the time increment, the characteristic element length, the kinematic viscosity, the advection velocities, and the element Reynolds number, respectively.
Substituting Eqs. (22) and (24) into Eq. (25), we obtain the following FE equations for the fluid phase: which constitutes the set of the momentum balance and the continuity (mass conservation) equations, respectively. These discretized equations are to be solved for the nodal fluid velocity and pressure vectors at nodes, provided that the nodal solid velocity vectors are given as data. Here, M, A, G, Q, S c , F and C are the mass matrix, advection matrix, pressure matrix, interaction matrix, shock-capturing matrix, force vector, and continuity matrix, respectively. Also, subscript 'i' indicates the component of a nodal velocity vector in the x i -direction, and subscripts 's' and 'p' indicate the SUPG and PSPG matrices, respectively. The specific forms of these matrices are written as follows: Here, subscripts ' j' and 'l' indicate the components of the nodal velocity vectors, 'α', 'β', 'γ ', and 'δ' denote the values of the grid nodes. Also, the volume fraction of the solid material points in an element, ϑ, has been introduced, because there is a possibility that the fluid phase and the solid skeleton do not overlap with each other. In these discretized equations, n, k, and ϑ are functions of space to be interpolated by use of the shape functions as where n α and k α are the nodal values of the porosity and permeability, respectively. On the other hand, substituting Eqs. (22)∼(24) into Eq. (21), we obtain the semi-discrete momentum balance equation for the solid skeleton as This equation is solved for the solid acceleration at nodes by using the pressure of the fluid phase as input a datum, which is the solution of the momentum balance equation for the fluid phase in Eqs. (30) and (31). Here, the mass matrix for the solid skeleton is lumped as where m p is the mass of the solid particle computed as m p = (1 − n p )ρ s Ω p . Also, the nodal values of the internal and external forces are written as, respectively, in which the fluid pressure must be interpolated as

Time integration
The semi-discrete equations derived above are discretized with a discrete set of times, and the numerical solution is obtained at a current time. In what follows, the discrete approximation at time t L is indicated by superscript L and the time increment is defined as The time discretized forms of the momentum balance and the mass balance equations for the fluid phase, Eqs. (30) and (31), can be written as which are simultaneously solved for the fluid phase velocities v f, L+1 i and fluid pressure p f, L+1 . Here, the range of θ is set to be 0 ≤ θ ≤ 1. Also, the nodal variables (time derivative, velocities, pressure) are approximated as where θ has been set at 1 in Eqs. (58) and (61), 1/2 in Eq. (59) and 0 in Eq. (60). To explicitly approximate the advection velocitiesv f i , we adopt the Adams-Bashforth method [12], as Also, the continuity term where θ has been set at 1 in Eq.  56) and (57) as from which the fluid phase velocity and pressure are determined for the next time step. The time-discretized form of the momentum balance equation for the solid phase (52) can be written as Using the nodal pressure, p f, L+1 , obtained from Eqs. (67) and (68), this is solved for a s, L i explicitly with θ = 0. After the nodal acceleration of the solid skeleton, a s, L i , is obtained, the nodal velocity vector in the and then the velocity and position vectors of each solid material point can be updated by the following formulae: Prior to updating the state valuables of each solid material point, we adopt the MUSL procedure [32] to "refine" the nodal velocity vector by an additional mapping process. Then, the deformation gradient of a solid material point can be computed as and its determinant J p, L+1 = det F p, L+1 is used to update its volume and the porosity as, respectively, According to Eq. (4), the permeability k p, L+1 can also be calculated using n p, L+1 as Here, the upper limit value of the solid porosity n p, L+1 has to be set to suppress numerical instability. In this study, it will be set at n p, L+1 = 0.99 in the numerical examples.

Mapping of state variables at node
Nodal values for the solid skeleton are calculated by the mapping procedures established in the previous studies [3,38].
As is the case when using the conventional solid-fluid coupled MPM, the weighted least square approximation is adopted to evaluate the nodal velocity vector and the permeability at each node as v s, L+1 The volume fraction of a grid node, ϑ L+1 α , is computed using Here, Ω α is the nodal volume defined as which is carried out element-wise with the Gaussian quadrature rule [41]. Similarly, we map the porosities of solid particles onto the value at a node such that to suppress numerical instability due to the interpolation of the porosity gradient [38].

Interface representation
It is assumed that the solid phase is homogeneous and that the fluid phase consists of liquid (water) and gas (air) phases, which are identified by superscripts " " and "g" on their variables and symbols, respectively. Under this circumstance, the fluid phase has an interface with the solid skeleton or another fluid phase within an element, as illustrated in Fig. 2, and moves with time. The moving interface must be carefully handled to determine the domains of the phases involved so that the momentum balance and continuity equations derived above can properly be evaluated.
There are two kinds of methods to determine the geometry of a free surface that is an interface between either a liquid and a gas or a fluid and a solid. One is the interface tracking method, which uses a Lagrange description to directly express the geometry using a moving grid. The other one is the interface-capturing method, which uses an Eulerian or spatially fixed grid to track the time evolution of the interface. Since our target problems involve wave generation and propagation, and complex free surfaces, we employ the phase-field method, which is an interface-capturing method in this study. The details of this capturing method are provided in Appendix B.

Numerical algorithm
The numerical algorithm for the proposed hybrid MPM-FEM is summarized as follows:

Some remarks
It is worth mentioning that both the FEM and MPM adopt weak formulations and therefore have almost the same mathematical structures. This is particularly advantageous for solid-liquid coupled problems, because the domains of solid and liquid phases overlaps with each other in most problems of interest. In fact, since the numerical algorithms of MPM are inevitably similar to that of FEM, the information exchanging between solid and liquid phases is easier than other CFD approaches such the finite difference method.
Nevertheless, surprisingly few studies have so far been made at combining FEM and MPM for solid-liquid coupled problems, for which liquid and solid phases are separately described in the Eularian and Lagrangian frames, respectively. Although there have been some studies to connect FEM and MPM with different goals in the literature [20,21], they employ a sequential coupling method to transition between a moderately large deformation state that can be dealt with FEM and an extremely large deformation state that can more suitably be handled by MPM. Thus, the novelty of this study is obvious.
A remark is also made on the computational cost of time integration. Since the explicit time integration scheme is adopted for the solid phase, the time step size controlled by the CFL condition must be sufficiently small. Even though the implicit time integration is employed for fluid phase, the overall time step size must be controlled by that for the solid phase. Nevertheless, the size needs not be as small as that of the existing solid-liquid MPM that assumes large bulk modulus of elasticity to represent the weakly compressible liquid phase. This is another advantage over the previous explicit MPM-MPM framework for solid-liquid coupled problems.
However, the proposed method is still time-consuming and therefore not suitable for large-scale and long-duration phenomena that might include submarine landslide-induced tsunamis. This is one of the major limitations and remains to be resolved in future studies. The other major limitation is that the formulation is only applicable to either fully saturated or dry soils but not to unsaturated ones. Therefore, some important phenomena such as rainfall-induced landslides cannot be handled by the present scheme and give us a significant challenge.

Representative numerical examples
In this section, several numerical examples are presented to demonstrate the capability of the proposed hybrid MPM-FEM. First, the formulation is verified by solving the onedimensional consolidation problem in comparison with the Terzaghi's analytical solution. Second, a model experiment of the wave generated by an underwater granular collapse is simulated to assess the performance of the proposed method. Also, the model in which half of the soil column soaked in water and the other half in air is simulated to further explore the capability of our method. Finally, a model experiment of the wave induced by subaerial landslide over inclined plane to demonstrate the potential of the method to be applied to a diversity of problems. What has to be noticed here is that soils are all the way assumed to be either fully saturated or dry ones.

One-dimensional consolidation
The one-dimensional consolidation theory is used to verify the proposed method. The same verification is carried out in the literature [3,38].
The Terzaghi's analytical solution is given as where T v is the virtual time that defined as Here, the normalized parameters are defined as where k is the intrinsic permeability in Eq. (4), and E and ν are the Young's modulus and Poisson's ratio of the solid skeleton, respectively. Fig. 3(a) shows the schematic diagram of the three-dimensional numerical model. The solid column  consists of 150 (= 3 × 50 × 1) cubic cells and 600 solid material points, and each cell has a side-length of 0.02 m and contains 4 (= 4 × 1 × 1) solid material points. An isotropic linear elastic material is assumed for the solid skeleton, and the material parameters are provided in Table 1. The initial pore water pressure and effective stress in the whole domain are set at zeros, while the effect of gravity is neglected. A load of p 0 = 10 kPa is applied to the solid material points at the top layer. For both of the phases, slip boundaries are applied to the sides and bottom of the grid to realize the one-dimensional state. The total time for this simulation is 2 second, and the time increment is set at 1.0×10 −5 s throughout the calculation.

Wave induced by underwater granular collapse
We target the model experiment [28,36] of a wave generated by the collapse of a granular mass that mimics a submarine granular collapse [26]. The schematic diagram of the experimental setup is shown in Fig. 4. A soil column is delimited by a removable vertical gate, which can be rapidly removed to simulate granular collapse, and the tank is partially filled with water. For the three-dimensional numerical model, the soil column is covered by a background grid consisting of 6, 400 (= 200 × 32 × 1) cubic cells and 12,288 solid material points, and each cell has a side-length of 0.0025 m and contains 16 (= 4 × 4 × 1) solid material points. The material parameters used in this simulation are provided in Table 2. For both of the phases, the non-slip condition is considered on the bottom surface, and the slip condition is applied on the sides. The initial stresses in the soil and water are calculated by linearly increasing the gravity acceleration for 10 seconds. In this process, the soil is assumed to be elastic by giving sufficiently large cohesion. At time t = 0.0 s, a cohesion of 0 kPa is instantaneously given to all the soil material points, while the vertical gate is removed upward within 0.1s in the experiment. The total time for this simulation is 2 second, and the time increment is set at 1.0 × 10 −5 s throughout the calculation.
The obtained interface functions φ indicating the profiles of the granular mass and water surface, and the water pressures at some representative times are shown in Fig. 5. At the initial stage, the sinking of the column causes the sinking of the water surface above the column, and vortex can also be observed at the upper right corner. Then the water level rises near the wall on the left, and the waves start to propagate to the right side. At the later stages, the free surface is disturbed during the wave propagation. Despite this, the water pressure remains stable throughout the collapse, thanks to the stabilization scheme implemented in the stabilized FEM for incompressible fluids. Fig. 6 shows the simulated deposition profiles at several representative times in comparison with the experimental measurements [28,36]. The profile obtained at 1 second is slightly different from the experimental one. This is probably due to the fact that the initial condition adopted in this numerical analysis is not fully consistent with the experimental setup. Nevertheless, the subsequent profiles exhibit similar trends to those obtained experimentally. The final geometry with the sliding distance is very close to the measured one. From these results, it can be concluded that the proposed method is capable of assessing complex interactions between soil and water.
In order to further explore the capability of the proposed method, another situation is simulated: no experimental data is available for this situation. The model is almost the same as that above, with only half of the soil column soaked in water and the other half in air. The schematic diagram of the model is shown in Fig. 7. The soil column consists of 7,680 solid material points with a background grid whose cell has a length of 0.005 m and contains 16 (= 4 × 4 × 1)  Fig. 6 Comparison of deposition configurations between numerical and experimental results [28,36]: Wave induced by underwater granular collapse material points. The parameters for this numerical analysis are provided in Table 3. For both of the phases, the non-slip condition is imposed on the bottom surface, and the slip condition is applied on the sides. The total time of this granular collapse simulation is 1 second, and the time increment is set at 1.0 × 10 −4 s throughout the calculation. The results are shown in Fig. 8. As can be seen from these figures, at the initial stage (t = 0.1 s), the soil column collapses and causes the free surface rise. From 0.2 seconds on,  Throughout the simulation, the moving free surface and interface with complex geometries are properly captured without suffering numerical instabilities. In summary, it is safe to conclude that the proposed method is applicable to this type of complex interaction between soil and water.

Wave induced by subaerial slide over inclined plane
In order to demonstrate that the proposed method can be applied to a variety of problems, we target the experiment [35], in which a granular mass drops into the water to generate impulse waves, mimicking a subaerial landslide-induced wave. The schematic diagram of the experimental setup is shown in Fig. 9. The tank is 2.2 m long, 0.3 m high and includes a 45-degree inclined plane on the right side. The initial water depth is 0.148 m. The granular mass is supported on the inclined plane with a vertical removable gate, and the bottom of the mass is located at the level of the free water surface. When the vertical gate is opened, the granular mass slides down the inclined plane, impacts the water surface and generates impulse waves propagating leftward. The evolution of the generated impulse waves and the granular flow shape are recorded with a high speed camera. Four wave gauges is set at 0.45 m (G1), 0.75 m (G2), 1.05 m (G3) and 1.35 m (G4) from the vertical removable gate to measure the time-varying amplitudes of the generated waves.
For the numerical simulation by the proposed hybrid MPM-FEM, the granular mass consists of 10,368 solid material points with a background grid, each of whose cell has 16 (= 4 × 4 × 1) material points and a length of 0.004 m. The material parameters used in the simulation are presented in Table 4. It is to be noted that the inclined plane is also regarded as a solid phase and made of an elastic material to realize the non-slip and undrained conditions on the interface with the granular mass. For both of the phases, the non-slip condition is considered on the bottom surface, and the slip condition is applied on the wall sides. The total time of simulation is 2 second, and the time increment is 1.0 × 10 −5 s throughout the calculation.
The numerical results are compared with the experimental ones as shown in Fig. 10, in which the snapshots taken at three representative times are presented. As can be seen form these figures, the simulated profiles of the granular mass and wave are in good agreement with the experimental observations. At the beginning (t = 0.23 s), the irruption of the granular mass into the water causes the upthrust of the water that subsequently wraps around the irruptive mass. At t = 0.41 s, the wave crest is raised in front of the forefront of the granular slide, and then the leading wave is generated and begins to propagate away from the impact source. Also, similarly to the experimental observation, the vortex is formed around the upper part of the granular flow. In addition, the obtained sliding distance is in close agreement with the experimental one, while the forefront has not yet touched the bottom surface. After the granular slide has experienced a large deformation and turned along the bottom surface, its vortex makes a backward flow prominent at t = 0.52 s. This simulated tendency is consistent with the experimental one. However, the simulated profile of the granular mass does not agree with the experimental one completely. Accuracy improvement must be a major challenge in future studies.    Figure 11(a) shows the waveforms in the tank caused by the subaerial granular slide after the granular mass begins to deposit on the bottom, but the corresponding experimental data are not available. As can be seen from the figures, several waves are generated right after the irruption of the granular mass and propagate leftward. Among them, the first wave has the largest amplitude and mimics a tsunami induced by subaerial landslide [27]. Additionally, Fig. 11(b) demonstrates that the water pressure has been successfully stabilized during the event thanks to the stabilized FEM adopted for flow simulations of incompressible fluids. Tsunamis caused by sediment movement of granular masses are hazardous to the opposite shore [24] and therefore are first choices for our future targets to simulate. Figure 12 shows the simulated waveforms in comparison with those measured at four different points. Here, the water surface identified with the air-water interface is determined by the narrow region with interface function φ = 0.5. As can be seen from the figures, the overall responses obtained by the simulation are in a reasonable agreement with the measurements. In particular, both the peak of the first wave level and its arrival time in each of the figures closely agree with the experimental ones. However, as time passed, the simulated subsequent wave is delayed and its amplitude declines more than the experimental one. A similar tendency is also reported in [30], in which by the multi-phase flow model is employed. The main reason for the discrepancy is probably due to the low spatial resolutions for both the MPM and FEM. By refining the FE mesh and increasing the number of solid material points, we expect that our simulations would be made more accurate than the present one. This must be open to further discussion. Nonetheless, the performance of the novel method of hybrid MPM-FEM has been satisfactorily confirmed in the presented numerical examples.

Conclusion
By combining MPM and FEM, we have proposed a new numerical method for solving granular mass-water interaction problems with a view to simulating tsunamis caused by submarine and subaerial landslides. Similar to the existing solid-liquid coupled MPM, the proposed method is formulated based on porous media theory, in which saturated soil is considered as a continuum composed of a soil skeleton and a liquid phase. The Lagrangian description is adopted for the governing equations of the solid skeleton to be handled with MPM, whereas the Eulerian one is for those of the fluid phase to be solved by the stabilized FEM. Also, the phase-field method is applied to represent the free surface. In addition, MPM and FEM use the same spatially fixed grid, to which B-spline basis functions are introduced for the spatial discretization. Four representative numerical examples were presented to demonstrate the promise and capability of the proposed hybrid MPM-FEM. In summary, the following conclusion can be drawn: -Instead of the Lagrangian description, the Eulerian description can gain high precision of numerical solutions with larger time steps and greatly improves computational efficiency compared to the conventional solidliquid coupled MPMs with the help of the stabilized FEM. -The results from the numerical simulation of the model experiment of a granular collapse in a tank shows that the proposed method was able to simulate large deformation of the solid skeleton accompanied by complex solidliquid interactions and complex free surface movement. Therefore, the proposed method has a great potential in understanding the mechanism of related events such as a tsunami induced by submarine or subaerial landslides. -The results from the last two numerical examples reveal that the proposed method can also be applied to a variety of hazardous phenomena, in which the interface function drastically changes to represent the interaction among three different phases; air, water and solid.
Although the simulation results were not in a perfect agreement with the experimental measurements, the overall behavior was fairly satisfactory. It is expected that the accuracy would be improved by making the analysis conditions coincident with the experimental ones, and increasing the numbers of grid nodes and material points. It may also be worthwhile to point out that the laboratory data used to validate the proposed method were all obtained in small scale experiments. Validation analyses with reference to larger-scale real data will be conducted in the future.

Appendix A: B-spline basis function
B-spline basis function are defined with the use of a knot vector = {ξ 1 , ξ 2 , · · · , ξ n+ p , ξ n+ p+1 }, ξ i is i-th knot with i = 1, 2, · · · , n + p, n + p + 1, n is the number of basis functions, and p is the polynomial order. The first and last knot are repeated p + 1 time to make the resulting B-spline interpolatory at both end points. The Cox-de Boor formula [5] defines B-spline basis functions recursively. The zeroth order basis functions ( p = 0) are defined as follow For p ≥ 1, the basis functions are defined by the follow recursion scheme.
where the 0/0 are defined as zero. Fig.13 shows quadratic basis functions with open knot vector = {0, 0, 0, 1, 2, 3, 3, 3}. B-spline basis functions satisfy the following properties: (i) They form a partition of unity: The gradient of N i, p (ξ ) also can be computed by a recursive formula as follow The gradient corresponding to the basis functions from Fig.13 are shown in Fig.14. The B-spline basis functions of order p by a uniform knot vector are C p−1 -continuous, they guarantee at least C 0 continuity of the gradients and hence, significantly reduce the cell crossing error.

Appendix B: Phase-field method for free surface
This appendix is devoted to explaining the method to capture a moving interface between separate domain by applying the phase-field (PF) method [14].
In the PF method, the movement of a free-surface is defined as the time-variation of the interface function φ, where φ takes 0.0 for gas, 1.0 for liquid and the interme-  Fig.15.
With the values of the interface function, the density ρ and the viscosity μ of the fluid can be defined as where ρ and ρ g are the densities of liquid and gas, and μ and μ g are the corresponding viscosity.
In this study, we use the following Allen-Cahn equation [33] for the free surface where a is the gradient coefficient, W is the energy barrier height, and M φ is the PF mobility. The weak form of the interface equation using the Galerkin method can be written as About the Time integration, The approximation of the variables(time derivative, interface function) are using the similar approaches (see Eqs.(58) and (59)), and φ e in P, P s are using the current time step. The final version of time discretized form of Eq.(97) can be written as From Eq.(106), we can obtain the interface function of the next time step.