Granular Material Point Method : unsaturated soil modelling

The paper proposed an enhancement to the Granular Material Point Method extending the original framework for unsaturated granular materials. It discusses briefly the key processes occurring during unsaturated flows and gives a set of assumptions for simplified analysis of such flows with the Granular Material Point Method. Then it proposes an algorithm for the Granular Material Point Method for unsaturated materials and shows how to adjust the original Granular Material Point Method code. Finally, the paper presents two sets of simulations involving unsaturated flows comparing outcomes from the previously proposed Granular Generalized Interpolation Material Point Method (Granular GIMP), original GIMP, as well as two newly proposed Granular GIMP simulations with an enhancement for unsaturated soils. The results show that the proposed formulation leads to improved and qualitatively different results and confirm that considering partial saturation is essential, also when the flow of material is of interest. The paper uniquely illustrates the developments not only with figures but also with videos, giving the whole extent of simulation instead of just a timestamped image.


Introduction
Granular flows occur both naturally (e.g. in landslides) and as a result of human activity (e.g.mixing and transport of granular materials in factories).The behaviour of granular flows is difficult to predict and replicate, even with modern numerical methods.Yet, having such a numerical method suitable for estimating landslide outreach or impact forces when the flow hits obstacles, man-made barriers and other structures would immensely help in safe design and construction.Similarly, a more accurate estimation of granular flows in the industry would help in the design of silos, pipelines and mixers.
The majority of granular flows are initiated from a quasistatic state.In that state, the behaviour of granular material is well replicated by constitutive models that rely on continuum mechanics, as in the quasi-static state, before failure, the assumption of continuity holds.These constitutive models average the behaviour of individual grains into well-defined macroscopic properties of granular material, allowing for parametrisation of the model based on simple laboratory tests.
As the continuum-based constitutive models and associated numerical methods are widely used in engineering, it is attractive to extend those to situations where the material behaviour changes and the initially neighbouring particles move away and are replaced by other particles, or when cracks in the material occur, invalidating the continuum assumption.For example, the continuum may be enriched with discontinuities in displacements (open/closed cracks) or discontinuities in stress/strain fields (e.g.Ref. 1).Granular Material Point Method (GMPM), Seyedan and Sołowski, 2 takes such concepts further to approximate the behaviour of the material during the flow.
Flows of unsaturated granular materials are difficult to model numerically as the flowing material behaviour is complex (e.g.Refs.3, 4) and involves large displacements and strains.GMPM replicates three regimes of dry granular material behaviour: quasi-static, moderate and disconnected flow, and disregards the collisional flow, which is uncommon and usually observed only in laboratories.In the quasi-static regime, GMPM assumes that material behaviour is following that established in geotechnical engineering, described by one of many constitutive models for granular materials.However, already in this regime relatively rapid deformation may take place.The material only transitions to the moderate flow regime when its density reaches the critical density, defined as the density at which volumetric deformations stop upon shearing.During moderate flow, the grains remain densely packed, volume changes are small and shear deformations dominate, creating a dependency between the shear rate and the shear stress.The GMPM framework further assumes a transition to disconnected flow, when granular material reaches the critical density, at which the grains lose contact with each other.As grains interactions are very infrequent in the disconnected flow, the stress of a material point is negligibly low and therefore assumed to be zero in GMPM, similar to the stress-free assumption of Dunatunga and Kamrin. 5GMPM also disregards any volume change in a disconnected state as the framework assumes that the material needs to be at critical density (or higher) to resume interaction.GMPM approximation of the material behaviour in the state of moderate flow or disconnected flow may be seen as of limited accuracy.However, the GMPM algorithm improves the accuracy of the simulations significantly helping them to replicate granular material behaviour much better than other existing Material Point Method algorithms.
This paper shows an enhancement of the GMPM framework so it is suitable for partially saturated granular materials while maintaining the simplicity of the original GMPM.In particular, the proposed enhancement could be suitable for simulations of flows of large masses of grains, for which the macroscopic constitutive material models based on continuum mechanics replicate the behaviour of the granular mass well.Also, for such large masses of grains simulations modelling individual grains of different shapes, contact between them and additional forces related to partial saturation due to water bridges between the grains (e.g.see Refs.6-9) are impractical, if not impossible.
The primary novelty of the paper is the Granular MPM framework enhancement for unsaturated soils described in Section 2. Furthermore, for qualitative validation of the framework, we show an option for adjusting the constitutive models close to the critical material density that qualitatively improve the replication of the material behaviour.Finally, the paper shows two simulations that illustrate how the proposed improvements help in the replication of unsaturated granular flows.
GMPM is the chosen approach for simulations, as during flow, granular material experiences very large strains and deformations.These large deformations are difficult to model with classical numerical methods based on continuum mechanics, such as the Finite Element Method (FEM), where the domain is discretised with a mesh.Even when using advanced versions of FEM, large deformations cause excessive deformations to the elements, leading to inaccuracies and stability issues.This can be remedied by re-meshing, but that leads to errors, due to the projection of variables from the old, deformed mesh to a new one.Such issues led to the development of numerical methods being not dependent on the initial mesh, such as the Material Point Method, with further enhancement (Generalized Interpolation Material Point Method, GIMP, Bardenhagen and Kober 10 ), used in this study.The Material Point Method does not discretise the continuum with a mesh of elements, but material points instead, which carry all the important material variables.The calculations in the Material Point Method are performed on a fixed grid.That means in each step the material data from the material points must be transferred to the grid, where the calculations are made in each step, and then transferred back to the points.

Overview of granular flows
Unsaturated granular flows are not that different to the flows of dry granular materials.However, the grains of unsaturated material can stick to each other more easily, thanks to suction induced intergranular forces.At a low degree of saturation, in the pendular state, grains are connected with water bridges, resulting in net force keeping the grains together. 11The net compressive force between the grains persists at higher degrees of saturation when the granular material enters funicular state with water bridges connecting multiple grains, [12][13][14] and in the capillary regime at a yet higher degree of saturation, to become zero at full saturation.In the quasi-static state, these extra forces between grains are taken into account in the continuum based constitutive models for unsaturated soils.Those models alter the granular material behaviour in line with the degree of saturation, or, through the means of a water retention curve, with suction.Therefore, to capture the unsaturated granular material behaviour in the quasi-static regime, GMPM only needs to use the correct constitutive model or models, if needed coupled with a water retention curve and laws capturing the transport of water in the material.
The moderate state of flow occurs in the dry granular material when the grains lose persistent contact with the neighbours (as the neighbouring grains frequently change), and the grains do not, on average, push each other away.As such, the moderate flow remains dense and occurs at an almost constant volume.Sometimes, however, the material point in the moderate flow state may end up with very low mean stress, leading to its volume above the critical and its density below the critical density.In a real flow situation, grains in the material volume corresponding to such a material point would lose contact with each other due to being above the critical volume.Hence, when such a situation occurs, GMPM assumes that the given material point has no volume change, limited shear deformations and zero stress.
Otherwise, in this paper, we do not use a special constitutive model for the moderate flow regime, hence we do not model the transition to this state, but treat moderate flow as an extension of quasi-static behaviour.However, we do propose cohesion reduction at a specific volume close to critical, showing mechanics that could be used to capture some of the changes in the material behaviour in moderate flow.
In unsaturated soils, the packaging of the grains is somewhat affected by the suction induced forces.However, those forces are very small in granular materials with grain sizes equal to or larger than those of medium sand (diameter above 0.25 mm), when compared to e.g. grain weight.Therefore, it is unclear whether the nature of the moderate stage of flow for unsaturated granular material is much affected.In the GMPM framework, it is assumed that the constitutive model for the material will capture the changes in the granular material behaviour depending on the suction and volume in the quasi-static state.
In the disconnected flow phase, due to suction induced intergranular forces, larger chunks of material are more likely to stay together.However, the grains will still not interact with each other, but rather just follow the kinematic trajectory.Hence, the granular material during disconnected flow will not interact with the other pieces of material and remain at a constant volume.As such, the assumptions valid for the disconnected flows seem to stay valid for the unsaturated disconnected flow.
During the flow, the particles of granular material will segregate.That means that the initially uniform granular material will partially separate into fractions, with grains of certain sizes likely to migrate towards the top or bottom of the flow.This is a complex process and difficult to capture, both in dry and partially saturated flows see Ref. 15. Furthermore, initially different granular materials will mix, leading to new granular materials with different properties (e.g.Ref. 16).The mixing processes are complex and framework.For those cases the Discrete Element Method simulations may be preferable, as long as the amount of particles is manageable (e.g.Refs.8, 9).
The unsaturated granular flows include an additional component-water.Therefore, they need a special set of assumptions related to water transport.Water transport occurs with water both in the liquid phase and vapour phase.In the quasistatic state, the water transport typically is approximated by a Darcian flow, which assumes a steady-state condition.When the granular flow starts, the material flows for seconds or minutes, which is a very short time when compared to the material permeability.Thus, the contribution of Darcian flow to the transport of water may be neglected.However, in granular materials with large grain sizes and close to full saturation, Darcian water transport may be significant enough to require its consideration.
Nonetheless, when partially saturated granular material flows, other processes related to water transport take place.The granular material may mix and segregate, with grains largely transporting the water attached to them.Furthermore, rapid acceleration may separate the water and solid phase.Finally, the material changes volume.This is a process neutral to water content but affecting suction value, hence material behaviour.

Assumptions of the granular material point method for simplified analysis of unsaturated granular flows
The paper presents a simplified analysis of granular flows, leading to a much simpler MPM algorithm and greatly reduced computational time.In the simulations, we assumed that: 1. Darcian flows of water in the granular material are insignificant.2. During flow, both water and grains are moving together at the same velocity.The relative velocity of the solid phase and water phase is neglected.3.In reality, during moderate flow, in the case of volume computed to be above the critical volume and hence assumed zero stress, the volume of the material can change.That, in the absence of a water source or sink, would lead to a change in the degree of saturation of a given material point.However, we neglect any such volume changes in calculations.
4. In reality, during moderate flow, in the case of volume computed to be above the critical volume and hence zero stress, the material may still have some shear strength.We fully neglect this shear strength of the material in the calculations.5.The water transport due to the mixing of material during flow is negligible and neglected.6.The material will not become fully saturated.Hence, we neglect the possibility of the development of positive pore pressure.7. Segregation of material can be neglected 8. Mixing of different materials leading to the creation of new materials with different properties is not taken into account.
As the framework uses an explicit Material Point Method formulation, the timeframe of the simulation is usually limited to seconds, or, at most, minutes.Therefore assumption 1 is typically not very limiting.Yet, Darcian flows should be accounted for in simulations of very coarse, almost fully saturated materials in which a significant amount of flow can happen in seconds.For such materials, the simulation should also take into account the possibility of the material reaching full saturation and developing positive pore pressures (see assumption 6).
Assumption 2 likely leads to some inaccuracies in most simulations.Simulated flows do experience very high accelerations occasionally, for example when the flowing material hits an obstacle or when it is deposited.In those cases, the water phase may move respectively to the granular phase, leading to a change in the degree of saturation of the material.The simulations neglect this effect, mainly due to the computational effort required to simulate the binding between the solid phase and liquid phase, as well as the acceleration of the phases correctly.
Assumptions 3 and 4 directly follow the framework of GMPM.Currently, water retention models, even when taking into account the changes in the void ratio (e.g.Refs.17, 18), are not suitable for conditions in which the material is below the critical density and the grains of the material are losing contact with each other.As such, when the grains are not in contact and cannot create force chains, which is the condition occurring during flow, the shear strength of the material at low shear rates is zero, as stated by assumption 4. Yet, thanks to the water bridges, in unsaturated flow some viscosity or cohesion during the flow may be present and is affected even by small volume changes.However, quantifying this effect is very difficult.
There is no doubt that water transport is happening during granular material mixing, perhaps primarily during shearing.We neglect that transport in assumption 4, as the rate of exchange of water content between the materials during shearing requires an introduction of a new constitutive law that is not easy to establish.Preliminary work on a GMPM algorithm that would allow for such transport also indicated rather high computational effort associated with identifying the material points between which the water content should be exchanged.
Finally, the framework allows for disconnected flow.The mechanics of such flow is different to that of completely dry material.Instead of a cloud of dry disconnected grains, interacting with each other only rarely, the disconnected unsaturated flow may consist of larger aggregates of grains connected by water bridges and held together due to suction induced forces.Although these partially saturated aggregates of grains may experience limited volume change due to suction-induced intergranular forces (dependent on the suction value and characteristics of the grains), GMPM neglects these volume changes and keeps the volume of these material points constant.

Formulation of granular MPM for unsaturated granular flows
The GMPM framework, as given in Ref. 2 enriches any explicit MPM scheme (such as GIMP) for solids with the ability to model moderate and disconnected flows as well as transitions between the quasi-static, moderate flow and disconnected flow states.These three states can co-exist in a computational domain, even within a single material domain.
The proposed set of assumptions allows for a simple formulation of GMPM for unsaturated granular flows.Most importantly, the water phase can be treated as attached to the solid phase, with the exchange of water between material points neglected.Hence, the formulation needs only minimal adjustments in the GMPM framework for dry granular materials.
The explicit Material Point Method assumes that the material points, with domains described by the shape functions, discretise the problem over a computational grid.MPM algorithm performs steps 1-11 in each time step (here Update Stress Last is assumed), on a single computational grid.GMPM performs those steps for particles associated with each parallel computational grid and adds extra steps G1-G5, while the algorithm for unsaturated granular materials adds a further extra step U1: G1. Determine the interacting material points.
In GMPM each material point has a domain of interaction.The shape of the domain of interaction is the same as the current shape of the material point domain (or square in the original MPM), while the size of the domain of interaction is corresponding to the critical density of the material point.The critical density is defined as the density at the largest volume the material point may have with grains still being in contact.
Interacting material points are points with overlapping or touching domains of interaction.The clusters of interacting points are assigned to the computational grids in step G2.The interacting material points are in a quasi-static state or moderate flow state.In the simulation, we use the same constitutive model for both states, hence the states are not differentiated.Such differentiation could be done e.g. by defining the volume at which the material should switch to the moderate flow constitutive relationship.
In case the domain of interaction of a given point is reaching a boundary, it is also treated as an interacting material point.
G2. Assign clusters of interacting material points to the same computational grids.Add/remove grids as necessary.
Run MPM algorithm (steps 1-11) for each grid: 1. Compute grid mass from material point masses 2. Compute grid momentum from material point masses and velocities 3. Compute internal forces from material point stresses 4. Compute external forces (e.g.body forces due to gravity and tractions) 5. Compute the rate of momentum at grid nodes 6. Update grid momentum from the rate of momentum at grid nodes 7. Update material point velocities from the rate of momentum at grid nodes 8. Update material point positions from the updated grid momentum 9. Compute grid velocity from grid mass and momentum 10.Compute material point deformation gradient and strains from grid velocity 11.Update material point stress from material point strains.
G3. Identify material points states: quasi-static, moderate flow or disconnected flow.Set the deformation gradient increment to zero or at least the volumetric strain increment to zero for the material points if the density after the timestep is below the critical value.The velocity of material points and their updated position are not affected.Set stress to zero for the particles in the moderate flow state and the disconnected flow state that had their density below the critical value The identification is first based on material points density/volume.Material points denser than critical are assumed to use the quasi-static constitutive law in the next timestep.When the material point density is lower than the critical density (plus some tolerance), the stress in the material is set to zero and the deformation gradient is reset to the one from the previous step.In the next timestep, those points are counted as being in a moderate flow or disconnected flow state.The material points with no domains of interaction overlapping (being at that stage also the actual points domains, as the points are at critical density) are considered to be in the disconnected flow stage, while other points will belong to grids of interacting points, some of which may be in quasi-static or moderate flow with non-zero stress condition (the calculations assume a single constitutive law for both states, hence there is no clear distinction made).

U1. Update suction values due to volume change
Steps 1-11 are generally the same for all explicit versions of MPM.However, more advanced formulations of MPM, e.g.GIMP and Convected Particle Domain Interpolation Material Point Method (CPDI), Sadeghirad et al. 19,20 will use more complex functions or algorithms to transfer the data from the material points to the grid nodes and back.Some other algorithms may add extra steps or slightly modify the presented steps (e.g. by moving the strain and stress update earlier in the algorithm).
In GMPM, the material points can behave differently depending on the state they are in (quasi-static, moderate flow or disconnected flow).The state of the material points is determined in step G3.As we do not have different constitutive laws for quasi-static and moderate flow stages, there is no density or other transition point defined, hence the material points in those states are treated in the same way.However, if the density is below the critical density, the material point is assumed to be stress-free.It is in a disconnected flow state when the material point has no neighbours within its domain of interaction and a moderate flow otherwise.
In unsaturated flow, the critical density depends on the degree of saturation and suction.Physically, it is the density at which the material particles lose contact with others.In unsaturated soils, neglecting the connections through extended menisci, this density is the same as for fully dry soils.As such, at first approximation, a single value of critical density is enough.However, reaching that density value in unsaturated material requires overcoming the suction-induced cohesive forces.Those forces are complex and computing them would require a carefully shaped constitutive law, realistically predicting soil behaviour at low stress and perhaps allowing for a limited amount of tension.Taking this into account step G3 stays the same in case of partially saturated flows.
Even though in the simplified algorithm presented the water content of material points remains constant, the value of suction can change when the volume of particles changes.As such, it is necessary to update the suction value in each round of the algorithm, as done in step U1.This requires an extra constitutive assumption-a water retention curve.
The algorithm can be further enriched by water transport, including Darcian flow, transport due to acceleration and material mixing.We suggest that the water transport could be performed as an extra algorithm in step G3, before resetting the deformations of the selected particles of the material in moderate flow as at that stage the material mixing could be quite easily taken into account.
Note that, in principle, the update of suction and water transport could be done earlier in the algorithm, e.g.before updating particles stress (i.e. between points 10 and 11).However, as we wanted to preserve the algorithm of MPM, that update is moved later.The benefit of this approach is that the MPM algorithm can be taken as an externally called 'black box' procedure, while as explicit MPM has typically very small time steps, the introduced inaccuracy is negligible.
Finally, the selected constitutive model should allow for easy transmission between the flow state and quasi-static state.It also should allow for computations of stress with the initial stress being zero.Here, we use the same constitutive model for the moderate flow and quasi-static state, but two separate models coupled with a switching condition between the quasi-static and moderate flow states could be introduced instead.
In summary, the proposed algorithm should be sufficiently accurate for simulations of unsaturated granular flows in which: -The granular material is not close to full saturation and neglecting the transition to fully saturated material with positive pore pressures does not introduce much error -The accuracy of predicted granular material saturation after deposition is not of great concern -The initial state of the material is relatively uniform so that the mixing of the material during flow will not cause much water transport, so it can be neglected in the simulation -The simulation mainly concerns the flow of the granular unsaturated material so neglecting the water transport in the quasi-static part of the simulation does not lead to significant inaccuracies Furthermore, the GMPM algorithm can be used with an advanced contact formulation.The accuracy of the contact formulation between the flowing mass of granular material and the other surfaces (e.g.silo) will affect results significantly.In some cases, the flowing material may also erode the existing surface, which can have much influence on the final results.
Additionally, as the framework does not allow for full saturation, material point density is bound both by the critical density ρ crit and the density corresponding to the full saturation ρ sat , when the degree of saturation S r reaches 1.

Suction update algorithm-step U1
Suction is the main extra variable to be additionally computed in computations for unsaturated soils.The computations require a water retention curve, linking the degree of saturation of the material and the suction value.The degree of saturation alters when the volume of material changes, without the transport of water in or out.As in the proposed framework, the amount of water in each material point is constant, any volume change leads to the corresponding degree of saturation change.
Material point of volume V has the degree of saturation: where V w and V s are the volumes of water and solid grains, respectively.The volume of the material point is routinely updated in the material point method algorithm based on the deformation gradient, hence available for the calculations.Having the degree of saturation is usually enough to compute suction based on the water retention curve.Some advanced water retention models which differentiate between wetting and drying branches may require also the values of the degree of saturation in previous time steps.That would require storage of historical values of degree of saturation in a material point.

Theory
For the quasi-static and moderate flow stages, the simulations use a Mohr-Coulomb constitutive model extended for unsaturated soils, with a van Genuchten water retention curve.The model uses linear elastic law and a dilation angle for a nonassociated flow rule, necessary to prevent excessive plastic deformation during shearing.It additionally includes the shear rate dependency, significantly affecting results for any flows with a significant shear component, as shown e.g. in Ref. 21.
Hence, the yield surface equation is: where ϕ is the friction angle and the modified cohesion c * is a sum of the cohesion of the material, the cohesion due to shear rate γ and the cohesion due to suction s: where η is a parameter quantifying the increase in cohesion due to shear rate, ϕ b is the friction angle relative to matrix suction 22 and c is cohesion.The calculations use the simplified van Genuchten model 23 for the water retention curve equation: where n, m and α are parameters.
The selected constitutive model is not well suitable for material close to the critical density at which grains are losing some contact with the neighbouring grains.As such, the cohesion c* from the initial value c initial computed for each point in each timestep from Eq. ( 3), in some simulations is additionally scaled linearly to c critical close to the critical density, that is for particles with specific volumes above ν start as: , where c critical is the value of cohesion at the critical density (in the simulations it was set to be 20% of the c* = c initial computed for each material point in each timestep from Eq. ( 3)), ν critical is the specific volume of the particle at critical density.In the simulations, the parameter ν start at which scaling begins was the specific volume corresponding to 98% of the ν critical .This leads to a linear reduction of cohesion in material points close to the critical density.Taking into account the parameters used in the calculations, particles that have their specific volume above 98% of the critical specific volume ν critical have their cohesion adjusted as: The proposed algorithm for GMPM for unsaturated material requires a definition of a water retention curve.For the illustration, we used a water retention curve based on data from Ref. 24.As we have used the simplified van Genuchten model (Eq.( 4)) which assumes full saturation at zero suction, the water retention data used was adjusted for full saturation (Fig. 1).Then, the simplified van Genuchten model parameters (see Table 1) were obtained using the least square minimisation procedure.Despite the adjustment, the part of the curve used, mostly below 3.5 kPa suction, is replicated reasonably correct and serves sufficiently well for the qualitative comparisons of the algorithms presented in the paper (Fig. 1).
The selected constitutive relationships are rather basic.However, they are sufficient to demonstrate the capabilities of the method.They are also robust, numerically simple to implement   and do not require many adjustments during the transition from quasi-static to flow state.More advanced constitutive models, such as Barcelona Basic Model 25 would require further adjustments in the framework.In particular, the elastic law relating the volume change and mean net stress may cause issues at volumes corresponding to the critical density and possibly inaccurate predictions of volume change.Furthermore, in the stress-free state during moderate flow, the material point interacts with other material points, which, in principle, may not be stress-free.If that is the case, the point develops some strains.Having the strains, we need to compute stress, which is not possible from the state of zero mean stress assumed.One possible solution would be to link the mean stress in the particle to the stress corresponding to the critical volume.That, however, would be not compatible with the zerostress assumption during flow and may cause problems in the simulation and affect the convergence of the results.Therefore, using the framework with constitutive models assuming linear elasticity in semi-logarithmic or log-log space requires further research.

Simulation of a silo discharge of unsaturated granular material
This section gives details of the illustrative simulations of a silo discharge, comparing the proposed method with the original GIMP and Granular GIMP.For simplicity, the simulation assumes that the problem is symmetric, plane strain and that the silo has a flat bottom.The silo contains 80 cm x 46 cm of unsaturated granular material parametrised as Toyoura sand, see Fig. 2 and Table 1.The material flows through the silo opening (14 cm wide) and falls on the rigid ground from 1.25 m height.To minimise the time of the calculations, the deposition space under the silo has the same width as the silo, which is 50 cm.The selected dimensions of the silo are small, allowing more easily to show the influence of unsaturation in the granular material.
The problem is discretised with a square grid of 10 mm, with one material point per grid element.The side of the silo is fixed, with the symmetry boundary condition on the right-hand side.There is a frictional contact [26][27][28] between the wall of the silo and the granular material, with 0.2 taken as surface friction.Granular material is approximated with the amended Mohr-Coulomb model (Eq.( 2)), with the parameters based on Toyoura sand, as selected by Seyedan and Solowski, 29 based on data from Ref. 30 with an initial void ratio of 0.917.The initial suction is 3.5 kPa, with the simplified van Genuchten model water retention curve parameters given in Table 1.The simulations use values of ϕ b = 12[deg], a value suitable for granular soils. 22ere is no special differentiation between the quasi-static and moderate flow.However, when the density of a material point falls below the critical density, the GMPM framework uses the stress-free approximation of material behaviour. 5The critical density, at which the material changes its behaviour to the stressfree approximation is taken as 1300 kg/m3.The silo is made from a linear elastic material.Table 1 gives the parameters of the constitutive model used in the simulations.The simulations are done using Uintah code amended with the GMPM algorithm.

Simulation of unsaturated granular avalanche
This example models the flow of unsaturated granular material and its impact on elastic barriers.The example demonstrates the need for the GMPM algorithm in a simulation of an unsaturated granular avalanche and compares the granular flow of partially saturated and dry material.
The problem is set as a 2D plane strain, see Fig. 3.The slopes and grid in the simulations are, however, horizontal, with the slope achieved by adjusting the gravity vector to capture the 45 • inclination.The simulation uses a square grid of 25 mm, one material point per grid element, and fixed boundary conditions on the sides.The interaction between the barriers and the flowing granular material is governed by a contact algorithm [26][27][28] with 0.2 taken as surface friction.The material parameters for the sand are the same as in the previous simulation (see Table 1 and Section 3.1).

Results of simulations of a silo discharge of unsaturated granular material
To investigate the influence of the unsaturation on the results, we run 4 simulations with different algorithms: dry Granular GIMP (i), and three unsaturated GIMP algorithms: original GIMP, as well as Granular GIMP with no change of cohesion close to the critical volume (iii) and Granular GIMP with a drop of cohesion close to the critical volume (iv).
Upon examining the results (Figs. 4 and 5) it is clear that the deposition of the material is not well captured by the original GIMP algorithm, as the initial 80 cm of material is deposited unevenly and with incorrect material density, unlike in the simulations using the granular algorithm.Also, due to the extra numerical cohesion demonstrated in Ref. 2, the chunks of material falling are larger than in any other simulation.In the case of unsaturated Granular GIMP simulations, the discharge takes approximately 30% longer when compared to the unsaturated GIMP simulation-close to 5 s, versus 3.75 s.When the cohesion is reduced close to the critical density, discharge is slightly quicker, and leads to a more even deposition surface, as the grains on top are deposited initially close to the critical density, with lower cohesion.The deposited material is, generally, almost as dense as in the case of granular GIMP for fully dry material.In that case, the emptying of the silo is much quicker, with discharge finishing only a little later than 2.25 s.The differences in discharge rate are therefore very significant and show the importance of taking into account the unsaturation of the material.Also, the use of granular MPM greatly increases the quality of the results and allow for much more correct deposition of the material, when compared with previous results 31 .

Results of simulations of unsaturated granular avalanche
As in the previous case, the 4 simulations differ by MPM algorithm and their treatment of cohesion close to critical volume: Granular GIMP for a dry material (i), unsaturated GIMP with no change of cohesion close to the critical volume (ii) unsaturated Granular GIMP with no change of cohesion close to the critical volume (iii) and unsaturated Granular GIMP with a drop of cohesion close to the critical volume (iv) are quite easy to discern.As in the previous case study, the original GIMP simulation leads to flow with the largest visible chunks of granular material, while the reduction in the overall cohesion leads to faster and more granular flow.Despite those differences, up to 1 s of time (see Fig. 6), results are relatively similar.It is though quite clear that cohesion reduction at close to the critical density leads to a similar outcome as for the fully dry material.That is because the material points after moving over the first obstacle usually have their density close to the critical density, with cohesion being the main factor preventing grain separation during flow.The GIMP simulation and the unsaturated Granular GIMP simulation without cohesion reduction are also looking similar at 1 s time, with points depositing at the second obstacle.
Ultimately (see Fig. 7), the flow in the unsaturated GIMP simulation finishes earliest, with the material stopped by the second obstacle at 1.66 s.The material points at that second obstacle have unnaturally large domains, and unphysically high volume, showing that GIMP is not well suitable for the simulation.Unsaturated Granular GIMP, with no cohesion reduction, has similar ultimate outcomes, though material points settle for a longer time and have their density correctly recovered.The Granular GIMP simulation for dry material leads to most deposition of material at the third obstacle.The outcome of the simulation run with the unsaturated Granular GIMP with cohesion reduction algorithm is similar to that of dry Granular GIMP, though the amount of material reaching the third obstacle is smaller, as expected for unsaturated material.

Discussion of simulation results
The simulations confirm that Granular Material Point Method improves the simulations in which the flow of the granular material occurs when compared to the original formulation of the method.In shown simulations, the main contributing factor is that in the Granular simulations, the volume of the given material point cannot raise beyond the defined critical volume.Additionally, the Granular MPM algorithm prevents the interaction of the material points that have their domains of interaction not overlapping.In the case of flow, this leads to material points quite often acting similarly to single grains of granular material, moving without any neighbours to interact.
The simulations, similar to those given in Ref. 2 were completed without major numerical problems, be it stress oscillations or convergence issues, even though we use no numerical damping.We believe that the lack of numerical damping (understood as the dissipation of part of velocity or energy of each or some material points in each timestep, typically used to smooth out any oscillations) is key in simulations involving gravitational flow and is strictly necessary to obtain physically correct results that are independent of the grid density and timestep.This good performance is partially due to the advanced formulations of the Material Point Method used in the simulations designed to improve accuracy.Our high-quality own code also contributes, as well as the knowledge gained from investigations related to the stability and accuracy of MPM (e.g.Refs.32, 33).We also need to mention the very high-quality Uintah code we build upon (e.g.Refs.34, 35).To further improve the stability of the calculations, the simulations used the Mohr-Coulomb model implemented based on analytical stress integration algorithm 36,37 coupled with a linear elastic law.The small-strain extension 38 was not used as it increases complexity and requires several extra parameters without affecting the results qualitatively.

Conclusions
The paper presented a simple framework for simulations of unsaturated granular flows.The framework extends the Granular Material Point Method framework so it can be used for partially saturated soils.The paper proposes a set of simplifying assumptions for the calculations.Those limit the applicability of the show algorithm to cases in which water transport is not dominant.However, it greatly simplifies the calculations, meaning that a relatively standard MPM code can be used.The proposed framework is used to simulate the flow of granular material out of a silo and its deposition, as well as a flow of an avalanche over obstacles.The outcomes of simulations show significant qualitative differences between the algorithms, with Granular GIMP algorithms being superior to the original GIMP.The shown two sets of assumptions related to the behaviour of the material during flow, one with constant cohesion, and the other with reduced cohesion, also lead to different outcomes.It is clear that considering unsaturation in the GMPM framework leads to visibly different results and is necessary for calculations involving partially saturated material flow.
In the future, the framework needs quantitative validation.That may lead to revisiting the assumptions, as well as adjustments related to the introduction of more accurate constitutive models for unsaturated material.The authors aim, in particular, to give more meaning to the critical density that currently is just a user-defined parameter.The framework should also be tested and adjusted for constitutive models for granular materials based on the critical state soil mechanics.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Van Genuchten water retention curve used in the study with indicated initial suction equal to 3.5 kPa.During the flow, the suction value tends to be lower than that, as the material density drops.

Fig. 4 .
Fig. 4. Comparison between Granular GIMP (top) for dry granular material and GIMP (bottom) for unsaturated material.The simulation results are shown at different times.The particle colour indicates its relative velocity (blue-low, red-high) (movies of simulations are available in Appendix A).

Fig. 5 .
Fig. 5. Comparison of simulation results at time 0,5 s-5 s between unsaturated Granular GIMP without cohesion change (top) and with cohesion decrease close to critical density (bottom).The particle colour indicates its relative velocity (blue-low, red-high).

Fig. 6 .
Fig. 6.Comparison of simulation results at time 1 s between Granular GIMP (top left), unsaturated GIMP (top right), unsaturated Granular GIMP without cohesion change (bottom left) and unsaturated Granular GIMP with cohesion decrease close to critical density (bottom right).The particle colour indicates its relative velocity (blue-low, red-high).

Fig. 7 .
Fig. 7. Comparison of simulation results at the end of simulation between Granular GIMP (top left), unsaturated GIMP (top right, t = 1.66 s), unsaturated Granular GIMP without cohesion change (bottom left) and unsaturated Granular GIMP with cohesion decrease close to critical density (bottom right).The particle colour indicates its relative velocity (blue-low, red-high).

Table 1
Constitutive model parameters used in the silo discharge simulations, approximating Toyoura sand.