An explicit stabilised material point method for coupled hydromechanical problems in two-phase porous media

This paper presents a single-point Material Point Method (MPM) for large deformation problems in two-phase porous media such as soils. Many MPM formulations are known to produce numerical oscillations and inaccuracies in the simulated results, largely due to numerical integration and stress recovery performed at non-ideal locations, cell crossing errors, and mass moving from one background grid cell to another. The same drawbacks lead to even worse consequences in the presence of an interstitial fluid phase, especially when undrained/ incompressible conditions are approached. In this study, an explicit stabilised MPM, based on the Generalised Interpolation Material Point (GIMP) method with Selective Reduced Integration (SRI), is proposed to mitigate typical numerical oscillations in (nearly) incompressible coupled problems. It includes two additional features to improve stress and pore pressure recovery, namely (i) patch recovery of pore pressure increments based on a Moving Least Squares Approximation, and (ii) two-phase extension of the Composite Material Point Method for effective stress recovery. The combination of components leads to a new method named GC-SRI-patch. After a detailed description of the approach, its effectiveness is verified through analysing various consolidation problems, with emphasis on the representation of pore pressures in time and space.


Introduction
Large deformation problems in two-phase porous media are of great importance in geo-engineering, for instance in the analysis of landslides or foundation installation processes. Traditional numerical methods such as the Finite Element Method (FEM), however, are often negatively impacted by the development of large deformations, which may cause numerical simulations to abort abruptly or provide misleading results. Several numerical methods, such as the Smoothed Particle Hydrodynamics (SPH) method (Gingold and Monaghan, 1977;Monaghan, 1994), the Material Point Method (MPM) (Sulsky et al., 1994;Sulsky et al., 1995), the element-free Galerkin method (Belytschko et al., 1995), the Particle Finite Element Method (PFEM) (Oñate et al., 2004;Monforte et al., 2017;Della Vecchia et al., 2019), and meshfree methods (Navas et al., 2016;Navas et al., 2018) have therefore been proposed to solve large deformation problems in porous media, with each one featuring a specific mix of advantages and drawbacks. Recently, MPM has been gaining recognition as an appropriate approach for this class of problem .
The present work focuses on soil-like two-phase media, in which the pores of a stiff solid skeleton are filled with a (nearly) incompressible fluidmost often water. In common with FEM poromechanical formulations (Zienkiewicz et al., 1999), the MPM solution of coupled problems usually builds on the so-called u − p and v − w formulations, named after their primary unknown variables. In the u − p formulation the relative acceleration between the solid skeleton and the pore water is neglected, so that the only unknown degrees-of-freedom (DOFs) are those associated with soil displacement (u) and pore pressure (p) (Abe et al., 2014;Higo et al., 2010Higo et al., , 2015Zabala and Alonso, 2011;Zhang et al., 2009;Zhao and Choo, 2020;Zheng et al., 2013). However, to address dynamic applications, recent MPM studies have been performed based on the v − w formulation (with the primary variables being the total velocity of the soil skeleton, v, and the discharge velocity of the pore water relative to the skeleton, w), which includes soil-water relative acceleration terms (Zhang et al., 2007;Jassim et al., 2013;Bandara and Soga, 2015;Soga et al., 2015;Yerro, 2015;Yerro et al., 2015;Yerro et al., 2017;Ceccato et al., 2016;Liu et al., 2017;González Acosta et al., 2019). In the framework of coupled MPM, both single-layer and twolayer approaches have been explored, i.e., the use of one or two sets of material points (MPs) to describe the response of distinct phasessee Soga et al. (2015) for more details. The high computational cost of twolayer MPM has so far determined the broader popularity of the singlelayer approach, which is also followed in this study. The core of MPM relates to the use of a background mesh to solve (the discrete version of) the relevant governing equations, while MPs serve as quadrature points and can freely move within the domain set by the mesh, see e.g., Sulsky et al. (1994Sulsky et al. ( , 1995. Low-order shape functions are often preferred in MPM, so as to avoid the numerical divergence possibly caused by the negative parts of higher-order polynomial shape functions. These play an especially important role in MPM, as quadrature takes place directly using the material point locations and often grid cells are only partially filled, so material points far away from the node (which usually take the negative part of the shape function) can strongly influence the quadrature. When applied to coupled hydromechanical problems, however, low-order MPM may suffer from numerical instabilities, typically in the vicinity of the so-called undrained-incompressible limit. Unstable/inaccurate results will normally be obtained under such conditions, due to MPM's low-order interpolation functions violating well-known inf-sup requirements, in a way similar to that widely observed for finite element calculations (Brezzi and Bathe, 1990;Bathe, 2001). Alternatively, the use of higher-order interpolation functions, such as B-spline functions, which do not have the negative parts that polynomial interpolation functions have, has also been considered (Steffen et al., 2008;Tielen et al., 2017;Gan et al., 2018;Tran et al., 2019). However, resorting to higher-order interpolation directly leads to a larger number of material points (MPs) and 2D or 3D implementations have not been developed. For this reason, the present work focuses on the improvement of low-order, coupled MPM.
Most existing literature on coupled MPM formulations concerns twophase materials with compressible components. By considering water to be more compressible, spurious oscillations can be reduced. In addition, unsaturated conditions (such as in Wang et al. (2018)) implicitly are much more compressible, and therefore result in less oscillations. Nonetheless, a few studies dealing with (nearly) incompressible problems and inf-sup-related instabilities have already been published. For example, a fractional time stepping method combined within an enhanced volumetric strain formulation was proposed in Jassim et al. (2013) to mitigate pathological locking and spurious oscillations in the pore pressure field. However, this time stepping method is not equally effective for all possible drainage conditions, nor straightforward to implement into existing coupled MPM codes. A stabilised implicit MPM has been recently developed by Zhao and Choo (2020), in which the mass balance equation is augmented with a stabilising term to make equal-order mixed interpolation stable under undrained conditions. Such a term is derived using polynomial pressure projection in a way specific to the adopted time integration algorithm.
In standard MPM, integration and recovery of strains, stresses and pore pressures always occurs at the MPs. Reducing the number of integration points and fixing the location within individual grid cells (i.e. at so-called Gauss points as used in finite elements) can be used for stabilisation purposes. The idea of benefiting from reduced integration in coupled-MPM has been previously introduced, for instance, by Abe et al. (2014), Bandara and Soga (2015), and Wang et al. (2018). Accordingly, pore pressures are evaluated at the central Gauss points (GPs) of the background mesh (instead of at MPs), which has been found to alleviate the aforementioned pore pressure instabilities. Additionally, as reduced integration is exclusively performed to evaluate pore pressure variations, computed results appear not to suffer from spurious hourglass modes . This approach can be readily implemented into existing explicit MPM codes and is further pursued in the present study. I1]t emerges from previous literature that pore pressures are most usually recovered from GPs to MPs by assuming uniform pore pressure increments within each grid cell of the background mesh. This determines the direct influence of grid cell size on MPM solutions, and can sometimes lead to pore pressure discontinuities (at inter-cell bound-aries) and inaccuracies (at MPs). When MPs cross grid cell edges, it can cause a sudden change in pore pressure at the MPs, and, as a consequence, spurious variations of nodal internal forces and stress oscillations, especially when a coarse background mesh is adopted. In this respect, some authors tested reduced integration in the Generalised Interpolation Material Point (GIMP) method (Bardenhagen and Kober, 2004), as a way to alleviate the stress oscillations related to cell-crossing (Abe et al., 2014). More recently, GIMP has also been introduced in the implicit MPM to mitigate cell-crossing inaccuracies in two-phase problems built on the simplified u − p formulationsee the aforementioned work (Zhao and Choo, 2020).
This study proposes an explicit stabilised MPM, named Generalised Interpolation Material Point method with Selective Reduced Integration, based on the patch recovery of pore pressure increments (GC-SRI-patch) and a complete dynamic formulation of the v − w type. The GIMP method (Bardenhagen and Kober, 2004) contributes to reducing (stress) oscillations promoted by grid crossing errors. To avert spurious pore pressures, selective reduced integration (SRI) has also been introduced for pore pressure recovery at central GPs. Patch recovery based on a Moving Least Squares Approximation (MLSA) (similar to the moving least squares technique used in Tran et al. (2020) for improved moving least square shape function construction) is proposed to map calculated pore pressure increments from central GPs to MPs in order to increase the accuracy of the results. As for the evaluation of effective stresses, the same approach adopted in the recent (one-phase) Composite Material Point Method (CMPM) (González Acosta et al., 2017) is herein extended to the proposed GC-SRI-patch to enhance stress recovery. This work is limited to elastic constitutive behaviour of the solid skeleton, and focuses on exploring the effectiveness of (the ingredients combined in) the proposed method. Further investigation is necessary to guarantee stable/accurate solutions in the presence of material plasticity (Coombs et al., 2018), as well as to explain fundamentally why the GC-SRI-patch approach is beneficial against inf-sup related instabilities.
The content of the paper is organised as follows. After providing the governing hydromechanical equations (Section 2), the technical details regarding the numerical formulation and implementation of the proposed GC-SRI-patch method are described in Section 3. Section 4 presents 1D/2D numerical examples for the verification of the proposed GC-SRI-patch method.

Coupled formulation for two-phase porous media
In line with most MPM literature, this work builds on a velocity formulation of the governing hydromechanical equations (Jassim et al., 2013;Wang et al., 2018;González Acosta et al., 2019). In particular, a fully dynamic formulation is used, in which the total velocities of the soil skeleton and pore water, v s and v w respectively, are used as primary variables. However, it should be noted that the relative discharge (Darcy) velocity w may be used in lieu of v w , as explained, e.g., by Zienkiewicz et al. (1999). Soil-like fully saturated porous media are considered in what follows. The density of the soil-water mixture ρ is obtained from the individual phase densities as ρ = nρ w + (1 − n)ρ s , where the subscripts s and w denote solid and water phases, respectively, and n is the volume porosity. For more convenient translation into their discrete versions, all governing equations are presented hereafter using a vector notation. In particular, the total stress acting on the mixture includes effective stress (σ ′ ) and pore pressure (p w ) components, σ = σ ′ + mp w , where m is the (2D) vector version of the Kronecker tensor.
Both stresses and pore pressures are assumed to be tension positive.

Governing equations
Momentum balance for the whole two-phase mixture is fulfilled if where S is the differential operator defined for 2D problems (Zienkiewicz et al., 1999): while v s and v w are the total (true) velocities of the soil and pore water, respectively, b is a body acceleration term (e.g., gravity acceleration), and dots above symbols indicate time differentiation. It is also worth noting that the relative discharge velocity w may be obtained from the individual true velocities as w = n(v w − v s ).
A similar equilibrium equation can be formulated for the water phase only: where R represents the drag force exchanged by the solid and fluid phases during water seepage, which is proportional to the relative discharge velocity w according to Darcy's law: where the hydraulic conductivity k is here assumed to be isotropic for simplicity, and g is the gravitational acceleration. For mass balance, the following equations ensure the conservation of the solid and water massesunder the assumptions of uniform densities and porosity, and incompressible soil grains: Density variations in a barotropic fluid obey the following relationship: where K w is the bulk modulus of the fluid. Substituting Eq. (5) into Eq. (6) and combining with Eq. (7) allows the pore pressure rate to be obtained: The constitutive relationship between strain (ε) and effective stress (σ ′ ) rates can be expressed as (Zienkiewicz et al., 1999) where D is the tangent stiffness matrix of the solid skeleton, and ε 0 is a strain rate term related to, e.g., thermal effects. Since emphasis is hereafter on the development and verification of the proposed GC-SRIpatch method, (i) isotropic linear elastic behaviour of the solid phase and (ii) linearised/infinitesimal definition of strain rates are exclusively considered. Fully general modelling of large deformations can be achieved by adopting well-established finite strain measures (Holzapfel, 2000), though with no expected detriment to the hydromechanical performance of the proposed method.

Boundary and initial conditions
Considering a fully saturated porous domain Ω, its boundary surface Γ can be decomposed into subsurfaces on which Dirichlet and/or Neumann boundary conditions are imposed. Surface decomposition should be such that Γ = Γ u ∪ Γ τ = Γ p ∪ Γ w and Γ u ∩ Γ τ = Γ p ∩ Γ w = ∅, so as to enable the enforcement of relevant conditions on the solid and water velocities: as well as on the (total) surface traction and water pressure: In Eqs. (12) and (13), G is a matrix containing components of the unit vector normal to the boundary surface Γ (Zienkiewicz et al., 1999), while ṽ s (t),ṽ w (t),τ(t) and p w (t) are prescribed boundary values of the solid and water velocities, surface traction and pore pressure, respectively. The full set of initial conditions are in which α equals s or w to indicate either the solid or water phase.

Integral weak formulation
Standard manipulation of the governing equations and boundary conditions allows the following integral/weak version of the momentum balance equations, Eqs. (1) and (3), to be obtained (Zienkiewicz et al., 1999): where δv is a vector of suitable test functions.

The GC-SRI-patch method: formulation and implementation
This section provides technical details regarding the formulation and implementation of the proposed GC-SRI-patch method. Emphasis is on the combined application of lessons learned from previous studies, with a view to mitigating the deficiencies of standard MPM with respect to cell crossing errors, pore pressure instabilities, and numerical quadrature and stress/pore pressure recovery performed at non-ideal locations.

Spatial discretisation
Primary unknowns in the adopted velocity formulation (Jassim et al., 2013;González Acosta et al., 2019) are the velocities of the solid (v s ) and water (v w ) phases. The same shape functions are used to approximate the velocities of both phases, as well as the test function vector δv: where v α and δv define vectors of nodal values. In regular MPM, N(x) contains linear shape functions of the same kind as those used in FEM. It is known that regular MPM may suffer from stress oscillations when MPs cross grid cell boundaries due to discontinuous shape function gradients. GIMP was thus proposed ( Bardenhagen and Kober, 2004) to alleviate such oscillations, with shape functions constructed by integrating linear FEM shape functions N i (x) over the MP support domain Ω p (as shown in Fig. 1). In one dimension, the GIMP shape functions S i,mp and their gradients ∇S i,mp are calculated as over the problem domain Ω, where V mp is the MP volume, and χ mp is the "particle characteristic function": The support domain Ω mp is assumed to be of size 2l p in each dimension, and can be computed by dividing the grid cell size by the initial number of MPs along the considered direction. Fig. 1 compares a shape function and its gradient for both GIMP and MPM for a 1D case. The support domain for a specific MP is also shown in Fig. 1. In 2D/3D problems, shape functions are obtained by multiplying the individual 1D functions for different directions.
After substituting the approximation Eq. (18) using GIMP shape functions S(x) into Eqs. (16) and (17), the discrete versions of Eqs. (1) and (3) are: where M s , M w , and M w are global nodal mass matrices which are diagonalised through "mass lumping" -see, e.g., al-Kafaji (2013). In addition to s and w, the subscript m is used to denote "soil-water mixture". The superscripts trac, body, int and drag denote global force terms associated with surface tractions, body forces, internal forces and soil-water drag, respectively. The global mass matrices in Eqs. (22) and (23) are obtained by assembling the following local mass matrices associated with individual grid cells: where subscript i defines the i th grid cell node, and N mp is the number of MPs in the grid cell whose spatial coordinates, mass, volume, density and porosity are denoted by x mp , m α,mp , V mp , ρ α,mp and n mp , respectively (with α = s, w). The remaining force terms are obtained as follows: where N bmp is the number of MPs near the domain boundary Γ on which the traction forces are applied. In Eqs. (27)-(33), p w,bmp and τ m,bmp are the prescribed pore pressure and traction force at MPs near the domain boundary, while p w,mp ,σ mp ,ṽ w,mp and ṽ s,mp are the respective pore

Time discretisation
An explicit algorithm is adopted to integrate Eqs. (22) and ((23) owing to the diagonality of the above lumped mass matrices. The nodal accelerations are then used to update phase velocities over the time increment Δt at the MP locations: where N n is the total number of nodes in the problem domain, and subscripts mp and i denote MPs and background mesh nodes, respectively. It should be noted that the same shape functions S i ( x mp ) are used to map kinematic information from MPs to nodes in the background mesh and vice versa. After mapping the velocities at MPs v t+Δt α,mp to background nodes, the new nodal velocities v t+Δt and the strain rates at MPs are evaluated as is the shape function gradient. Finally, the strains, (effective) stresses, pore pressures and porosity values are updated as where J is the Jacobian of the deformation gradient tensor, i.e. J (

Mitigating pore pressure instabilities in MPM
Due to its apparent similarity to FEM, MPM suffers from inf-suprelated instabilities when low-order (linear) interpolation is adopted. This can also be the case for hydromechanical incompressible problems in porous media, giving rise to undesired oscillations in the pore pressure field (Bathe, 2001;Belytschko et al., 2013;Chen et al., 2018).
Currently, a number of techniques, such as multi-field variational principles, fractional step time integration, high-order interpolation and selective reduced integration are employed in FEM (Li et al., 2003;Zienkiewicz et al., 2005;Bathe, 2006;White and Borja, 2008;Pisanò and Pastor, 2011;Belytschko et al., 2013;McGann et al., 2015) in order to mitigate this type of pore pressure instability. These techniques may also be applied to MPM. Herein, the performance of the proposed method is improved by combining GIMP, which partially enhances the order of the interpolation, with Selective Reduced Integration (GIMP-SRI) of the pore pressures. In its original conception (Bardenhagen and Kober, 2004), GIMP was used to integrate stresses in one-phase media, resulting in a significant improvement due to reduced cell-crossing errors. Nevertheless, in the two-phase case, large pore pressure oscillations remain inside grid cells even using GIMP. Further benefits can be achieved by adopting reduced integration (GIMP-SRI) for pore pressure recovery, as illustrated in Fig. 2. However, instead of directly calculating incremental pore pressures at MPs, reduced integration requires pore pressure increments to be computed at central integration GPs in each grid cell (e.g. GP1, GP2, GP3, and GP4 in Fig. 2) as where ε t s,gp and ε t w,gp are volumetric strain rates of the solid/water phases, and n t gp is a weighted porosity. They are evaluated at the central GP position x gp by and where V t mp,e is the intersection volume between the MP support domain and the current grid cell where the considered GP is located. As reduced integration is only performed to evaluate pore pressure variations, computed results are found not to suffer from spurious hourglass modes .
After obtaining incremental pore pressures through Eq. (43) at the central GPs, the key issue is to recover them back to the MPs. It is wellknown that stresses calculated at the centre of low-order rectangular elements in FEM are of high accuracy and convergence order. As the calculation phase in MPM is a FEM calculation, with modified integration, it can be concluded that this also holds for an MPM grid cell. Following Zienkiewicz and Zhu (1992a), it is here proposed to calculate the pore pressure increments at the MPs by so-called patch recovery based on a moving least squares approximation (MLSA). As shown in Fig. 3, a patch of four quadrilateral cells can be identified for each internal node i. Within such a patch, a rectangular area can be delimited around the node by using the central GPs in the four grid cells. It is thus possible to introduce for the pore pressure increments (Δp w ) the following polynomial approximation of order p in the considered rectangular domain Ω i (bounded by the red dashed lines in Fig. 3): where (x, y) is the location of GPs in Ω i , while Q and a are vectors containing polynomial basis functions and interpolation degrees-offreedom. In general, different shape functions may be chosen to approximate the incremental pore pressure field. In this study, a linear version of Q ( is chosen, giving rise to the interpolation plane in Fig. 3 after the determination of the coefficients in a = [a 0 a 1 a 2 ] T . Based on a posteriori error estimator, the relative error at sampling GPs is calculated as where N gp is the total number of GPs in the approximation domain Ω i , and (x i , y i ) are the coordinates of the GPs. Minimising the error with respect to a leads to the following linear system: Finally, the pore pressure increments at the MPs located in the approximation domain Ω i can be obtained as and these can be used to derive the pore pressure at time t + Δt. For MPs near the domain boundary there are insufficient grid cells to form a complete patch. For these cases, the pore pressure increments are determined by extending internal patches up to the MP position. Similar concepts for determining stresses at the boundary nodes in FEM can be found in previous studies (Zienkiewicz and Zhu, 1992a;Zienkiewicz and Zhu, 1992b;Zienkiewicz et al., 2005).

CMPM stress integration/recovery
In general, MPM also suffers from oscillations and inaccuracies due to performing numerical integration and stress recovery at non-ideal (MP) locations. The recently developed Composite Material Point Method (CMPM) (González Acosta et al., 2017Acosta et al., , 2020, which extends the solution domain for each grid cell through considering the influence of neighbouring cells, can significantly alleviate stress oscillations and help to recover stresses for one-phase problems.
In CMPM, new grid cell shape functions are established based on an  extended influence domain (i.e., a patch) using Lagrangian interpolation. All nodal displacements within this extended domain are used for stress recovery, which can lead to improved stress values at MPs. The constructed shape functions are summarised in Appendix A, while more details can be found in González Acosta et al. (2017Acosta et al. ( , 2020. Here, CMPM is firstly extended to the case of coupled two-phase problems and then exploited to improve the recovery of effective stresses at MPs in selected verification examples.

Numerical implementation
Each step in the proposed GC-SRI-patch method is explicitly solved according to the following sequence of sub-steps (see also the flow chart in Fig. 4): (1) initialise all variables at the nodes of the background mesh (Eqs. (24)-(33)); (2) calculate the water nodal accelerations v t w using the discrete equilibrium equation for the water phase (Eq. (34)); (3) substitute the water nodal accelerations v t w into the discrete equilibrium equation for the soil-water mixture and calculate the soil nodal accelerations v t s (Eq. (35)); (4) update both soil and water nodal velocities using explicit forward Euler integration; (5) update the velocity and position of all MPs; (6) update the nodal velocities for both the soil and water by mapping variables back from the MPs; (7) calculate the effective stresses at the MPs by using GIMP combined with CMPM (Eq. (40)); (8) calculate the pore pressure at the GPs via Eq. (43), and then recover the pore pressures at the MPs from the GPs using Eqs.

Verification examples
This section presents three (plane strain) verification examples confirming the suitability of the proposed GC-SRI-patch method. In all examples the considered porous medium (soil) is fully saturated, with a solid skeleton modelled as isotropic linear elastic. Square background meshes are used in all cases, with each grid cell initially hosting four, equally-spaced material points. Given the emphasis of this work on the development of the GC-SRI-patch method, only relatively simple boundary conditions are considered in these analyses; further work will be devoted in the future to tackling more complex hydro-mechanical boundary conditions.

1D consolidation of a soil column
The first example is Terzaghi's 1D consolidation problem, which is commonly used to verify numerical methods for coupled poromechanical problems (Jeremić et al., 2008;Bandara and Soga, 2015). Fig. 5a illustrates the problem geometry and boundary conditions. The width and height of the problem domain are 0.1 m and 1.0 m, respectively. The pore water is allowed to drain through the top surface, whereas all other boundaries are impermeable. The displacement boundary conditions are a fixed mesh base and rollers at the two vertical boundaries allowing only vertical displacement.
A uniformly distributed static load p a , of either 1.0 kPa or 200.0 kPa (with no gravity loading), has been applied to the top surface to test GC-SRI-patch's performance for small and large deformations, respectively. The MPM discretisation is shown in Fig. 5b. The problem domain is discretised by ten 4-node quadrilaterial grid cells of size 0.1 m × 0.1 m, while a time-step size of Δt = 1.0 × 10 6 s has been used for time marching. It should be noted that, instead of applying external tractions on the top layer of MPs, they are directly applied on the movable top surface of the column , as illustrated in Fig. 6. For each MP, the support domain is defined by 2l px and 2l py in the horizontal and vertical directions, respectively. Thus, the location of the top boundary can be determined by the coordinates of the uppermost layer of MPs combined with the value of l py . At each time step, the support domain of the MPs at time t +Δt is updated as  X. Zheng et al. position of the top surface, the applied external distributed load τ m,ts is mapped from the top surface to the surrounding background nodes τ m,i using regular shape functions: where x ts is the position of top surface, and Δh is the grid size of the background square mesh. It should be pointed out that a linear shape function N i (x ts ) is used here for the mapping of the applied external load. Fig. 7 compares the GC-SRI-patch solution to the corresponding analytical solution for different values of the (dimensionless) time factor T v , defined as

Small deformation analysis
where H v is the drainage path length (here equal to the thickness of the soil layer), and c v is the coefficient of consolidation defined by 1 The small deformation analytical solution, reported by Terzaghi (1943), relies on the assumption that the layer thickness (H v ), hydraulic conductivity (k), and Young's modulus (E) of the soil layer remain constant during the consolidation process. For a clearer comparison, a dimensionless pore pressure p and a dimensionless current layer thickness H are introduced: The analytical solution for the considered initial/boundary conditions can be expressed in terms of pore pressure as a function of dimensionless depth and time (Fig. 7a): where M = (m − 1 2 )π. The average degree of consolidation, as shown in Fig. 7b, is defined as The results in Fig. 7 confirm the excellent agreement between the analytical and GC-SRI-patch solutions in terms of both excess pore pressure and average degree of consolidation. The proposed smooth distributions of pore pressure further demonstrate the advantage of the GC-SRI-patch even when using a relatively coarse background mesh. In comparison, the GIMP solution shows piecewise constant pore pressures over each cell. Similar results can also be found in Bandara and Soga (2015), where a much finer background mesh is needed to obtain a satisfactory solution (note that in Bandara and Soga (2015) a mesh ten times finer than in the present CG-SRI-patch case is used).

Large deformation analysis
The reference large deformation theory of soil consolidation was developed by Gibson et al. (1967). Among other aspects, the presence of large deformations makes it no longer appropriate to consider a constant H v and k, as their values may change significantly due to soil deformation and reduction in porosity (n). This aspect is captured in the analytical solution presented in Xie and Leo (2004), which builds on Gibson's theory and the assumption of porosity-dependent hydraulic conductivity k t (n), given as where k 0 and n 0 are the initial hydraulic conductivity and porosity, respectively, which are the same as used above for the small deformation case. According to Xie and Leo's large deformation solution (Xie and Leo, 2004), the dimensionless pore pressure varies in space and time as follows: where m vl = 1/E is the 1D compressibility 2 and p a is the applied external load.
The analytical and GC-SRI-patch solutions, which are both based on the relationship in Eq. (57), are compared in Fig. 8. The numerical and analytical excess pore pressure isochrones at different average degrees (b) Average degree of consolidation (a) Excess pore pressure isochrones  (Xie and Leo, 2004), m vl is assumed not to vary with soil porosity as a first approximation. of consolidation U s , from U s = 0 to U s = 0.9, are compared in Fig. 8a. In the presence of large deformations, U s is obtained as where S t is the top surface settlement at time t, calculated analytically as and S ∞ is the asymptotic value of S t as T v →∞: The comparison in terms of top surface settlement is given in Fig. 8b.
Overall, Fig. 8 shows that the GC-SRI-patch results compare well with the analytical large-deformation solution, notwithstanding the simplified representation of the strain field mentioned in Section 2.1. However, slight pore pressure oscillations are visible near the top boundary in Fig. 8a, which may be related to the external load being applied at the top of the MPM domain and then transferred to the background mesh nodes using shape functions.
It should be noted that, in the large deformation case, large pore pressure oscillations (at MPs) near the top domain boundary lead the explicit GIMP simulation to abort after significant displacement of the MPs. For comparison purposes, Fig. 9 shows pore pressure profiles corresponding to U s = 0.15 for the GC-SRI-patch, GIMP, and analytical solutions. The comparison further demonstrates the applicability of the explicit stabilised GC-SRI-patch method proposed in this study.

Pressurised hollow cylinder
In this benchmark, a two-phase hollow cylinder subjected to an internal pressure is studied. The problem geometry and boundary conditions are shown in Fig. 10. The inner and outer cylinder radii are r i = 0.20 m and r e = 1.20 m, respectively, giving a cylinder wall thickness of 1.0 m. The height of the cylinder H is equal to 1.0 m, while the problem domain is discretised using grid cells of dimensions Δr = Δy = 0.20 m. The boundary conditions are that the nodes at the top and bottom of the domain are only allowed to move in the radial direction, while the nodes at the outer boundary are fully fixed. The pore water is not allowed to flow in/out of the cylinder, so as to replicate (globally) undrained conditions. The soil and water properties are the same as in Section 4.1.
The benchmark has been solved by applying an internal total pressure p i = 100 kPa. Since drainage is not allowed, the pressure applied to the cylinder is transferred onto the (nearly incompressible) pore water. The near incompressibility also implies that the MPs do not displace significantly from their original positions. Fig. 11 compares the normalised radial excess pore pressure distribution through the cylinder wall obtained using MPM, GIMP and GC-SRI-patch. Because the locations of the MPs hardly change due to the pressurisation, the pore pressures obtained via MPM and GIMP are almost identical. The results obtained using GC-SRI-patch correctly show a constant pore pressure, equal to the applied pressure p i , through the cylinder wall. In contrast, large pore pressure oscillations are observed in both the MPM and GIMP solutions, with values near the pressurised boundary being significantly smaller than the applied pressure.

2D slumping block (self-weight consolidation)
In this section, the 2D large-deformation consolidation of an elastic slumping block (of width and height equal to 4 m and 2 m, respectively) under the sole action of gravity is studied (Zhao and Choo, 2020).  Taking advantage of symmetry, only the right half of the domain is considered as shown in Fig. 12a. Both the top and right boundaries are unconstrained and freely draining, while the left and bottom boundaries are impermeable and supported by rollers. No surface loads are applied, so that the consolidation process is exclusively driven by gravity loading applied with a ramp-like time history to avoid dynamic oscillations (Fig. 12b). The gravitational force gives rise to pore pressure build-up, the dissipation of which promotes gradual deformation of the block.
Time-domain simulations were performed with a time-step size equal to Δt = 1.0 × 10 6 s. For comparative purposes, standard MPM, GIMP, and GC-SRI-patch methods were used to analyse the problem. Fig. 13 shows the excess pore pressure distributions at t = 0.05 s obtained in the MPM, GIMP and GC-SRI-patch analyses. Due to only limited displacement experienced by the MPs up to that time, the MPM and GIMP analyses return very similar results. However, large pore pressure oscillations in a typical checkerboard pattern are visible in Fig. 13a. Because of the lack of MP grid crossing during the short time considered, the observed oscillation may be attributed to incompressibility and related instabilities. The extent of the oscillatory behaviour increases as the consolidation process evolves, and causes the explicit GIMP simulation to abort soon after the end of the loading ramp. This confirms the ineffectiveness of the low-order, non-stabilised GIMP scheme for incompressible problems (González Acosta et al., 2019).
The excess pore pressure response resulting from the proposed GC-SRI-patch is presented in Fig. 13b. In contrast with the GIMP checkerboard pattern shown in Fig. 13a, the GC-SRI-patch solution appears to be oscillation-free with compressive pore pressures everywhere. Due to the relatively quick gravity loading, water cannot rapidly drain and the numerical simulation develops under approximately undrained conditions. Therefore, the applied gravity loading is mostly translated into   pore water pressure increase, as can be observed in Fig. 13b. Unlike a uniform pore pressure in each grid cell, the pore pressures show a continuous distribution both within each grid cell and at inter-cell boundaries. The visible smoothness of the pore pressure field further confirms the suitability of the proposed GC-SRI-patch, even very near the undrained-incompressible limit.
Additionally, Fig. 14 shows the simulated excess pore pressure distributions at different times during the self-weight consolidation, from t = 0.1 s until t = 0.5 s. As observed in Ma et al. (2010) with regard to the GIMP method, if MPs are located at the far sides of a cell, the masses of some nodes (typically those close to the domain surface) may become small while shape function gradients and nodal forces may not. As nodal accelerations are explicitly obtained by dividing total nodal forces by nodal masses, large accelerations are obtained, which can in turn cause numerical instability. To alleviate acceleration inaccuracies, a distribution coefficient algorithm to deal with small nodal masses in MPM  was proposed by Ma et al. (2010), and applied in this study in combination with the proposed GC-SRI-patch method. Following the distribution coefficient algorithm, a part of the force acting on a node with a small mass is transferred to neighbouring nodes with a larger mass, so that mass and momentum conservation laws continue to be fulfilled. However, in the authors' opinion, such an algorithm cannot fully resolve acceleration inaccuracies near domain boundary nodes. Moreover, the same issue is particularly problematic in two-phase coupled problems, as it leads to spurious pore pressure oscillations in the presence of nearly incompressible pore water. For this benchmark, it was observed that pore pressure increments Δp t w,gp at the GPs were usually very small (i.e., Δp t w,gp ⩽1 × 10 − 5 kPa), but occasionally very large values (i.e., Δp t w,gp ⩾ 1.0 × 10 − 3 kPa) occurred for some grid cells at the domain boundary, especially those cells not containing MPs but influenced by the particle support domain within the framework of GIMP. These spurious large increments are directly related to large nodal accelerations caused by the small nodal mass issue that typically occurs near the domain surface. These large pore pressure increments can lead to inaccurate pore  pressure recovery at a small number of surrounding MPs (less than 1% of the total number of MPs). Without special treatment, these inaccurate pore pressures may propagate to the whole problem domain and result in a misleading pore pressure distribution across all MPs. For this reason, if the calculated pore pressure increment Δp t w,gp of a grid cell at the domain boundary is larger in absolute value than a specified threshold ζ d (for this benchmark, ζ d = 1.0 × 10 − 3 kPa), then it is set to zero. Even though this treatment was rarely used in this analysis and influenced only a small number of MPs, it was found to be generally effective in suppressing spurious pore pressure oscillations originating from MPs near the domain boundary.
Figs. 14a-e show how pore water drainage takes place gradually through the permeable boundaries of the block and promotes mechanical deformation of the solid skeleton during consolidation. As is apparent in Figs. 14a-c, the pore pressure dissipation is not monotonic in time, an occurrence associated in 2D problems with the so-called Mandel-Cryer effect (Mandel, 1953;Cryer, 1963). The same characteristic is more clearly illustrated in Fig. 15, where the time evolution of the excess pore pressure at two selected MPs (A and B) is plotted for both GIMP and GC-SRI-patch simulations. Stable GIMP results are only available up to shortly after the end of the gravity ramp loading. Conversely, GC-SRIpatch provides stable results for the entire duration of the hydromechanical analysis. In addition, it can be observed that the results of GC-SRI-patch with a much coarser mesh are quite similar to the solutions provided in Zhao and Choo (2020). The deviatoric stress (defined as , where σ 1 , σ 2 and σ 3 are principal stresses) distributions shown in Fig. 16 demonstrate the applicability of CMPM in coupled problems, and further validate the performamce of the proposed GC-SRI-patch. Fig. 17 displays the evolution in time, during and after the gravity ramp, of the excess pore pressure at the middle section of the slumping block (i.e., along the column of MPs highlighted in Fig. 15a). It is worth noting that the proposed GC-SRI-patch method captures correctly the gradual build-up of pore pressures during the ramp loading, as well as the downward propagation of a pore pressure wave when the increase in gravity is suddenly arrested at t = 0.1 s. Such propagation occurs simultaneously with global pore pressure dissipation, and is a natural outcome of the complete dynamic formulation.
To determine the influence of space discretisation on the numerical  solution, two additional space discretisations are used to simulate the consolidation of the slumping block (i.e., 400 and 2500 MPs initially placed on 100 and 625 square grid cells, respectively). Fig. 18 shows the calculated excess pore pressures at different times using GC-SRI-patch for both cases. The GC-SRI-patch returns stabilised solutions for both additional discretisations, with encouraging convergence performance upon mesh refinementcompare with the results in Fig. 14.

Conclusions
This study has presented an explicit, stabilised two-phase material point method named GC-SRI-patch for application in coupled poromechanical problems. The Generalised Interpolation Material Point (GIMP) method with a single set of MPs was adopted to alleviate cellcrossing errors and reduce the computational burden. To avert pore pressure instabilities, a Selective Reduced Integration (SRI) was used for the calculation of pore pressure increments at central GPs, which are of high(er) accuracy and convergence order. Such increments are then mapped to MPs using the proposed linear patch based on a Moving Least Squares Approximation (MLSA). Further improvement of effective stress recovery was achieved through the recently proposed Composite Material Point Method (CMPM), here applied for the first time to coupled two-phase problems. Other practical issues, including application of a surface traction on a movable boundary and the mitigation of 'small mass' issues near the domain boundaries, were also investigated.
Numerical verification examples supported the conclusion that the proposed GC-SRI-patch method can effectively be used to analyse relevant hydromechanical processes over a wide range of loading/drainage conditions. Instead of piecewise constant pore pressures over each cell, the proposed pore pressures return continuous distributions both within grid cells and at inter-cell boundaries, even in the presence of a coarse background grid. In particular, pore pressure instabilities were greatly mitigated by the new method, as is clearly demonstrated by the benchmark numerical solutions in terms of pore pressures and effective stresses. Future work will be devoted to further testing more challenging large-deformation analyses, more complex hydromechanical boundary conditions and more sophisticated constitutive models for the soil skeleton.

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.
The CMPM shape functions for a cell with only one neighbour cell (i.e., MPs located in a boundary cell) can be written as and are shown in Fig. 20. It should be mentioned that the above shape functions can only be used for a structured background mesh. In 2D/3D problems, shape functions are computed by multiplying the individual 1D functions in the different directions.