Computer Modeling in Engineering & Sciences |

DOI: 10.32604/cmes.2022.019112

ARTICLE

A Flux Based Approximation to Simulate Coupled Hydromechanical Problems for Mines with Heterogeneous Rock Types Using the Material Point Method

1Institute of Mine Seismology and University of Tasmania, Kingston, 7050, Australia

2School of Natural Sciences–University of Tasmania, Sandy Bay, Hobart, 7001, Australia

3School of Engineering–University of Tasmania, Sandy Bay, Hobart, 7001, Australia

*Corresponding Authors: Gysbert Basson. Email: Gys.Basson@IMSeismology.org

Received: 03 September 2021; Accepted: 21 October 2021

Abstract: Advances in numerical simulation techniques play an important role in helping mining engineers understand those parts of the rock mass that cannot be readily observed. The Material Point Method (MPM) is an example of such a tool that is gaining popularity for studying geotechnical problems. In recent years, the original formulation of MPM has been extended to not only account for simulating the mechanical behaviour of rock under different loading conditions, but also to describe the coupled interaction of pore water and solid phases in materials. These methods assume that the permeability of mediums is homogeneous, and we show that these MPM techniques fail to accurately capture the correct behaviour of the fluid phase if the permeability of the material is heterogeneous. In this work, we propose a novel implementation of the coupled MPM to address this problem. We employ an approach commonly used in coupled Finite Volume Methods, known as the Two Point Flux Approximation (TPFA). Our new method is benchmarked against two well-known analytical expressions (a one-dimensional geostatic consolidation and the so-called Mandel-Cryer effect). Its performance is compared to existing coupled MPM approaches for homogeneous materials. In order to gauge the possible effectiveness of our technique in the field, we apply our method to a case study relating to a mine known to experience severe problems with pore water.

Keywords: Fluid; permeability; Material Point Method; mining; faults

Numerical modelling is an essential tool for understanding the induced stress changes in underground hard rock mines. Mines rely on these tools for numerous purposes, of which arguably the most important is the planning of mining operations that mitigate seismic risk. Since engineers cannot physically observe the solid rock mass, numerical techniques must be capable of providing reliable predictions of the complex behaviour of the material. Many of the available commercial software packages simulate models that only allow for a single solid phase and frequently neglect groundwater’s influence on the stress regime in the rock mass [1–4]. Studies investigating the effect of groundwater principally concentrate on mine dewatering or the environmental impact in coal mines [5–8]. The coupled effect of the pore water pressure on the stress in the solid of hard rock mines is often overlooked when undertaking numerical stress modelling experiments. For example, it is expected that the effective stress of rock decreases as the water pore pressure increases. Therefore, it means that a geological fault that is filled with water may become unclamped under high pore pressure, making it much more likely to slip and consequently more seismically hazardous. Experiments that exclude the effect of the fluid may misleadingly suggest that the fault is stable. It is this prospect that motivates much of the work in this paper, in which we present a new numerical tool to aid hard rock mining applications which takes explicit account of the fluid phase, and thus making our tool a multiphase solver. We emphasize that we use the term multiphase to refer to the solid and fluid phases of the different materials in our context. This contrasts with the term multiphase’s standard meaning in petroleum engineering, where it is typically used when discussing two different fluids such as oil and water.

A standard modelling technique that is often used in the mining industry is based on boundary elements. The Boundary Element Method (BEM) [9] can represent large-scale mining configurations but assumes a medium that is linear elastic. This restricts its use to heterogeneous rock situations, making it less suited to the simulation of more complex multiphase problems. On the other hand, Finite Element Method (FEM) based solvers tend to be very reliable for analysing a wide range of mining-related problems. It has the distinct advantage that it can be used to model complex mining geometries in heterogeneous mediums. Furthermore, many extensions of the basic FEM, like the control volume FEM (CVFEM) [10], can simulate the coupled effect of multiphase problems. However, FEM can struggle when handling large deformations, and complex remeshing techniques are often essential to deal satisfactorily with this effect. To alleviate this problem, we suggest and investigate the properties of another family of methods often referred to as the Material Point Method (MPM) [11]. This strategy is attractive because the MPM can describe large deformations without invoking expensive and intricate remeshing. In mining simulations, grid locations near excavations (empty spaces) are susceptible to undergoing large deformations if high stresses exist. Moreover, the MPM can construct layouts with complex geometries with relatively little effort; this property is a direct consequence of the fact that the technique relies on a regular grid and so avoids the remeshing complications typically encountered within FEM. Furthermore, it has the property that it can simulate coupled multiphase problems using a single solver, therefore by-passing the need for a hybrid approach. One drawback of the MPM is its computational cost since particle, as well as grid states, have to be accounted for. High-density grids of entire mining layouts may necessitate significant memory resources. Basson et al. [12] describe a way to alleviate this effect by making use of a layered MPM background grid to simulate large mining problems.

The original formulation MPM was proposed by Sulsky et al. [11]. It is a hybrid Lagrangian and Eulerian Particle-In-Cell method that uses Lagrangian material points (otherwise known as particles) to carry history-dependent properties and employs an Eulerian background mesh as a communication tool for transferring information between the particles. Physical properties are first interpolated onto the background mesh, then the equations of motion are solved, and finally, the updated states are interpolated back to the particles themselves.

The MPM has been successfully used to study many problems, and its application to geotechnical questions [13–15]. More recently, the MPM has been extended to account for multiphase problems that embody fully coupled soil-fluid interaction [16–19]. The treatment of multiple phases within hydro-mechanical geotechnical problems using a single solver is particularly interesting to our current study. This method, used in the investigations [20,21], eliminates the need for complicated hybrid approaches that would otherwise be necessary to handle the coupling of the various phases.

Many coupled MPM formulations use the Biot theory approach [22] in which the momentum balance equation for the solid phase is solved on the Eulerian grid, and a Darcy velocity is assumed to describe the relative motion of the fluid with respect to the solid. The theory of mixtures [23] has also been successfully applied; this approach is founded on momentum balance equations solved independently on the Eulerian grid, while the coupling is modelled by introducing drag terms for both phases. Bandara [24] has successfully applied both the theory of mixtures and the Biot’s theory approach and used a two-point formulation to describe the coupled MPM (see Fig. 1). In the two-point formulations of the multiphase MPM, the medium consists of distinct groups of particles for each phase; solid particles and fluid particles are treated separately. This leads to the doubling up of particles within the grid and requires extensive computational resources. In the single-point formulation of the MPM, each particle carries information about both phases simultaneously.

Both Al-Kafaji [16] and Bandara et al. [17] have used a single-point formulation to solve problems within saturated porous media. Yerro et al. [19] extended the single-point formulation to allow for three phases in which solid, fluid and gas components are coupled. Each formulation has its advantages and disadvantages, and arguably none are demonstrably superior to the others in all instances. Instead, the choice of the particular implementation is best guided by the details of the specific problem at hand. For example, when a problem requires the presence of free water (see Fig. 2 as well as pore water, then the use of a two-point formulation based on Mixture Theory is recommended. However, attempting to apply this to quasi-static problems requires additional mechanical damping to both phases simultaneously, which may be difficult as this could lead to pressure oscillations [24]. It significantly increases the computational overhead associated with the requirement to account for both solid and fluid particles. If there is no free water present, then it is much easier to implement a single-point formulation as the position of the fluid is redundant and focus can be directed on the fluid pressure and velocities at the location of the solid particle instead. Therefore, for our study, the single-point MPM was chosen since our interest is to assess the response of the solid phase due to pore pressure changes in the fluid. Our investigations do not require any free water, which is ignored henceforth.

The remainder of the paper is structured as follows. In Section 2 we first outline the problems associated with using existing coupled MPM approaches to simulate materials with heterogeneous permeability. We then propose a modification that enables MPM to cater for mediums with variable permeability; this is the subject of Section 2.1. We then describe the details of a full computational cycle of our modified MPM. The performance of our new method is validated in Section 3 against two well-known analytical expressions and a comparison with the performance of a conventional coupled MPM as commonly found in the literature. Having benchmarked the new method against some test examples, we then use it to model an actual active mine in Section 4; in particular, we investigate the effect of pore pressure changes of a fluid-filled fault on the surrounding stress. We conclude the paper with some observations and final remarks.

In a single-point formulation of the two-phase MPM, the medium is represented by a solitary particle layer which means that each particle must carry information about both the solid and fluid phases. Most existing literature that uses this approach assumes that the continuum is defined as a homogeneous porous solid material with constant permeability and porosity [16,17,24]. This can be problematic when trying to simulate heterogeneous materials in which the permeability is heterogeneous. It is essential for what follows, so it is worth exploring this claim a little further. The schematic shown in Fig. 3, sets out, in simple terms, the main elements needed to advance the coupled MPM through one time step. For definiteness, let us suppose that the background grid is a regular Cartesian mesh; we need to bear in mind that what happens in this instance equally applies to more intricate nonstructured grids. The grid is initially populated with particles in some predetermined state, as depicted in Fig. 3a. The particle properties are then mapped onto the nearby cell nodes; this interpolation is performed using suitably chosen shape functions (Fig. 3b). In the original formulation of MPM, the shape functions are standard tri-linear forms similar to those often used in FEM methods. In more advanced versions of the implementation, higher-order shape functions can be used if greater transfer tiers are required outside the host cell (see, for example [25], or a variant of the MPM known as the Generalized Interpolated Material Point Method (GIMP) [26]). This intentionally smears information over a much larger volume than the cell occupied by the particle; this is a complicating factor when simulating nonhomogeneous materials as usually a clear boundary between different materials cannot be distinguished. Several adaptions that purport to address this problem in the MPM have been proposed [27–30]. In Fig. 3c, the grid is stepped forward in time by solving the equations of motion at the grid nodes. Once complete, the particles are updated in Fig. 3d. These include position, velocity and stress, among others. Finally, in Fig. 3e, the grid is reset in preparation for the next time step.

To delve a little deeper into the source of the problem, we focus on the process by which the fluid phase properties are mapped onto the nearby nodes, see Fig. 3b. First, we define the generalised Darcy’s law [31–33] for a fluid phase

where

The summation is taken over all Np particles that are located in the cells that share this node and depends upon the gradient of the shape functions

where n is the solid phase porosity and

so that the total nodal fluid phase mass is expressed as

On the other hand, when the material is heterogeneous, the situation changes as

2.1 Application of the Two Point Flux Approximation within the MPM

When using domain discretisation approaches to solve partial differential equations, a relatively large domain is divided into smaller elements, otherwise known as finite volumes (FVs). Finite Volume methods are usually the industry-preferred method for solving problems involving porous media as they are based on the conservation of continuity of fluxes [35]. These methods can be divided into two main classes, namely Two Point Flux Approximations (TPFA) [36] and Multi Point Flux Approximations (MPFA) [37]. TFPA methods are widely used in commercial software packages to study fluid flow in porous media; it is an attractive approach because it is simple to implement and computationally efficient. However, these methods struggle to accurately resolve the flux for grids that are not K-orthogonal; that is, instances when the principal directions of the permeability tensor do not coincide with the orientation of the input grid. The standard TPFA must be adapted to deal with such cases [38]. On the other hand, MPFA methods are adept at solving problems involving unstructured grids. Furthermore, as mentioned previously, when dealing with the full permeability tensor, one may choose to solve the problem using either FEM, Mimetic Finite Difference or the MPFA [39]. However, these implementations are typically much more complicated than the TPFA, which leads to a computationally expensive algorithm. In what follows, we attempt to adapt the TPFA to fit into the computational scheme of the MPM, and since the standard MPM grid is K-orthogonal, it follows that the TPFA is the preferred choice. If, however, unstructured grids are to be considered, then the MPFA should be preferred.

The total fluid flux through the surface of a finite volume can be written as

where

where the sum is taken over the six faces of the cube. The integral term on the right-hand side of Eq. (6) is usually referred to as the flux through face

To demonstrate how the flux at the interface of two cells are calculated in practice, consider two cells A and B, each with length l, as depicted in Fig. 5. Suppose that cell A has permeability k1, cell B has permeability k2 and that the shared face q of the two cells is a distance l/2 from the centre of each cell. The flux between the centre of A and q can be written

where pq is the fluid pressure at the interface and pA is the fluid pressure at the centre of A. An analogous expression for the other cell gives

where pB is the fluid pressure at the centre of B. Since the fluxes

where Tj = kj/l with j = 1, 2. If we use Eq. (9) in expression (7) we obtain the flux expression

which shows that the flux through q is a function of the harmonic average of the permeabilities of the two cells. If we were to apply Eq. (10) directly within the MPM framework, we would need to conduct additional calculations to determine the flux across all the cell interfaces in the background grid. Although the MPM scheme can be modified to include this step, it is undesirable as it adds extra computational overhead to the algorithm. Furthermore, since the nodes of the background mesh are used directly in the method, it would be advantageous to resolve the Darcy velocities at the nodes of the background mesh rather than across the interfaces of the grid cells. A schematic of conventional mapping scheme is depicted in Fig. 6 and requires the direct application of Eq. (2) to each one of the nodes denoted by the black squares. The Darcy velocity depends on the pore pressure and conductivity of all particles in the cells that share the node. As already mentioned, this method is only accurate if the conductivity of the particles contributing to the node is isotropic.

To improve matters, here we propose a modification of Eq. (2) to ensure that each node only uses information gleaned from particles located in cells that share a common interface. This will allow fluid flow to only occur across these boundaries, as is the case in the TPFA. To explain the process, we direct the reader to the second diagram in Fig. 6. Here we have drawn nine square cells, with the central one denoted C. Now the goal is to set the Darcy velocities at each of the four nodes of C such that it describes the total fluid flux through the surfaces of C. For definiteness, let us focus on the upper-left (north-west) node of C, which we have labelled Q; this choice of the node is not essential, and analogous considerations apply for the other three nodes of C.

We base many of our design decisions guided by the strategy employed in the usual TPFA. First, we remark that in examining node Q, we will only consider contributions from particles located in cells 1, 3 and 5; that is, cell C together with those cells that have Q has a node and which share an interface with C (i.e., the cells to the north and the west of C). For each of these, the cell pore pressure pc is set to be the average value of the individual particle pressures so that

where Np denotes the number of particles contained in the cell and

where

The standard shape functions of an isoparametric hexahedron can be calculated using the schematic idea that is sketched in Fig. 7. The illustration shows how to compute the shape function for a one-dimensional problem, and this is then easily extended to three dimensions by calculating the product of three one-dimensional functions

Next, the conductivity term is evaluated independently for each cell that shares an interface with C and has Q as one of its nodes. In this example, there are two such cells, and we denote these interfaces

With these key changes in mind, we can now derive the appropriately modified version of Eq. (2). To start the process, it is convenient to imagine that

where the terms Vi, pi and

As has been kindly pointed out to us by an anonymous referee, if

With these changes taken into account we can modify (14) to

where

We have developed result (16) with the specific node Q in mind. However, this balance can easily to extended so as to apply for any of the nodes of cell X: specifically

where subscript C represents a value evaluated at cell C and subscript c represents a value evaluated at an adjacent cell that shares a common interface with C as well the the node.

Given this modified expression for the Darcy velocity, we need to embed it within the cycle of calculations that are required to complete a time step within the modified MPM.

In this section, we will explain how we deploy the dynamic MPM in our modelling experiments. The Lagrangian particle layer present in our solid skeleton and the momentum balance equations for the fluid-solid mixture are solved on a regular Eulerian background mesh composed of trilinear shape functions S(x). Isothermal conditions hold as we shall assume that there is no mass exchange between phases, the solids are always fully saturated, and there is no acceleration of the fluid with respect to the solid. Moreover, we suppose that the densities of the two phases remain constant, the porosity of the solid phase does not change, the compressive stress is negative, while the pore pressure is positive.

Our medium is heterogeneous, which means we can not assume constant conductivity as commonly adopted in the majority of MPM applications [16,17,24]. We, therefore, need to carefully consider our treatment for determining the Darcy velocity of the fluid phase at the nodes of the background grid in the presence of this complication.

To initiate our model, suppose that at time t a particular particle

With this particle state specified, we can derive:

1. the volume of the mixture

2. the fluid phase mass

3. the particle mass

The method then proceeds as follows. For each particle p, we can identify the particular grid cell that contains p and then links are created between p and the eight nodes of that cell. Reciprocal links from the nodes back to the particle are formed, and the combined information is stored in a bi-directional map

For each grid node l, we use the information within

respectively. The force exerted due to the solid-fluid mixture is

where q = [1 1 1 0 0 0]T, b is a body force (such as gravity) acting on the particle p, and ts is an external traction force acting on the solid phase of the particle. The linear momentum balance equation for the mixture is solved at the node using

where the frequency independent damping term is defined to be

and is added to enhance the numerical stability of the process and to encourage the convergence to a quasi-static solution [40]. The constant

and boundary conditions for

At this point in the process, our novel MPM formulation deviates from many of the established two-phase MPM strategies. In these standard versions, the nodal Darcy velocities are calculated by appealing to Eq. (2), but, as we have already discussed, this does not take into account any anisotropy in the conductivity of particles. We have shown that Eq. (17) captures this phenomenon, and we use this balance to execute nodal Darcy velocity calculations.

Having updated the state of each node, we can revert to the particles. For each particle, we use

while the strain increments for the solid phase are

The fluid phase strain increment for each particle requires the calculation of the average pore pressure and conductivity of the cell which are set to be

Each cell is then revisited to deduce the associated nodal Darcy velocities

and the pore pressure is advanced using

where Kf is the fluid bulk modulus. The appropriate constitutive relation can then be applied to update the particle stress. Once all particle information has been advanced, the whole system has been moved through one timestep, and the process can then be repeated.

The explicit nature of the proposed iteration scheme means that the timestep

Fortunately, for a simulation that is kept near quasi-static equilibrium, this stringent constraint on the timestep can often be relaxed somewhat. Basson et al. [12] used the maximum nodal unbalanced force ratio and maximum particle velocity of the MPM to determine whether a system is in approximate equilibrium. If it is, the time step can be safely increased by adapting Eq. (20) to

where Dt is a mass density scaling factor. Here we set Dt = 1 if model is not in equilibrium but set

Before using our technique in practical contexts, it is essential to validate its performance against well-known and standard examples. For this purpose, we selected two problems [24] for which analytical results exist. The case studies were specifically chosen to test our new MPM and compare it with the results of existing two-phase MPM implementations. These benchmarks assume the medium has homogeneous permeability and thus, unfortunately, will not test the performance of our method when handling heterogeneous cases since analytical expressions for these cases are very sparse. On the other hand, one of the key features of our technique is that its treatment of the fluid-phase coupling is very different to that implicit in other coupled MPM formulations. In this way, our tests compare the performance of our methodology against analytical knowledge and give a way of assessing its usefulness and advantages over existing MPM implementations. In our experiments, we examined the values of the two sums in Eq. (17)

in the two subsections that follow.

3.1 Geostatic Stress Distribution Due to Gravity

The first example taken from Bandara [24] explores the geostatic consolidation problem presented by a fully saturated one-dimensional soil column. The initial solid-phase effective stress and fluid phase pore pressure are set to be zero at time t = 0. Gravity is then applied to the material, and the system is allowed to evolve to quasistatic equilibrium. The final solution is known to be

where

Our experiment was conducted on a grid of dimensions

We remark that in the course of the calculation, we noted the presence of oscillations in both the pore pressure and effective stress profiles during the system’s evolution and before settling to equilibrium. This phenomenon is illustrated in Fig. 9 and has been widely documented in literature [42,43] for other coupled methods. Further probing of the system revealed that these oscillations arise owing to the presence of the drag term

In this second example, we attempt to validate the performance of our proposed technique by applying it to the standard Cryer sphere problem [44]. This simulation requires tracking the pore pressure in response to a constant increase in the pressure within the sample medium followed by a steady decrease back to zero. This type of behaviour is commonly referred to as the Mandel-Cryer effect [45]. Cryer developed an analytic expression for the pore pressure evolution at the centre of a poroelastic sphere subjected to external pressure. The sphere has no initial stress or pore pressure, and the fluid is allowed to drain freely from its surface. Formal analysis predicts that the pore pressure at the centre begins to increase, reaches the level of the applied pressure and then continues to rise to some peak value. Then, as the fluid is allowed to drain, the pore pressure decreases again. Since this problem is fully dynamic, it is a good test for assessing the accuracy of the component T1 in Eq. (29). Notice, however, that this problem is not purely concentrated on T1 as the drag term

An MPM model was constructed to investigate this effect. A two-dimensional half-circle can visualise the inherent symmetry of the whole three-dimensional sphere (see Fig. 10. The grid has dimensions

Fig. 11 shows a comparison between our numerical results and the analytical prediction of the pore pressure at the centre of the sphere. We have encouragingly good agreement with the analytic results as well as against previous simulations reported in the literature [24,46]. We notice that at some point, our solution deviates slightly from the analytical expressions. This is especially evident for the case where the Poisson’s ratio was set to 0.1, as shown in the top panel of Fig. 11. Further investigations suggested that this effect can be ascribed to the observation that the regular MPM background grid cannot faithfully capture the perfect boundary of the sphere. A grid made of rectangular cells can only approximate the smooth boundary of the disk; an effect that is inherently clear and which is emphasised in Fig. 10. This effect can be mitigated by making the cell sizes smaller.

4 A Case Study: The Ernest Henry Mine

Ernest Henry is a copper and gold producing mine situated approximately 40 km North-East of Cloncurry in Queensland, Australia. Initially, the mine design was open pit; however, it introduced underground operations; presently, the deepest part is at about 1.2 km. The mine is seismically active, and its largest recorded event had a moment magnitude of 2.1. Groundwater is a significant problem for engineers working at the mine as they constantly have to drain the deeper levels to prevent flooding. Even though the host rock contains some pore water, it is thought that most of the water is contained in faults that act as the fluid’s principal transport mechanisms. Thus the porosity and conductivity are significantly more than that of the host rock.

On 28 September 2019, the mine drilled an exploration hole in the deepest part of the mine. The hole encountered a geological fault, at which point fluid was ejected from the hole. Several hours later, a large seismic event occurred on the same fault, but several hundreds of metres higher and near development areas. A schematic of the hole and the event location is shown in Fig. 12. This sequence of events raised two particular questions:

• Was the event triggered by the reduction in pore pressure in the fault?

• If pore pressure is released in the fault, how does it affect the nearby geology?

The fluid that is trapped within the faults is at high pressure, but this does not imply that it necessarily moves at high speed. The high pressure arises owing to the overburden weight of the fluid, and the porous nature of the rock means that Darcy’s law is applicable. The Institute of Mine Seismology (IMS) was tasked to investigate these questions. An in-house developed version of the standard coupled MPM was used to explore the issues raised. IMS has successfully used their version of the MPM in various three-dimensional numerical studies owing to its ability to simulate large deformation, dynamic problems twinned with the relative ease of building complex mining geometries. It is the prospect of this occurrence at Ernest Henry mine that led to the investigation of a way to simulate the coupled solid-fluid interaction of an heterogeneous medium as it is assumed that the permeability/conductivity of the faults is noticeably different than that of the host rock. We assume that the host rock permeability is very low compared to the rock mass that describes the fault zones and so, to this extent, can be legitimately regarded as a heterogeneous permeability rock experiment. This is a circumstance in which current MPM implementations will likely fail, and we try to overcome this problem using our new method.

A high-resolution single-point coupled MPM model was constructed for the mine. The scheme incorporated an adaptive grid of the type used in the work of [12]; this has the advantage of using various cell sizes with larger ones in areas far from the mine/faults but allowing for finer resolution elsewhere. This refinement facilitates optimal use of computational resources; on the one hand, it dispenses with dense cell packing in regions of little activity while, on the other, it enables high accuracy to be retained near important zones. In our simulation, the largest cells were cubic of side 10 metres while the smallest cubes had an edge of 1.25 metres. The cells that overlapped with the faults were of this smallest dimension, and thus it is remarked that the faults are not thin cracks but instead have a volume consistent with the assumption that they represent intrusions that are infilled by porous rock.

Fig. 13 shows a schematic of the model as well as the material properties. Insitu stresses used in the simulation are listed in Table 3. The model was composed of three materials, one being the host rock while the others were located inside the two faults: the first, identified as fault 6, is the one that intersected with the drill hole, and the second fault, known as the Angry Man Fault, is oriented nearly parallel to fault 6 directly above it. In some places, the two faults are quite close, being within 20 metres of each other, and for numerical purposes, both were ascribed a width of 3 metres and assumed to be porous, so they provided the transport mechanism for the fluid. The total model consisted of approximately 6 million grid cells, with each cell consisting of eight regularly spaced particles.

Several assumptions were made in the model:

• All materials were assumed to be poroelastic, however, it was assumed that only the two faults were filled with water and that the host rock was dry. This means that fluid flow was confined to be only within the faults.

• Initial pore pressure in the faults was set to approximately 20 MPa. This is a typical value seen in practice, but precise pressure could not be used as it was not measured by the mine.

• Draining in the model started approximately 150 metres below the lowest mining level. To reduce the computational complexity, no restriction was placed on the fluid’s outflow rate at the location where the hole intersected the Fault 6.

• To reduce pressure oscillations caused by the drag term (see Section 3.1) the system was damped such that fd was non-zero in Eq. (21). This ensured that the system remained in a near-equilibrium solution.

Initially, the model was solved for quasi-static equilibrium by assuming no mining has taken place and using the material properties and stress conditions stated. It was assumed that all materials were dry, which set the rock mass’s undisturbed state. The model was then resolved for the quasi-static equilibrium by removing all particles from the system within the mining excavations. This step represents the current mining state, and it was again assumed that all materials were dry. These two initial steps are common in stress-strain modelling for mining simulations for single-phase problems. Several common stress parameters were selected to quantify the faults as:

• The major principal stress magnitude

• The normal stress to the fault plane

• The in plane shear stress

• The excess shear stress, otherwise known as the Coulomb stress [47], defined as

The excess shear stress (ESS) is often used in elastostatic modelling to determine whether a fault may be at risk of slipping within the context of the elastic model. Areas on the fault that display negative (or zero) ESS values are considered stable, whereas areas of positive ESS are liable to slip, or may have already slipped at a prior point. Fig. 14 shows the results of the ESS on both faults at the end of the initial phase. Positive values of ESS are observed on both faults near the location of the seismic event, showing that this area of the fault is at risk of slip even before draining the F6 fault. It is envisioned that since the model predicts that this area is near failure, only minor changes in stress in this area could be sufficient to act as a triggering mechanism for the seismic event. Therefore, for the next phase, we are interested in changes in the selected parameters on the faults, rather than the absolute values as used in Fig. 14.

The draining phase was initiated by modifying the particles that characterise the faults. It was assumed that these were fully saturated, and a 20 MPa pore pressure was assigned to each particle. To preserve quasi-static equilibrium, the diagonal terms of the particle stress matrix was reduced by 20 MPa, thus ensuring that the total particle stress remains unaltered and Eq. (19) is unchanged. The draining process was then commenced by identifying several particles close to the drain point location and gradually reducing their pore pressures to zero. These pore pressures then remained at zero throughout the remainder of the experiment. Since the simulation was maintained near quasi-static equilibrium, we applied the less stringent time-stepping scheme as discussed in Section 2.3 enabling the calculation to be completed in a reduced time.

Fig. 15 shows the pore pressure inside the fault at two instances, namely at 5 sec and 41 sec after the commencement of draining. Figs. 16 to 19 summarise the form of four of the key parameters evaluated at the faults. The left panels of each of these figures show the Angry Man’s changes, whereas the right insets show the changes on the F6 fault. The values of the major principal stress

The results as outlined in Figs. 16 to 19 illustrate the importance of taking pore water into account in numerical simulations for mines such as Ernest Henry. Simulations that do not allow for this may show that both faults are stable and unlikely to fail. We have shown this in Fig. 14. Our models suggest that pore water is indeed essential in assessing the stability of the two faults. Even though pore pressure reduction in the F6 improved its stability, one may imagine that the opposite could occur; if pore pressure were to increase, then the fault becomes less stable. Furthermore, the effect of reduction in pore pressure may affect surrounding geological structures. In this case, the Angry Man fault conditions worsened, and even though the result was less pronounced than for the F6 fault, one may speculate that the reduction of pore pressure of the F6 fault could indeed be a triggering mechanism for the seismic event.

Video clips of the simulations of each modelling parameters used in this experiment can be obtained from:

• http://downloads.imseismology.org/numericalmodelling/mpm/ess.m4v Link to video clip showing Excess Shear Stress changes as the F6 drains

• http://downloads.imseismology.org/numericalmodelling/mpm/shear-stress.m4v Link to video clip showing

• http://downloads.imseismology.org/numericalmodelling/mpm/sigma1.m4v Link to video clip showing

• http://downloads.imseismology.org/numericalmodelling/mpm/normal-stress.m4v Link to video clip showing

In each of the clips we show how the draining of the F6 fault changes the properties of both itself and the nearby Angry man fault. These clips are free to download and redistribute.

A fully coupled two-phase three-dimensional MPM has been developed to model stress changes caused by fluid-filled rock in mining. This has been accomplished by modifying existing coupled MPM approaches to incorporate the description of mediums with heterogeneous permeability. Two verification examples were used to test the accuracy of the numerical model, and good agreement was observed in both cases. These benchmarks provide confidence that our methodology has value when modelling real mining cases in which a seismic event is likely triggered due to pore pressures changes caused by draining of a nearby fault. In this experiment, we have shown that even though the effective stress changes on the fault were relatively small, it was nonetheless still sufficient to play a role as a provocative mechanism for the event.

The case study at Ernest Henry mine highlights the importance of accounting for the effect of pore water in mining simulations, an effect that is often neglected. However, as mining methods become increasingly complicated, so the need for more advanced multi-phase numerical simulations grows. This is especially true for simulations of mine tailings dams for which the presence of fluid can play a central role in determining the outcome of experiments.

Our method can benefit from further improvement in two areas; a more advanced soil model might be implemented for dealing with problems involving partial saturation, and the numerical fluctuations in pore pressure could be reduced when dealing with dynamic problems. Further research in both these aspects must be pursued.

An intriguing avenue for future work would be an investigation of the relationship between our formulation and the innovative technique known as peridynamics. This is a novel non-local method that has been suggested as a vehicle for simulating hydraulic fractures [48,49]. It would be of interest to see how our method might be adapted to hydraulic fracturing scenarios as this is typically a multi-phase problem. This can not be done immediately since our current technique assumes that the faults are of finite thickness intruded by a porous material as opposed to being thin cracks. Further work is required to see if our method can be extended to address this important issue.

In summary then, our new model appears to show promise for modelling multi-phase coupling mechanisms in realistic mining situations. What we have described here is a first step in this direction, and extensions to further real case studies is the obvious way forward. As the catalogue of simulations increases, this would enable us to predict those circumstances in which mine sites might benefit from our ideas. One particular attraction of our suggestions is the fact that they can be incorporated into conventional MPM methods with only marginal extra computational cost.

Acknowledgement: We would like to thank Ernest Henry mine for allowing the use of their data in our case study. We are grateful to the reviewers for numerous suggestions and additional references that have significantly improved the quality of the manuscript.

Funding Statement: The authors received no specific funding for this study.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

- Itasca (2019). Fast lagrangian analysis of continua. https://www.itascacg.com/software/FLAC.
- Itasca (2021). Distinct element modelling of joint or blocky material. https://www.itascacg.com/software/3DEC.
- Simulia (2019). Abaqus finite element analysis. http://www.3ds.com.
- Wiles, T. (2021). Map3D. https://www.https://www.map3d.com.
- Surinaidu, L., Gurunadha Rao, V., Srinivasa Rao, N., & Srinu, S. (2014). Hydrogeological and groundwater modeling studies to estimate the groundwater inflows into the coal mines at different mine development stages using modflow, Andhra Pradesh, India.
*Water Resources and Industry*,*7–8*, 49-65. [Google Scholar] [CrossRef] - Rapantova, N., Grmela, A., Vojtek, D., Halir, J., & Michalek, B. (2007). Ground water flow modelling applications in mining hydrogeology.
*Mine Water and the Environment*,*26(4)*, 264-270. [Google Scholar] [CrossRef] - Atkinson, L. C., Keeping, P. G., Wright, J. C., & Liu, H. (2010). The challenges of dewatering at the victor diamond mine in Northern Ontario, Canada.
*Mine Water and the Environment*,*29(2)*, 99-107. [Google Scholar] [CrossRef] - Szczepinski, J. (2019). The significance of groundwater flow modeling study for simulation of opencast mine dewatering, flooding, and the environmental impact.
*Water*,*11(4)*, 1-2. [Google Scholar] [CrossRef] - Brebbia, C. (1973). The boundary element method for engineers. London: Pentech Press Limited.
- Kandelousi, M. S., Ganji, D. D. (2015). Control volume finite element method (CVFEM). In: Hydrothermal analysis in engineering using control volume finite element method. San Diego, US: Elsevier Science Publishing Co., Inc.
- Sulsky, D., Chen, Z., & Scheyer, H. (1994). A particle method for history-dependent materials.
*Computational Methods. Applied Mechanical Engineering*,*118(1)*, 179-196. [Google Scholar] [CrossRef] - Basson, G., Bassom, A. P., & Salmon, B. (2021). Simulating mining-induced seismicity using the material point method.
*Rock Mechanics and Rock Engineering*,*54(9)*, 4483-4503. [Google Scholar] [CrossRef] - Anderson, A., Mikkel, S., & Vabbersgaard, L. (2010). Modeling of landslides with the material point method.
*Computational Geosciences*,*14(1)*, 137-147. [Google Scholar] [CrossRef] - Larroca, F. (2015). The material point method in slope stability analysis (Master Thesis). University PolitÃယnica de Catalunya, Delft University of Technology, Netherlands.
- Wang, B., Vardon, P., Hicks, M., & Chen, Z. (2016). Development of an implicit material point method for geotechnical applications.
*Computers and Geotechnics*,*71(1)*, 159-167. [Google Scholar] [CrossRef] - Al-Kafaji, I. (2013). Formulation of a dynamic material point method (MPM) for geomechanical problems (Ph.D. Thesis). University of Studgard, University of Stuttgart, Germany. DOI 10.18419/opus-496. [CrossRef]
- Bandara, S., Ferrari, A., & Laloui, L. (2016). Modelling landslides in unsaturated slopes subjected to rainfall infiltration using material point method.
*International Journal for Numerical and Analytical Methods in Geomechanics*,*40(9)*, 1358-1380. [Google Scholar] [CrossRef] - Bandara, S., & Soga, K. (2015). Coupling of soil deformation and pore fluid flow using material point method.
*Computers and Geotechnics*,*63(2)*, 199-214. [Google Scholar] [CrossRef] - Yerro, A., Alonso, E., & Pinyol, N. (2015). The material point method for unsaturated soils.
*Géotechnique*,*65(3)*, 201-217. [Google Scholar] [CrossRef] - Vermeer, P., Sittini, L., Beuth, L., Wieckowski, Z. (2013). Modeling soil-fluid and fluid-soil transitions with applications to tailings. Tailings and Mine Waste 2013. Alberta, Canada. http://www.mpm-dredge.eu/documents/Vermeer,%20Sittoni,%20Beuth,%20Wieckowski%20-%202013%20-%20Modeling%20soil-fluid%20and%20fluid-soil%20transitions%20with%20applications%20to%20tailings%20-%20Tailings%20and%20Mine%20Waste%202013.pdf.
- Fatemizadeh, F., Stolle, D., & Moornann, C. (2017). Studying hydraulic failure in excavations using two-phase material point method.
*Procedia Engineering*,*175(1)*, 250-257. [Google Scholar] [CrossRef] - Biot, M. (1941). General theory of three-dimensional consolidation.
*Journal of Applied Physics*,*12(2)*, 155-164. [Google Scholar] [CrossRef] - de Boer, R. (2003). Reflections on the development of the theory of porous media.
*Applied Mechanics Reviews*,*56(6)*, R27-R42. [Google Scholar] [CrossRef] - Bandara, S. (2013). Material point method to simulate large deformation problems in fluid-saturated granular medium. London: University of Cambridge. https://books.google.com.au/books?id=uXTeoQEACAAJ.
- Steffen, M., Wallstedt, P., Guilkey, J., Kirby, R., & Berzins, M. (2008). Examination and analysis of implementation choices within the material point method.
*Computer Modeling in Engineering & Sciences*,*31(2)*, 107-127. [Google Scholar] [CrossRef] - Bardenhagen, S., & Kober, E. (2004). The generalized interpolation material point method.
*Computer Modeling in Engineering & Sciences*,*5(6)*, 477-495. [Google Scholar] - Nairn, J. (2003). Material point method calculations with explicit cracks.
*Computer Modeling in Engineering & Sciences*,*4(6)*, 649-663. [Google Scholar] [CrossRef] - Guo, Y., & Nairn, J. (2006). Three-dimensional dynamic fracture analysis using the material point method.
*Computer Modeling in Engineering & Sciences*,*1(1)*, 11-25. [Google Scholar] [CrossRef] - Jairn, J. (2007). Numerical implementation of imperfect interfaces. Computational Material Science, (40), 525–536. DOI 10.1016/j.commatsci.2007.02.010. [CrossRef]
- Nairn, J. (2013). Modeling imperfect interfaces in the material point method using multimaterial methods.
*Computer Modeling in Engineering & Sciences*,*1(1)*, 1-15. [Google Scholar] [CrossRef] - Coussy, O. (2003). Problems of poroelasticity. In: Poromechanics, USA: John Wiley & Sons, Ltd.
- Zhang, Q., Yan, X., & Shao, J. (2021). Fluid flow through anisotropic and deformable double porosity media with ultra-low matrix permeability: A continuum framework.
*Journal of Petroleum Science and Engineering*,*200*, 108349. [Google Scholar] [CrossRef] - Zhao, Y., & Borja, R. (2020). A continuum framework for coupled solid deformation fluid flow through anisotropic elastoplastic porous media.
*Computer Methods in Applied Mechanics and Engineering*,*369(4)*, 113225. [Google Scholar] [CrossRef] - Abe, K., Bandara, S., & Soga, K. (2014). Material point method for coupled hydromechanical problems.
*Journal of Geotechnical and Geoenvironmental Engineering*,*140(3)*, 1-16. [Google Scholar] [CrossRef] - Khattri, S. K. (2007). Analyzing finite volume for single-phase flow in porous media.
*Journal of Porous Media*,*10(2)*, 109-123. [Google Scholar] [CrossRef] - Aziz, K., Settari, A. (1979). Petroleum reservoir simulation. London: Elsevier Applied Science.
- Aavatsmark, I. (2002). An introduction to multipoint flux approximations for quadrilateral grids.
*Computational Sciences*,*6(3)*, 405-432. [Google Scholar] [CrossRef] - Zhang, W., Kobaisi, M. (2016). A pseudo two point flux approximation for general, non-orthogonal quadrilateral grids with improved monotonicity. 15th European Conference on the Mathematics of Oil Recovery. Amsterdam, Netherlands. DOI 10.3997/2214-4609.201601846. [CrossRef]
- Salama, A., Sun, S., & Wheeler, M. F. (2014). Solving global problem by considering multitude of local problems: Application to fluid flow in anisotropic porous media using the multipoint flux approximation.
*Journal of Computational and Applied Mathematics*,*267(4)*, 117-130. [Google Scholar] [CrossRef] - Cundall, P. (1987). Analytical and computational methods in engineering rock mechanics. In: Distinct element models of rock and soil structure. London: Allen & Unwin.
- Courant, R., Friedrichs, K., & Lewy, H. (1967). On the partial difference equations of mathematical physics.
*IBM Journal of Research and Development*,*11(2)*, 215-234. [Google Scholar] [CrossRef] - Haga, J., Osnes, H., & Langtangen, H. (2012). On the causes of pressure oscillations in low permeable and low compressible porous media.
*International Journal for Numerical and Analytical Methods in Geomechanics*,*36(12)*, 1507-1522. [Google Scholar] [CrossRef] - Murad, M. A., & Loula, A. F. D. (1994). On stability and convergence of finite element approximations of Biot’s consolidation problem.
*International Journal for Numerical Methods in Engineering*,*37(4)*, 645-667. [Google Scholar] [CrossRef] - Cryer, C. (1963). A comparison of the three-dimensional consolidation theories of Biot and Terzaghi.
*Quarterly Journal of Mechanics and Applied Mathematics*,*16(4)*, 401-412. [Google Scholar] [CrossRef] - Verruijt, A. (2015). Theory of problems of poroelasticity(Unpublished). http://geo.verruijt.net/software/PoroElasticity2015.pdf.
- Wong, T. (1998). A numerical study of coupled consolodation in unsaturated soils.
*Canadian Geotechnical Journal*,*35(6)*, 926-937. [Google Scholar] [CrossRef] - Ryder, J., Jager, A. (2002). A textbook on rock mechanics for tabular hard rock mines. Johannesburg: Creda Communications.
- Wang, Y., Zhou, X., & Xu, X. (2016). Numerical simulation of propagation and coalescence of flaws in rock materials under compressive loads using the extended non-ordinary state-based peridynamics.
*Engineering Fracture Mechanics*,*163(4)*, 248-273. [Google Scholar] [CrossRef] - Zhou, X. P., & Wang, Y. (2021). State-of-the-art review on the progressive failure characteristics of geomaterials in peridynamic theory.
*Journal of Engineering Mechanics*,*147(1)*, 3120001. [Google Scholar] [CrossRef]

This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. |