A Eulerian–Lagrangian Coupled Method for the Simulation of Submerged Granular Column Collapse

A two-fluid Eulerian–Lagrangian coupled model is developed to investigate the complex interactions between solid particles and the ambient water during the process of submerged granular column collapse. In this model, the water phase is considered to be a Newtonian fluid, whereas the granular column is modeled as an elastic–perfectly plastic material. The water flow field is calculated by the mesh-based Eulerian Finite Volume Method (FVM), with the free surface captured by the Volume-of-Fluid (VOF) technique. The large deformation of the granular material is simulated by the mesh-free, particle-based Lagrangian Smoothed Particle Hydrodynamics method (SPH). Information transfer between Eulerian nodes and Lagrangian particles is performed by the aid of the SPH interpolation function. Both dry and submerged granular column collapses are simulated with the proposed model. Experiments of the submerged cases are also conducted for comparison. Effects of dilatancy (compaction) of initially dense (loose) packing granular columns on the mixture dynamics are investigated to reveal the mechanisms of different flow regimes. Pore water pressure field and granular velocity field are in good agreement between our numerical results and experimental observations, which demonstrates the capability of the proposed Eulerian–Lagrangian coupled method in dealing with complex submerged water–granular mixture flows.


Introduction
Submerged granular column collapse is a typical water-grain interaction problem which serves as a model to investigate many underwater natural and hazardous processes, such as debris flows [1], landslides [2], submarine avalanches [3], etc. A large amount of granular material (mud, sand, gravel or rocks) that starts flowing underwater can create a tsunami [4] and damage undersea cables, pipelines, etc. A most recent example is the 2018 Sunda Strait tsunami. On 22 December 2018, a tsunami devastated several coastal regions of Sumatra and Java, Indonesia. Hundreds of people are reported dead with thousands injured. The tsunami is believed to have been caused by the collapse of Anak Krakatau that followed an eruption of the Anak Krakatau volcano in the Sunda Strait. Due to its harmful impacts to the coastal regions and underwater structures, submerged granular column collapse has long been a research concern for geophysicists, hydrologists and underwater engineers.
Ref. [5] developed a depth-averaged version of two-phase flow equations to describe the initiation of underwater granular avalanches. The rheology of the granular phase is based on a shear-rate-dependent critical state theory and a rheological model proposed for immersed granular flows. Using these phenomenological constitutive equations, the model can describe both the dilatancy effects experienced by the granular skeleton during the initial deformation and the rheology of wet granular media when the flow is fully developed.
Ref. [6] derived a set of two-phase continuum equations for studying a compressible granular flow composed of homogenous solid particles and a Newtonian ambient fluid. They found that the hydrodynamic fluid pressure helps the solids transit from densecompacted to dense-suspended granular flow, whereas the drag forces counteract the solids movement, especially within the near-wall viscous layer. Ref. [7] studied the dynamics of a dense granular flow composed of a binary mixture of small and large grains immersed in an ambient fluid. Segregation of small and large grains with the presence of ambient water was investigated.
Ref. [8] conducted the experimental research of the collapse of a granular column in a viscous liquid. The role of the initial volume fraction was investigated. They found that the morphology of the deposits is mainly controlled by the initial volume fraction of the granular mass and not by the aspect ratio of the column, an observation which differs from dry granular collapse. For initially loose packing columns, the collapse is over within the order of seconds, whereas for dense packing columns, the collapse could take as long as tens of seconds. Understanding how these complex media flow underwater remains a challenge.
Ref. [9] investigated the effect of an ambient fluid on the dynamics of collapse and spread of a granular column by means of the contact dynamics method and computational fluid dynamics. They found that the effects of fluid in viscous and fluid-inertial regimes are to both reduce the kinetic energy during collapse and enhance the flow by lubrication during spread. Hence, the runout distance of the granular material in a fluid may be below or equal to that in the absence of fluid due to compensation between these effects. They also found that for a given aspect ratio and packing fraction, the runout distance may be similar in the grain-inertial and fluid-inertial regimes but with considerably longer duration in the latter case.
Ref. [10] developed a mixture model to predict various flows involving high concentration liquid-particle mixtures. Compared with the experiments of Ref. [8], their approach gives reasonably good predictions for the case of initially loose columns, but it was unsuccessful in simulating the collapse of an initially dense columns. Ref. [11] also proposed a Eulerian-Eulerian two-phase model based on a collisional-frictional law for the granular stress to describe the underwater granular flows. However, free surface wave was not considered in their paper.
Recently, a two-fluid SPH mixture model to analyze the water-soil interactions was proposed by the authors [12][13][14][15]. With this mixture model, it is possible to investigate the temporal and spatial evolutions of the volume fractions of both constituents. Dilatancy or compaction of granular materials can also be taken into consideration in this model.
In this paper, a Eulerian-Lagrangian coupled method for the simulation of underwater granular column collapse is presented. There are some reasons for developing coupled methods. First, for the water phase, it often provides an environment for the granular flow. Using Lagrangian method to solve the environmental water flow problem is always difficult due to the time-consuming process of particle methods. In fact, large number of numerical particles are needed for the water phase and the interaction pair particles should be searched in a fast and efficient way. However, with a mesh-based method, the cell size can be much larger than the particle interval and the node connectivity remains fixed in most of the cases. Second, the water pressure term in the Lagrangian methods is difficult to calculate correctly due to the weak compressibility of the model. On the contrary, the mesh-based method can result in a very smoothed pressure field. Thus, the coupled method should be much more efficient and accurate than the particle methods, as illustrated by [16] where a novel methodology that combines the smoothed discrete particle hydrodynamics (SDPH) and the finite volume method (FVM) was developed to enhance the effective performance in solving the problems of gas-particle multiphase flows.
In this paper, the coupled framework of SDPH-FVM of [16] is extended to the water-granular mixture flow. For the water phase, FVM is employed to discretize the water flow field on a stationary grid to simulate the fluid motion. At the meanwhile, the SPH method is employed to simulate the large deformation of the granular phase. The proposed two-fluid mixture Eulerian-Lagrangian coupled model is applied to analyze the water-grain interactions during submerged granular column collapse. Experiments of the same scale are also conducted to verify the numerical results. Dilatancy or compaction effects of the granular materials on the formation of different regimes are analyzed, which helps reveal the mechanisms of granular material collapse in the presence of an ambient fluid.
This paper is organized as follows. In Section 2, the mathematical formulation of the mixture theory is given. Section 3 is devoted to the numerical implementation of the Eulerian-Lagrangian coupled method. Model validation and numerical simulation of underwater granular column collapse are presented in Section 4. Conclusions and remarks are given in Section 5.

Water-Soil Mixture Model
In the two-fluid mixture model, the grains and the fluid are described as two continuum phases characterized by individual velocities, stresses and interaction through hydrodynamic forces. For each phase we define partial densities ρ η , partial velocities v η , and partial stresses σ η , with η = l, s standing for fluid and solid, respectively. Each phase must satisfy individual balance laws for the conservation of mass and momentum where ⊗ is the dyadic product, g the gravitational acceleration, D η (·)/Dt the material time derivative along the path of particles of η phase with D η (·) Dt = ∂(·) ∂t + v η · ∇(·), f η the inter-phase force exerted on phase η by the other constituent, satisfying the Newton's third law, i.e., f l + f s = 0. Making use of Equation (1), the momentum Equation (2) can be rewritten as The above governing Equations (1) and (3) are expressed with partial variables, but the constitutive relationships are written with intrinsic (or true) variables as if they were single-phase flows. Thus, it is necessary to postulate the relationship of both kinds of variables. Here it is assumed that where variables with a tilde (˜) denote intrinsic variables, and φ η is the volume fraction of phase η, satisfying the saturation condition of φ l + φ s = 1.
For the stress, we assume that whereσ s is the intrinsic stress of the granular material, p the pore water pressure, andτ l the intrinsic deviatoric stress tensor of the water. The volume fractions φ s and φ f in (5) imply that stresses are only acted on part of each constituent due to the porosity of the granular material. The interaction force f s (i.e., − f l ) is assumed to be in the form Here the first term on the right-hand side can be identified as a buoyancy force and the second term is simply an inter-phase resistance term, with C d being the drag coefficient. C d can be calculated according to the Darcy's law [13].
In this study, the water phase is calculated in the Eulerian framework, thus it is convenient to rewrite the momentum equations of water as where use has been made of The shear stress tensorτ l for water phase is given bỹ where µ is the dynamic viscosity of the water, ν = μ ρ l the kinematic viscosity, and the strain rate tensorε αβ l is defined as in which δ αβ is the Kronecker delta symbol, and γ is a dummy index applicable to Einstein's summation convention. The liquid phase is allowed to have limited compressibility, in which case Equation (1) for water can be written as where the constant c l is the speed of sound in the fluid. For incompressible fluids 1/c 2 l is set to zero. In the limited compressibility model density changes are assumed to be small (say less than 1%) and the ρ l appearing in the pressure gradient terms in Equation (11) can be treated as constant. An artificial sound speed in water is set to a value lower than the physical one to obtain a higher time step size according to the Courant number. The staggered grid method is employed to discretize the pressure and velocity fields. A weakly formulation of the continuous equation is applied to achieve the pressure iteratively.
A finite volume method (FVM) is used to solve the equations of water phase. Underwater granular column collapse can result in free surface waves which can be captured using the volume of fluid (VOF, Hirt and Nichols [17]) method. It is assumed that near the free surface the granular phase vanishes and there exists only the water and air phases separated by the free surface. The time dependence of the VOF function F, defined as the volume fraction of the water in a computational grid cell across the water-air interface, is governed by For an incompressible fluid, Equation (12) may be combined with Equation (1) to yield As we have seen, there are two volume fraction functions: the first one is the volume fraction F to distinguish the free surface between the water and the air; the second is the volume fraction φ l denoting the water content in the liquid-granular mixture. The first volume fraction is introduced by the VOF method. The second volume fraction is introduced by the mixture theory.

Constitutive Laws for the Solid Phase
Consider the granular material as an elastic-perfectly plastic material with a Drucker-Prager yield criterion where I 1 is the first invariant of the total stress tensorσ αβ s , and J 2 is the second invariant of the deviatoric stress tensorτ αβ s , defined by For plane strain problem, α θ and k c in Equation (14) are determined by (see Bui et al. [18]) where c and θ are the cohesion and the friction angle of granular materials, respectively. A plastic flow rule is needed to describe the development of the plastic deformation. The flow rule is characterized by the plastic potential function H. If H is chosen to be equal to the yield function Y, we then use the so-called associated flow rule, with which the constitutive relation can be written as [18] whereė αβ s is the deviatoric strain rate tensorε αβ s , G the shear modulus, K the bulk modulus. Here the Jaumann rate has been applied, with the rotational rate tensorω αβ s defined aṡ The rate of change of plastic multiplierλ is calculated bẏ The second type of flow rule is the non-associated flow rule, in which the plastic potential function H takes the form where ψ is the dilatancy angle. In this case, the constitutive relation for the granular material iṡσ where the rate of change of plastic multiplierλ is given aṡ In this study, a critical state theory [5,19] is employed to determine the dilatancy angle ψ where φ eq is the granular volume fraction obtained in the steady regime, K 3 is a constant, |γ| is the magnitude of the shear rate tensor ( , defined by its second invariant, i.e., |γ| = 1/2γ ijγij . Equation (23) can be derived from the mass conservation equation and the definitions of the dilatancy angle ψ and the strain γ [20]. It implies only that dilatancy can cause variation of porosity which in turn affects the pore pressure and efficient solid stress, as well as the inter-granular friction. Equation (24) is a closure relation, obtained by fitting the experimental measurements. In this study, K 3 is set to 4.09, according to Pailha and Pouliquen [5]. Because Equations (23) and (24) act as a feedback control law for φ s , this phenomenon is known as "pore pressure feedback" [21,22] and will be shown in our numerical simulations later.
Dilatancy (or compaction) also contributes to the inter-particle friction due to the geometrical entanglement. To reflect this, the friction angle θ in (16) should be replaced by θ + ψ, see [15].

FVM and SPH Coupling
In the numerical simulation, the domain of water phase is discretized by Eulerian mesh, while the granular phase is discretized by material particles, as seen in Figure 1. Both fictitious cells and virtual boundary particles are used to implement the boundary conditions. The VOF code used is an adaptation of the SOLA-VOF code developed in the Los Alamos National Laboratory (Hirt and Nichols [17]). An in-house SPH code is developed to undertake the analysis of the granular phase, for the problem of submerged granular column collapse. In this section, the VOF and SPH implementations are not described in detail. Interested readers are referred to [17] for VOF method and [12,13] for SPH implementation. The key issue is the transfer of the information between Eulerian grid nodes and Lagrangian particles. The interpolation process is depicted in Figure 2. In our approach, this transfer is accomplished using the smooth kernel of the SPH method. In fact, the foundation of SPH method is based on the following integral representation of a field function f (x): where W(x − x , h) is called the kernel or smoothing function, and h is the smoothing length which defines the influence domain Ω of the kernel W(x − x , h). In SPH method, the computational domain is discretized into a finite number of particles. The continuous integral representation of the field variable f (x) at position x i , in Equation (25), is then approximated by summation over the neighboring particles in the support domain, as where W ij = W(x i − x j , h); m j and ρ j are the mass and density of particle j at position x j , respectively; N is the number of particles in the support domain of particle i. A Wendland quintic kernel function is chosen in this study to perform the SPH interpolation approximation (Wang et al. [12]). SPH method suffers from the so-called boundary deficiency problem due to the truncation of the support domain of the kernel function W ij by the boundary. To overcome this deficiency, a re-normalization technique is introduced to perform the SPH approximation (Chen et al. [16]): where N is the number of particles in the support domain of particle i. With this renormalization technique, information transfer from the solid particles to the grid nodes can be implemented by where x g is the node position, N p is the number of solid particles in the support domain of grid node x g . Analogously, information transfer from the grid nodes to a certain particle can be done by where N g is the number of grid nodes in the support domain of the solid particle at x i ; m g and ρ g are, in principle, the mass and density of the water particles at grid node x g , respectively. However, due to the Eulerian nature of the FVM method for the water phase, no water particles exist in the present method. In fact, m g /ρ g represents a volume of fluid at the position x g , thus the volume of the cell centered at x g can be used as m g /ρ g to perform the interpolation in (29). Usually, the cell size is larger than the particle interval. The support domain of the interpolation function should therefore be enlarged to contain sufficiently many grid nodes. Interpolation strategy between Eulerian grid nodes and Lagrangian particles. Staggered grid for the Eulerian variables is illustrated. Left: Grid information transfer to a certain particle; Right: Transfer of particles information to the grid nodes. A square-type cell system (same ∆x, ∆y sides) is used here. kh is the radius of the support domain of a solid particle or a grid node, h is the initial distance between solid particles and k is a constant normally in the range from 2 to 4. A larger value of k implies more particles in the support domain.

Special Treatment of the Pore Pressure
The interaction force f s in (6) contains the buoyancy force −φ s ∇p. Using the SPH approximation rule (26), this term for the solid particle i would be approximated by where subscript g means the nearby grid nodes. However, this approximation results often in numerical instability (shown later). Thus, the buoyancy force −φ s ∇p at solid particle i is replaced by where p i is the pore water pressure at the location of solid particle i, obtained using the re-normalized SPH interpolation approximation (29). Water volume fraction φ g on a grid node x g is calculated, according to the saturation condition φ l + φ s = 1, by where the summation term gives solid volume fraction φ s at the position of grid node x g .

Time Integration and Boundary Conditions
For the FVM simulation, a forward Eulerian integration scheme is employed to integrate Equations (7), (11) and (13) for the fluid phase, while for the SPH calculation, a second-order accurate numerical integration scheme, namely the Leap-Frog (LF) method is adopted here to integrate Equations (1), (3), (17) or (21), and (23) for the granular phase (see Wang et al. [14]). The time step size is controlled by the so-called Courant-Fredrich-Levy (CFL) condition and additional constraints due to the interaction forces and viscous diffusion where ∆t is the time step size; h is the smoothing length of particles; c l and c s are the sound speed of water and grains, respectively; f l and f s are the drag forces acting on water particles and solid particles, respectively; ν is the water viscosity. The time step ∆t for a simulation should be chosen to be the minimum value of the above constraints. Although Leap-Frog integration is a second-order method, in contrast to Eulerian integration scheme, which is only first order, the accuracy of the time stepping scheme of the coupled method is only first order.
In FVM simulation, in addition to the free surface boundary conditions treated by the VOF technique, it is necessary to set conditions at all mesh boundaries. At the mesh boundaries, a variety of conditions may be set using layers of fictitious cells surrounding the mesh. The basic idea is to set values for the dependent variables in the fictitious cells such that the desired boundary conditions are met at the boundaries, see [17]. For the SPH simulation, the wall boundary is treated by the dynamic particle method, proposed by [23].

Dry Granular Column Collapse
We first apply the proposed SPH model to the problem of single-phase dry granular column collapse in the air. This problem has been studied by the present authors [12], using SPH method. In [12], we compared our numerical results with experimental and computational results of [18]. An excellent agreement was observed. Thus, in this section, only the effects of some material and geometrical parameters, such as the internal friction angle θ, the material stiffness E and the aspect ratio a of the column, on the granular column collapse, are investigated.
We first investigate the case that a 2D granular column of 0.1 m wide and 0.1 m high is released at t = 0. The initial configuration of the granular column is shown in Figure 3. Material properties for the granular medium are listed in Table 1. The granular material is assumed to be cohesionless. In the simulation, a total of 12,800 solid particles are used, with an initial particle spacing of 0.00125 m. The sound speed is chosen c s = 150 m/s, which is needed to control the time step size according to the CFL condition (33). The time step size is set to ∆t = 5 × 10 −6 s.  In this study, computations of the plastic behavior are based on the non-associated flow rule which results in the phenomenon of dilatancy. For the non-associated flow rule, the significance of dilatancy can be characterized by the dilatancy angle ψ. In the simulation for the dry granular collapse, the dilatancy angle ψ is set to be constant, while for submerged cases, this variable is a dynamic quantity computed from the mixture rheology relationships (23) and (24), as seen in the next section, i.e., although in the dry case, a granular medium obeys a relaxation towards the critical state, we do not use the same dilatancy law as used in the submerged section with the dilatancy angle coupled with the volume fraction. We point out that although dilatancy can adjust the friction coefficient to some extends, its effect on the collapse dynamics of a granular medium is prominent only if it is interactive with the pore pressure feedback mechanism. Here for the dry case, a dilatancy angle ψ = 2.5 • is taken, which normally results in a relatively high dilatancy. Figure 4 shows the influence of the internal friction angle θ on the morphology of the deposit, where two friction angles, namely θ 1 = 20 • and θ 2 = 30 • , are considered, respectively. We can see that the runout distance of the granular column decreases with the increase of the internal friction angle. This is reasonable if we bear in mind the fact that the internal friction increases the strength of the granular material, as seen in Equation (14).    The influence of the aspect ratio a of the column on the morphology of the deposit is shown in Figure 6. The aspect ratio a is defined as a = H i /L i , where H i and L i denote the initial height and width of the granular column, respectively. We can see that for a granular column with high aspect ratio, say a > 1.0 for example, the height of the deposit is lower than the initial height of the column, and at the end of the collapse, a pile with a sharp top is developed. However, in the case of low aspect ratio, a < 1.0, the height of the column is almost unchanged during the collapse, leading to a pile of flat top at the end of the collapse. This result agrees with the experimental observations conducted by [24]. Ref. [18] pointed out that the use of the continuity Equation (1) in the current SPH application to dry granular material is optional and can be removed from the governing equations if the bulk density is kept constant. Retaining the continuity equation with the assumption of constant density corresponds to resolving the variation of void ratio or porosity of a granular material. The idea is illustrated in Figure 7, where numerical simulations with and without continuity Equation (1) are presented. We can see that the numerical results in these two cases are nearly the same, thus the continuity Equation (1) is not important for the simulation of dry soil column collapse. However, in the case of submerged granular column collapse, the continuity equation should be employed to resolve the variation of volume fraction of granular material, as seen in the next section.

Submerged Granular Column Collapse
In this section, the submerged granular column collapse is simulated using the proposed coupled Eulerian-Lagrangian model for two-fluid mixture flows. Both initially dense and loose packing columns are considered in succession. Figure 8 shows a rectangular tank of 0.50 m in length and 0.15 m in height. The initial water level in the tank is 0.1 m. The submerged granular column is L i long and H i high.
At t = 0, the gate is suddenly removed, and the column then collapses. The gate can be removed within 0.03 s. Although granular column collapse is known to be quite sensitive to the way the gate is removed, the effects of the gate movement on the granular motion and the surrounding fluid motion are ignored. Material properties used in the calculation are shown in Table 2. Model experiments of the same scale are also conducted to validate our numerical simulations. The experimental setup is the same as that in [13], but with a pressure sensor mounted on the side wall of the tank to measure the pore water pressure at point P 1 . The position of the pressure sensor is shown in Figure 8. The sensor used is an HM91 micro pore pressure sensor and transducer produced by HELM SENSOR. The measuring range of this diffused silicon piezoresistive pressure sensor is from 0 to 10 kPa, with a 24 V DC power supply. The overall accuracy is ±0.1%FS (Full Scale). The loose packing is made by gently pouring the spherical glass beads into the reservoir delimited by the removable gate. To create dense packing columns, we gently tap on the tank wall. By doing so, the solid volume fraction reached in our setup is 0.55 ± 0.005 for the immersed loose packing and 0.60 ± 0.003 for the immersed dense packing. The solid volume fraction is measured and calculated by the displacement method. The diameter of the spherical glass beads ranges from 150 µm to 300 µm, with an average diameter of 225 µm and standard deviation of 40 µm.  There are a total of 9600 Eulerian cells and 3072 SPH particles for water and solid phase in the tank, respectively. The initial particle spacing is 0.00125 m. The grid size is twice as that of the solid particle interval, which means there are nearly four solid particles in one cell. Thus, the support domain of a solid particle should be doubled in order to contain more grid nodes. Numerical parameters are chosen as follows: Virtual speed of sound for water c f = 10 m/s; speed of sound for granular material c s = 150 m/s; artificial viscosity parameters α = 0.1, β = 0; the artificial diffusive parameter in δ-SPH method δ = 0.1 for SPH simulation of solid phase; for the repulsive force, D = 0.01 and r 0 = 0.00125 m; time step size ∆t = 2.5 × 10 −6 s, which satisfies the constraints in (33). The non-associated flow rule with dilatancy angle ψ determined by (23) and (24) is implemented in the calculation. For the meaning of some of the above parameters, see [12].

Dense Packing Column Collapse
Submerged granular column collapse can cause free surface waves. Figure 9 shows the free surface indicated by the water volume fraction φ f , where φ f = 0 in the air domain, φ f = 1 in the pure water domain and 0 < φ f < 1 in the water-grain mixture. Here it should be pointed out that although the VOF function F is used to indicate the water free surface in the Eulerian domain, we combine F and φ f together and show the water free surface and solid surface simultaneously using φ f only, as seen in Figure 9. It is seen that the collapse of the submerged granular column leads to the sinkage of the water free surface above. The magnitude of the sinkage of the water free surface is in the order of 1cm, which agrees well with the experimental observations. Thus, the proposed numerical model has the capability of simulating free surface flows caused by underwater granular column collapse.  The granular column profile indicated by the water volume fraction φ f is also shown in Figure 9 which resolves the water volume fraction within the column and at its surface. Here, it should be pointed out that the granular surface is, in fact, a surface where φ f changes rapidly across it. In our numerical model, the water volume fraction φ f is calculated by interpolating the water volume fraction φ a (i.e., Equation (32)) on the Lagrangian solid particles to the Eulerian grid points using the SPH re-normalized interpolation approach (i.e., Equation (27)). Usually, in SPH simulation, a transition zone of φ f at the granular surface develops, due to the incompleteness of the support domain of the SPH interpolation function at the interface. However, as seen in Figure 9, due to the re-normalized interpolation technique used in this study, the transition zone is limited in a narrow space, namely the abrupt change of the volume fraction φ f at the interface is maintained.
We investigate the effectiveness of the special treatment of the pore water pressure mentioned in (31). Figure 10 shows some solid particles expelled out from the granular column by the improper treatment of the pore water pressure, i.e., Equation (30). Particle expelling results in the enlargement of the particle intervals. Consequently, a large transition zone of the water volume fraction φ f is formed and the pore water pressure in this zone cannot be predicted correctly. This instability, previously reported by [25], can be avoided by introducing p in the approximation of the pore pressure gradient, as seen in (31). The main idea behind this treatment is to address the importance of the pressure difference, instead of the pressure itself, in the approximation of the pore water pressure gradient. With this treatment, the movement of the solid particles and the evolution of the water volume fraction φ f are stable, as seen in Figure 9.  Figure 11 shows the snapshots from both experimental observations and numerical simulations of an initially dense packing submerged granular column collapse. The granular velocity fields at representative times from the experiments and simulations are given for comparison. The experimental velocity fields of the solid phase are obtained by the software package Photoinfor which can capture the deformation patterns and analyzing the evolution of the deformation quantitatively (Li et al. [26]) by means of the so-called digital speckle correlation method (DSCM). The snapshots on the left column of Figure 11 are from experimental observations with arrows representing the velocity vectors of the granular particles they start on at a certain instant. It is seen that upon removal of the gate, particle movement takes place mainly at the upper left corner where the particles fall off the pile with a large velocity. After that, the particles mainly move along the granular surface layer and a "hydraulic-like granular jump" can be observed due to the stoppage of the slow-moving particles at downstream. The same phenomena are observed from the numerical simulations, except the discrepancy on the runout distance which may be due to the sparseness of particles at the surge front. Both experimental observation and numerical simulation show that the collapse of the initially dense packing column is finished within about 3 s. Figure 12 shows the water pressure field at some representative times for an initially dense packing granular column collapse. Dense packing granular material exhibits dilation when sheared. As a result of the dilation, the interstitial pore water is sucked into the column, giving rise to a negative pore pressure field. In Figure 12, the hydrostatic pressure has been subtracted and only the dynamic deviation of the pore pressure is shown. That would make the low-pressure zone much more obvious. It is seen from Figure 12 that a lower pressure zone in the shearing zone develops at the initial stage of the dense column collapse. During the later stage of the collapse, this negative pressure field dissipates and recovers to the hydrostatic pressure due to the pore water infiltration. To investigate the time evolution of the pore water pressure quantitatively, pore water pressure at point P 1 is investigated experimentally and numerically. Figure 13 shows the predicted relative pore water pressure at point P 1 compared with the experimental measurements. Here, "relative" means the pressure difference from the hydrostatic pressure of the surrounding fluid at the same depth. Experimental measurements show that at the initial stage, say 0 < t < 1, a negative pressure is measured, and then it is recovered to the hydrostatic pressure during the later stage of the collapse. The fluctuations in the experimental pore pressure signal are due to the sloshing of the water free surface in the tank caused by the removing of the gate. Numerical pore water pressure shows the same trend, but the negative pressure appears later than in the experiments. The reason may be also due to the gate removing duration.

Loose Packing Column Collapse
Now we investigate the loose packing case. The initially loose packing column collapse shows some different scenarios. In Figure 14, the granular velocity fields obtained for an initially loose packing column is shown. It is seen that upon removal of the gate, the upper part of the column drops quickly, resulting in a surge of the front. After that, a thin layer of surface particles flows down a mild slope. Both experimental observation and numerical simulation show that the collapse of the initially dense packing column of this case is finished within about 1.5 s, much shorter than in the case of dense packing collapse. The experiments of [8] revealed that the collapse of an initially loose granular column can be finished within a few seconds, while for an initially dense granular column, the collapse can last for a much longer time, e.g., tens of seconds. Moreover, the runout distance of loose columns is twice longer than that of dense columns. Please note that [8] used a kind of very sticky liquid (µ = 0.012 Pa·s) in the experiments. Thus, the difference between our results and the experimental observations of [8] may be due to the fluid viscosity, which should be a good topic for the future. Figure 15 shows the relative pore pressure field evolving during the collapse of an initially loose packing granular column. It is seen that due to the compaction of the loose packing, a high-pressure zone develops in the granular column at the initial stage of the collapse. For loose column collapse, a positive pore pressure can be measured below the column. Here, positive pressure is relative to the hydrostatic pressure of the surrounding fluid at the same depth. After that, the pore pressure dissipates and restores to the hydrostatic pressure in a short period. Figure 16 shows the predicted relative pore water pressure at point P 1 compared with the experimental measurements. Experimental measurements show that at the initial stage, say about 0 < t < 0.2 s, a positive pressure is observed, and then it is recovered to the hydrostatic pressure during the later stage of the collapse. Numerical pore water pressure shows the same trend, but the positive pressure lasts a longer time than in the experiments. Interestingly, overshoot and oscillation of the pore pressure can be seen from both numerical and experimental results. Please note that in the dilatancy model (24), K 3 is a rheology parameter which should be calibrated in practice. Thus, a well calibrated K 3 might be useful to obtain a more realistic pore pressure feedback mechanism. However, this should be clarified in the future.

Conclusions
In this paper, numerical simulation of submerged granular column collapse by means of a Eulerian-Lagrangian coupled method is presented. Conservation equations of mass and momentum for both phases (water and solid particles) are formulated, based on the mixture theory. The solid is considered to be an elastic-perfectly plastic material. A Drucker-Prager yielding function with associated and non-associated plastic flow rules, respectively, is employed to describe the elastic-plastic behavior of the solid particles. The critical state theory is combined to consider the dilatancy or compaction of the granular material when subjected to shear.
A Eulerian-Lagrangian coupled method is employed to solve the proposed two-fluid mixture model. In this method, the continuous phase is solved using the Eulerian finite volume method (FVM) with the free surface captured by the volume-of-fluid (VOF) technique. The dispersed phase is solved by the Lagrangian smoothed particle hydrodynamic method (SPH). The Eulerian grid is composed of rectangular uniform Cartesian cells. The granular is decomposed into many material particles. The key issue is the information transfer between the Eulerian grid nodes and the Lagrangian particles. In the present method, this transfer is accomplished using the SPH kernel function. A re-normalized interpolation technique and an alternative treatment of the pressure term in the momentum equation for solid are employed to reduce the numerical instability of the SPH simulation.
The dry granular column collapse is studied first using the proposed approach. We investigate the influence of the stiffness, internal friction angle and aspect ratio on the morphology of the deposit. It is found that the morphology of the deposit is mainly dependent on the internal friction angle and the aspect ratio. The stiffness of the solid particles, which is identified by the bulk modulus of the solid material, has little effects on the granular deposit.
Compared with the dry granular column collapse, submerged granular column collapse exhibits many different features. For loose packing, the collapse process takes within a second; whereas dense packing can retard the flow significantly. Dilatancy (compaction) effects on the formation of different regimes is investigated. Agreement between our numerical simulations and experimental observations is good.In many situations, such as the underwater landslides, it is necessary to consider the effects of the ambient water on the mixture movements. The proposed coupled approach in this paper can deal with such problems widely occurred in underwater engineering and science.
Author Contributions: Conceptualization, C.W., X.M. and C.P.; writing-Original draft preparation, C.W.; writing-review and editing, Y.W.; coding and data analysis, C.W.; software for experimental data analysis, G.Y. All authors have read and agreed to the published version of the manuscript.