*Article* **Applicability of Constitutive Models to Describing the Compressibility of Mining Backfill: A Comparative Study**

**Ruofan Wang \*, Feitao Zeng and Li Li**

Research Institute on Mines and Environment (RIME UQAT-Polytechnique), Department of Civil, Geological and Mining Engineering, École Polytechnique de Montréal, C.P. 6079 Succursale Centre-Ville, Montréal, QC H3C 3A7, Canada; feitao.zeng@polymtl.ca (F.Z.); li.li@polymtl.ca (L.L.)

**\*** Correspondence: ruofan.wang@polymtl.ca

**Abstract:** The compressibility of mining backfill governs its resistance to the closure of surrounding rock mass, which should be well reflected in numerical modeling. In most numerical simulations of backfill, the Mohr–Coulomb elasto-plastic model is used, but is constantly criticized for its poor representativeness to the mechanical response of geomaterials. Finding an appropriate constitutive model to better represent the compressibility of mining backfill is critical and necessary. In this paper, Mohr–Coulomb elasto-plastic model, double-yield model, and Soft Soil model are briefly recalled. Their applicability to describing the backfill compressibility is then assessed by comparing numerical and experimental results of one-dimensional consolidation and consolidated drained triaxial compression tests made on lowly cemented backfills available in the literature. The comparisons show that the Soft Soil model can be used to properly describe the experimental results while the application of the Mohr–Coulomb model and double-yield model shows poor description on the compressibility of the backfill submitted to large and cycle loading. A further application of the Soft Soil model to the case of a backfilled stope overlying a sill mat shows stress distributions close to those obtained by applying the Mohr–Coulomb model when rock wall closure is absent. After excavating the underlying stope, rock wall closure is generated and exercises compression on the overlying backfill. Compared to the results obtained by applying the Soft Soil model, an application of the Mohr–Coulomb model tends to overestimate the stresses in the backfill when the mine depth is small and underestimate the stresses when the mine depth is large due to the poor description of fill compressibility. The Soft Soil model is recommended to describe the compressibility of uncemented or lightly cemented backfill with small cohesions under external compressions associated with rock wall closure.

**Keywords:** mining backfill; compressibility; constitutive models; numerical modeling

### **1. Introduction**

Backfill is being considered as an integral part of several underground mining methods. It is used as working platform in overhand cut-and-fill mining method or for creating safer working space in underhand cut-and-fill mining method. Using mine waste as underground mining backfill helps to minimize the surface disposal of mine waste [1–4]. However, the main objective of backfilling the mined-out spaces is to effectively control the rock wall closure and maintain the regional ground stability [5–9]. The compressibility of backfill plays an important role in resisting the closure of surrounding rock walls associated with adjacent extraction or/and creep behavior. Previous studies showed that a backfill of low compressibility can carry significant stresses generated by walls convergence and provide considerable support to surrounding rocks [6,10–13]. In underground mines, especially with overhand cut-and-fill or open stoping methods, however, the commonly used backfill is uncemented or has a low-cement content. The compressibility is large under low compression state and decreases as the compression increases. Understanding

**Citation:** Wang, R.; Zeng, F.; Li, L. Applicability of Constitutive Models to Describing the Compressibility of Mining Backfill: A Comparative Study. *Processes* **2021**, *9*, 2139. https://doi.org/10.3390/pr9122139

Academic Editor: Haiping Zhu

Received: 25 October 2021 Accepted: 22 November 2021 Published: 26 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and properly describing the compressibility of mining backfill is thus of specific interest for mining industry to evaluate fill performance and stability of underground structures.

Numerical modeling provides an efficient and cost-effective method to study the complex mechanical behavior of backfill. Nonetheless, the reliability and applicability of a numerical model largely depends on the capability of applied constitutive model. There are many constitutive models of geomaterial proposed in literatures with various levels of complexity [14–17]. Until now, the Mohr–Coulomb elasto-plastic model is the most used one to simulate the mining backfill due to the simplicity and clear physical meaning of model parameters [18–25]. The justification of the Mohr–Coulomb model is usually attributed to the good fit with the shear strength of backfill [26–30]. However, it is wellknown that the Mohr–Coulomb elasto-plastic model suffers from its linearly elasticity and the neglect of volumetric yielding. In reality, geomaterial can have a nonlinear behavior before yielding while a mining backfill can become denser upon a large compression generated by wall closure. Given the restrictions of the Mohr–Coulomb elasto-plastic model, it remains unknown which constitutive model should be applied to better represent the compressibility of mining backfill.

There are few research studies devoted toward identifying a constitutive model applicable to describing the compressibility of mining backfill. Among these studies, Oliver and Landriault [31] investigated the convergence resistance of backfill by simulating an oedometer test of dense sand with the Mohr–Coulomb model and the strain-softening model. Different values of Young's modulus and Poisson's ratio were applied. Results show that numerical model remains elastic over the full strain range except when a very small Poisson's ratio is applied. The predicted compressive stresses obtained with the two constitutive models are almost identical and significantly smaller than the experimental results. Clark [32] reproduced the non-linear stress–strain response of the dewatered tailings backfill in uniaxial compression tests with a cap model in FLAC. The input parameters for the cap model were obtained by applying the curve fitting on all the experimental results. The good agreements between the numerical model and experimental results do not mean that the calibrated numerical model can be used to correctly predict the mechanical behavior of the backfill under an untested stress condition. Fourie et al. [33] performed a finite element analysis by making use of the linearly elastic, Mohr–Coulomb, Drucker-Prager, and modified Cam-Clay models to simulate a backfilled stope at depth of 2000 m below the ground surface. The hanging wall convergence for the modified Cam-Clay model was found to be 11% larger than the results obtained with the other three models due to the plastic volumetric strain of backfill. The numerical results were not further compared with physical test data. Lagger et al. [34] used the double-yield model in FLAC*3D* to simulate the oedometer test of pea gravel as a filling material. The cap pressure of the double-yield model was calibrated based on the experimental results in a stress range of 0 to 6 MPa. Within this range, numerical results reasonably correlate with the test data, but the comparisons for the higher load stage was not shown. Therefore, more works are necessary to identify a suitable constitutive model to describe the compressibility of mining backfill, particularly by analyzing its predictive capability and comparing with physical results.

In this study, the Mohr–Coulomb elasto-plastic, double-yield, and Soft Soil models were recalled. Their abilities to describe the compressibility of backfill were assessed by comparing the numerical results obtained using FLAC*3D* with the experimental results of one-dimensional consolidation and consolidated drained triaxial compression tests of lowly cemented backfill available in the literature. Some unknown model parameters were calibrated based on part of the experimental results, and the calibrated models were applied to predict the other part of the experimental results. The identified model was then benchmarked with the Mohr–Coulomb model in simulating a backfilled stope, overlying a sill mat at different mine depths. The applicability of the identified model and the significance of modeling fill compressibility will be shown and discussed. In addition, validations of FLAC*3D* against analytical solutions of a cylinder hole in the linearly elastic

and Mohr–Coulomb material and the sensitivity analyses of the numerical models are provided in the Appendices A–C.

#### **2. Commonly Used Constitutive Models in Geotechnical Engineering**

For the sake of completeness, a few constitutive models commonly used in geotechnical engineering, including the Mohr–Coulomb elasto-plastic, double-yield, and Soft Soil models, are briefly recalled. Compression stresses are positive and tension is negative. All strength parameters are in terms of effective stresses.

#### *2.1. Mohr–Coulomb Elasto-Plastic Model*

The Mohr–Coulomb elasto-plastic model considers a material as linearly elastic and perfectly plastic once the stress state reaches a state of yield [35]. It is the most commonly used constitutive model in modeling the mechanical behavior of mining backfill.

Figure 1 shows the envelope of the Mohr–Coulomb model in *p*–*q* space (Figure 1a) and the typical stress–strain relation (Figure 1b). In the figure, *φ* and *c* are the friction angle and cohesion, respectively; *ε<sup>q</sup>* is the deviatoric strain; *p* and *q* are the mean and deviatoric stresses, respectively, expressed as:

$$p = \frac{(\sigma\_1 + \sigma\_2 + \sigma\_3)}{3} \tag{1}$$

$$q = \sqrt{\frac{\left(\sigma\_1 - \sigma\_2\right)^2 + \left(\sigma\_2 - \sigma\_3\right)^2 + \left(\sigma\_3 - \sigma\_1\right)^2}{2}}\tag{2}$$

where *σ*1, *σ*2, and *σ*<sup>3</sup> are the major, intermediate and minor principal stresses, respectively.

**Figure 1.** Schematic presentation of the Mohr–Coulomb elasto-plastic model: (**a**) yield surface in *p*–*q* space; (**b**) stress–strain relationship.

The linear stress–strain relation in the elastic regime (below the envelope) is described using a constant Young's modulus *E* and Poisson's ratio *ν* and assumed to follow Hooke's law. Once the stress state reaches the yield envelop defined by the Mohr–Coulomb criterion, infinite plastic shear strain occurs under a constant load.

The stress–strain relation in the elastic regime is expressed as:

$$
\varepsilon\_x = \frac{1}{E} \left[ \sigma\_x - \nu \left( \sigma\_y + \sigma\_z \right) \right] \tag{3}
$$

$$
\varepsilon\_{\mathcal{Y}} = \frac{1}{E} \left[ \sigma\_{\mathcal{Y}} - \nu (\sigma\_{\mathcal{X}} + \sigma\_{\mathcal{Z}}) \right] \tag{4}
$$

$$
\varepsilon\_z = \frac{1}{E} \left[ \sigma\_z - \nu \left( \sigma\_x + \sigma\_y \right) \right] \tag{5}
$$

$$\gamma\_{xy} = \frac{2(1+\nu)}{E} \mathfrak{r}\_{xy}, \ \gamma\_{yz} = \frac{2(1+\nu)}{E} \mathfrak{r}\_{yz}, \ \gamma\_{xz} = \frac{2(1+\nu)}{E} \mathfrak{r}\_{xz} \tag{6}$$

where *σx*, *σy*, *σz*, *τxy*, *τyz*, *τxz* are the components of stress tensor; *εx*, *εy*, *εz*, *γxy*, *γyz*, *γxz* are the components of the strain tensor.

The well-known Mohr–Coulomb criterion is expressed as follows, to relate shear strength *τ* and the corresponding normal stress *σ* [36,37]:

$$
\pi = \sigma \cdot \tan \phi + c \tag{7}
$$

or as follows, in terms of principal stresses:

$$\frac{\left(\sigma\_1 - \sigma\_3\right)}{2} - \frac{\left(\sigma\_1 + \sigma\_3\right)}{2} \cdot \sin\phi - c \cdot \cos\phi = 0\tag{8}$$

A 3D generalization of Equation (8) in terms of stress invariants is given as [38]:

$$q + \left(\frac{3\sin\phi}{\sin\theta\_l\cdot\sin\phi - \sqrt{3}\cos\theta\_l}\right)p + \left(\frac{3\cos\phi}{\sin\theta\_l\cdot\sin\phi - \sqrt{3}\cos\theta\_l}\right)c = 0\tag{9}$$

where *θ<sup>l</sup>* is the Lode angle.

The Mohr–Coulomb criterion correlates well with the shear strength of backfill, but tends to overestimate the tensile strength [39]. A tensile strength *T* smaller than that calculated by applying the equation of Mohr–Coulomb, called tension cut-off and, thus, usually applies for mining backfill. An associated or non-associated flow rule can be applied by varying the value of dilation angle *ψ* to model the plastic volume change due to shearing. The Mohr–Coulomb model does not capture the plastic volumetric strain under isotropic compression.

#### *2.2. Double-Yield Model*

The double-yield model was built in FLAC [40], which involves shear and tensile yield criteria of the Mohr–Coulomb model, and a volumetric yield surface. Its stress–strain relation in the elastic region is described by Hooke's law. Figure 2 illustrates the schematic yield surface and stress–strain behavior of the double-yield model.

**Figure 2.** Schematic presentation of the double-yield model: (**a**) yield surface in *p*–*q* space; (**b**) stress–strain relationship.

The volumetric yield surface (or cap) of the double-yield model shown in Figure 2a is independent on the deviatoric stress, and is defined as:

$$p - p\_{\mathcal{L}} = 0 \tag{10}$$

where *pc* is the current cap pressure (or preconsolidation pressure).

The double-yield model has an associated volumetric flow rule and its hardening rule relates to the plastic volumetric strain *ε p <sup>v</sup>* through a defined piecewise-linear function. The prescribed piecewise-linear function is flexible, but needs to be calibrated based on the results of physical tests. The bulk modulus *K* in the double-yield model is proportional to the derivative of cap hardening function as:

$$K = R \frac{\mathbf{d}p\_c}{\mathbf{d}\varepsilon\_v^p} \tag{11}$$

where *R* is a constant.

Equation (11) indicates that the elastic modulus of the double-yield model is dependent on a piecewise-linear function of the cap pressure, which explains the varied slope of the stress–strain curve in the elastic region, as shown in Figure 2b. Compared with the Mohr–Coulomb model, the double-yield model involves a volumetric yield surface, which enables accounting the plastic volumetric strain due to the mean stress. The doubleyield model has been adopted in some studies to simulate the mechanical performance of considerably compressible backfill material [34,41].

#### *2.3. Soft Soil Model*

The Soft Soil model is an advanced Cam-Clay type model [14,42] based on the critical state concept [43] and captures the irreversible void change accompanying the soil deformation. Figure 3a shows the relation between the volumetric strain *ε<sup>v</sup>* and mean stress *p* in the Soft Soil model. It is postulated that *ε<sup>v</sup>* linearly reduces with the increase of *p* along a normal consolidation line (NCL) in the semi-logarithmic space. The NCL has a slope of *λ*\*. For unloading and reloading, *ε<sup>v</sup>* varies following an elastic swelling line (SL) with a slope of *κ*\*. In the figure, *p*<sup>0</sup> is a reference value of mean stress. *ε<sup>n</sup> <sup>v</sup>* is the reference volumetric strain corresponding to (*p*<sup>0</sup> <sup>+</sup> *<sup>c</sup>*·cot*φ*) on the NCL and *<sup>ε</sup><sup>s</sup> <sup>v</sup>* is the reference volumetric strain corresponding to (*p*<sup>0</sup> + *c*·cot*φ*) on the SL. The equations for NCL and SL in the Soft Soil model are given as:

$$\varepsilon\_{\upsilon} = \varepsilon\_{\upsilon}^{n} - \lambda^{\*} \cdot \ln \left( \frac{p + c \cdot \cot \phi}{p\_0 + c \cdot \cot \phi} \right) \tag{12}$$

$$
\varepsilon\_v = \varepsilon\_v^s - \kappa^\* \cdot \ln \left( \frac{p + c \cdot \cot \phi}{p\_0 + c \cdot \cot \phi} \right) \tag{13}
$$

**Figure 3.** *Cont*.

**Figure 3.** Schematic presentation of the Soft Soil model: (**a**) logarithmic relation between the volumetric strain and mean stress; (**b**) yield surface in *p*–*q* space; (**c**) stress–strain relationship.

Figure 3b,c show the yield surface and stress–strain relation of the Soft Soil model. The yield surface consists of the envelope of the Mohr–Coulomb model and an elliptical cap. The elliptical cap has an apex on a critical state line as shown in Figure 3b and its formulation is expressed as:

$$\frac{q^2}{M\_s^2(p+c\cdot\cot\phi)} + (p-p\_c) = 0\tag{14}$$

where *Ms* is the slope of the critical state line in *p*–*q* space, which can be calculated based on the flow rule of the modified Cam-Clay model [42,44] as:

$$M\_{\rm s} = 3\sqrt{\frac{(1-K\_0)^2}{\left(1+2K\_0\right)^2} + \frac{(1-K\_0)(1-2\upsilon)\left(\frac{\lambda^\*}{\kappa^\*}-1\right)}{(1+2K\_0)(1-2\upsilon)\frac{\lambda^\*}{\kappa^\*}-(1-K\_0)(1+\upsilon)}}\tag{15}$$

where *K*<sup>0</sup> is the coefficient of earth pressure at-rest in normally consolidated condition.

The Soft Soil model employs an associated flow rule for the volumetric yield surface. The hardening of the volumetric yield surface is attributed to the plastic volumetric strain *ε p <sup>v</sup>* and is defined as:

$$\frac{\mathrm{d}p\_{\mathcal{c}}}{\mathrm{d}\varepsilon\_{\upsilon}^{p}} = \frac{p\_{\mathcal{c}}}{\lambda^{\*} - \kappa^{\*}}\tag{16}$$

The elastic modulus in the Soft Soil model is mean stress-dependent expressed as:

$$K = \frac{p + c \cdot \cot \phi}{\kappa^\*} \tag{17}$$

The Soft Soil model can be considered as a combination of the Mohr–Coulomb criterion and the modified Cam-Clay model. The critical state line in the Soft Soil model controls the shape of yield surface while the shear strength is defined by the Mohr–Coulomb envelope. The Soft Soil model has been used in a few studies to analyze the compressibility of soft clay [45,46], but has rarely been applied for mining backfill.

#### **3. Comparisons between Numerical Models and Laboratory Tests**

The applicability of constitutive models presented above to describing the compressibility of mining backfill, is identified by modeling one-dimensional consolidation and consolidated drained triaxial compression tests made on lowly cemented backfill available in the literature using FLAC*3D* [40].

#### *3.1. Comparison with One-Dimensional Consolidation Tests*

Pierce [26] conducted one-dimensional consolidation tests on Golden Giant paste backfill. The backfill samples were made by mixing the tailings with 3% binder by weight and cured for 28 days. Samples were casted in a rigid metal cylinder, which also acted as a confining ring in the consolidation tests. The backfill samples have a diameter of 75 mm and a height of around 37.5 mm. The measured properties include a density *ρ* of 2013 kg/m3, a porosity *n* of 49%, a cohesion *c* of 40 kPa and a friction angle *φ* of 41◦. During the consolidation tests of Pierce [26], a porous stone was put under the cylinder to allow drainage and a platen was placed on the top of the cylinder for the incremental load. Figure 4 shows the physical model and the laboratory results of the applied stress– compressive strain curve of one-dimensional consolidation tests conducted by Pierce [26]. The loading path has four stages, involving an increase from 0 to 4 kN in an increment of 0.5 kN, a decrease from 4 to 1 kN in an increment of 1 kN, an increase from 1 to 6 kN in an increment of 1 kN, and an increase from 6 to 12 kN in an increment of 2 kN.

**Figure 4.** Physical model and stress–strain curve of one-dimensional consolidation tests on Golden Giant paste backfill with a binder content of 3% and a curing time of 28 days conducted by Pierce [26].

Figure 5 shows a numerical model of a one-dimensional consolidation test constructed with FLAC*3D*, which involves a backfill sample, a cylinder (cell) and a loading platen. The platen is built to allow loading on a displacement boundary of the top surface. The numerical model shown in Figure 5 has the identical dimensions as the physical model of Pierce [26]. The mesh size of the numerical model is 3 mm based on the sensitivity analyses. Both the cylinder and platen are modeled as linear elastic material with a Young's modulus of 200 GPa and a Poisson's ratio of 0.3. The Mohr–Coulomb, double-yield, and Soft Soil models are applied for the backfill sample to make comparisons. Interface elements are applied between the cylinder and backfill. The normal and shear stiffness of the interface are determined based on an equation recommended in the FLAC*3D* manual [40]. The interface friction angle *φ<sup>i</sup>* is taken as 2/3 of *φ*, while its cohesion *ci* is assumed equal to *c* of backfill [47,48]. The displacements on the bottom of the model are restricted and other boundaries are set free. The same loading path in Pierce [26] tests is applied in the numerical simulations while the average normal displacement on the top of platen is recorded to calculate the compressive strain for different applied stresses.

**Figure 5.** A numerical model of one-dimensional consolidation tests built with FLAC*3D*.

In the numerical simulations, the Poisson's ratio of backfill is related to its friction angle through *ν* = (1 − sin*φ*)/(2 − sin*φ*) by considering a unique value of at-rest earth pressure coefficient *K*<sup>0</sup> [49,50].The tensile strength *T* of backfill is taken as 1/10 of its uniaxial compressive strength (UCS) [39]. The initial value of void ratio *eini* is calculated based on the measured porosity. However, some parameters remain unknown including *K* and shear modulus *G* for the Mohr–Coulomb model, *R* and the piecewise-linear function of *pc* for the double-yield model, *λ*\*, *κ*\*, and *pc* for the Soft Soil model. These parameters are determined by calibrating the numerical results based on part of laboratory results of Pierce [26] associated with the loading paths 1 and 2 as shown in Figure 4. Table 1 summarizes all parameters applied for different constitutive models. Numerical models with the calibrated parameters are then called the calibrated models, which are further applied to predict the other part of laboratory results associated with loading paths 3 and 4.


**Table 1.** Parameters of different constitutive models applied for backfill in numerical simulations of one-dimensional consolidation tests with *ρ* = 2013 kg/m3, *φ<sup>i</sup>* = 27◦, *ci* = 40 kPa.

Note: *Kmax* and *Gmax* denote the upper limits of the bulk and shear modulus. *K* and *G* for the Mohr–Coulomb model, *R* and piecewise-linear function of *pc* for the double-yield model, *λ*\*, *κ*\*, and *pc* for the Soft Soil model are calibrated based on the experimental results.

> Figure 6 shows the comparisons between the laboratory results of one-dimensional consolidation tests reported by Pierce [26] and the numerical results by applying the Mohr–Coulomb (Figure 6a), double-yield (Figure 6b), and Soft Soil (Figure 6c) models for backfill. In Figure 6a, one sees the compressive strain of the Mohr–Coulomb model linearly

increases as the applied stress increases. For the unloading stage, the stress–strain curve of the Mohr–Coulomb model is almost parallel to that in the loading stage. The minor scatter between the curves of loading and unloading is attributed to the yield of fill-wall interface. However, the experimental strain shows nonlinear relation with the applied stress in the test while only a small component of compressive strain is reversible at the unloading stage. The poor agreement between numerical and laboratory results is explained as that the Mohr–Coulomb model simulates a constant elastic modulus and does not capture the volumetric yield of backfill.

**Figure 6.** Comparisons between the experimental results of one-dimensional consolidation tests reported by Pierce [26] and the numerical results by applying the (**a**) Mohr–Coulomb; (**b**) double-yield; and (**c**) Soft Soil models for backfill with parameters in Table 1.

Figure 6b shows that the calibrated results of the double-yield model correlate well with the laboratory results, but the predicted compressive strain steeply increases as the applied stress further increases. It is because that the prescribed piecewise-linear function of the cap pressure in the double-yield mode is calibrated, based on part of the laboratory results. The prescribed function is flexible and can result in a good fit between the numerical and test results. However, when the applied stress exceeds the range of prescribed function, the double-yield model demonstrates infinite plastic volumetric strain as shown in Figure 6b. The predictive capability of the double-yield model is thus limited when the test data for calibration are insufficient. For the Soft Soil model, Figure 6c illustrates that both the described and predicted numerical results agree well with the laboratory results. Minor difference between the described results of the Soft Soil model and test results is seen when the applied stress is smaller than 230 kPa. This is because that the cementation in the backfill increases its primary stiffness at a small stress level. As the applied stress increases, the cement bond yields as shown by a drop of fill stiffness in Figure 6c. The mechanical behavior of cemented backfill then approaches an uncemented condition. In the Soft Soil models, the large primary stiffness caused by the cement bond at the small stress level can be pseudo-simulated by using an overconsolidation state, though their mechanisms are different. The Soft Soil model is thus deemed capable of describing the compressibility of lightly cemented or uncemented backfill in a confined compression condition.

#### *3.2. Comparison with Consolidated Drained Triaxial Compression Tests*

Rankine [28] conducted consolidated drained triaxial compression tests on Cannington paste backfill. The backfill samples have a diameter of 38 mm, a height of 76 mm, a cement content of 2%, a solid content of 74% by weight, and were cured for 28 days. The density of backfill is 2091 kg/m<sup>3</sup> and the porosity is 51.2%. Figure 7 shows the physical model and deviatoric stress–strain curve of consolidated drained triaxial compression tests under confining pressures of 100, 200, 500 kPa performed by Rankine [28].

**Figure 7.** Physical model and deviatoric stress–strain curve under different confining pressures of consolidated drained compression triaxial tests on Cannington paste backfill, with a cement content of 2% and a curing time of 28 days performed by Rankine [28].

Figure 8 illustrates a numerical model of backfill sample built with FLAC*3D* in simulations of consolidated drained triaxial compression tests. The numerical model has same sizes as the samples of Rankine [28], while the mesh size of the numerical model is 2 mm based on the sensitivity analyses. The Mohr–Coulomb, double-yield, and Soft Soil models are applied for backfill. The normal displacements on the bottom of the numerical model are restricted. The initial state is modeled by applying the confining stress normal to the surface of the sample after which the displacement is reset to zero. A normal velocity of <sup>1</sup> × <sup>10</sup>−<sup>7</sup> m/step is then applied on the top surface to simulate the compression. The normal stress and the axial displacement on the top surface of the numerical model are recorded during the calculation.

**Figure 8.** A numerical model of consolidated drained triaxial compression tests built with FLAC*3D*.

In the simulations, the Poisson's ratio is related to the friction angle of backfill considering a unique value of the at-rest earth pressure coefficient *K*0. The ratio of tensile strength *T* to UCS of backfill is taken as 0.41 according to the experimental results of Rankine [28]. *eini* is calculated as 1.05 based on the measured porosity. Calibrations based on the laboratory results under confining pressures of 100 and 200 kPa are performed to obtain some unknown parameters involving *c* and *φ*, *K* and *G* for the Mohr–Coulomb model, *R* and the piecewise-linear function of *pc* for the double-yield model, *λ*\*, *κ*\*, and *pc* for the Soft Soil model. Table 2 shows material parameters used for numerical simulations. The calibrated numerical models are then applied to predict the laboratory results of Rankine [28] under a confining pressure of 500 kPa.


**Table 2.** Parameters of different constitutive models applied for backfill in numerical simulations of the consolidated drained triaxial compression tests with *ρ* = 2091 kg/m3.

Note: *c* and *φ*, *K* and *G* for the Mohr–Coulomb model, *R* and the piecewise-linear function of *pc* for the double-yield model, *λ*\*, *κ*\*, and *pc* for the Soft Soil model are calibrated based on the experimental results.

> Figure 9 shows the comparisons between the laboratory results of consolidated drained triaxial compression tests conducted by Rankine [28] and numerical results under different confining pressures by applying the Mohr–Coulomb (Figure 9a), double-yield (Figure 9b), and Soft Soil (Figure 9c) models for backfill. Figure 9a illustrates that the calibrated and predicted strength of the Mohr–Coulomb model are close to the laboratory results. However, the elastic modulus of the Mohr–Coulomb model is constant while

the stiffness of mining backfill increases as the confining pressure increases. The Mohr– Coulomb model thus largely underestimates the stress magnitude at a given strain under a confining pressure of 500 kPa. Meanwhile, it overestimates the strain at failure under a confining pressure of 500 kPa by predicting a value of 38% while the experimental result shows a value of 20%.

**Figure 9.** Comparisons between the experimental results (dash lines with points) of consolidated drained triaxial compression tests reported by Rankine [28] and the numerical results (solid lines) under different confining pressures by applying the (**a**) Mohr–Coulomb; (**b**) double-yield and (**c**) Soft Soil models for backfill with parameters in Table 2.

Figure 9b shows that the double-yield model can reasonably describe the laboratory results for the confining pressures of 100 and 200 kPa. However, the predicted strength and stiffness of the double-yield model under a confining pressure of 500 kPa are very different from the experimental results. It is explained as that the infinite volumetric plastic strain occurs once the applied stress exceeds the upper bond of prescribed piecewise-linear function of the cap in the double-yield model. The predictive capability of the double-yield model is thus limited. Figure 9c shows that the described and predicted results of the Soft Soil model reasonably agree with the laboratory results. Based on the comparisons between numerical results and laboratory tests, the Soft Soil model is identified superior to the Mohr–Coulomb and the double-yield model in describing the compressibility of mining backfill with slight cementation (or uncemented backfill). In order to further exhibit the applicability of the Soft Soil model, it will be benchmarked with respect to the Mohr– Coulomb model in modeling a typical backfilled stope overlying a sill mat at different mine depths.

#### **4. Simulations of Backfilled Stope Overlying a Sill Mat**

In underhand cut-and-fill mining, uncemented or lightly cemented backfill is used to fill the mined-out upper stope overlying a sill mat. During the extraction of underlying stope, the sill mat will act as an artificial roof, which makes the stress distribution within the overlying backfilled stope significant for its stability. Sobhi and Li [11] analyzed this problem with PLAXIS*2D* only using the Mohr–Coulomb constitutive model to simulate the backfilled stope. The compressibility of backfill under the rock wall closure associated underlying extraction was thus not properly considered by Sobhi and Li [11]. In this section, the problem of a backfilled stope overlying a sill mat at different mine depths *D* of 200 and 1000 m are numerically investigated with FLAC*3D*. Emphasis is placed on the comparisons between numerical results predicted by applying the Mohr–Coulomb and Soft Soil models. Figure 10 shows a physical model and a plane strain numerical model (*D* = 200 m) of the problem. The symmetry plane (*x* = 0) is taken into account by considering half of the model. The excavations have a width *B* of 6 m. The overlying stope has a height *H* of 10 m and is filled with uncemented backfill. A gap of 0.5 m is left on the top of the backfill to represent the poor contact between fill and stope roof. The sill mat has a height *Hs* of 3 m, while the underlying stope is 13.5 m in height. The domain size of the numerical model is a distance from the origin to the boundaries of the model. Based on the sensitivity analyses, the numerical model is constructed with the optimal domain and mesh sizes of 300 and 0.25 m.

**Figure 10.** (**a**) Physical model and (**b**) numerical model built with FLAC*3D*of an undercut below a sill mat with overlying backfill.

The rock mass and sill mat obey the Mohr–Coulomb model while the overlying backfill is modeled with different constitutive models. The rock mass is characterized by unit weight *γ<sup>R</sup>* = 27 kN/m3, Young's modulus *ER* = 42 GPa, Poisson's ratio *ν<sup>R</sup>* = 0.25, cohesion *cR* = 9.4 MPa, friction angle *φ<sup>R</sup>* = 38◦, and dilation angle *ψ<sup>R</sup>* = 0◦. The sill mat is characterized by unit weight *γ<sup>s</sup>* = 20 kN/m3, Young's modulus *Es* = 1.5 GPa, Poisson's ratio *ν<sup>s</sup>* = 0.3, cohesion *cs* = 5 MPa, friction angle *φ<sup>s</sup>* = 35◦, and dilation angle *ψ<sup>s</sup>* = 0◦. Table 3 provides the material parameters for the overlying uncemented backfill in which the same parameters are applied in the Mohr–Coulomb and Soft Soil models, where possible. The

shear strength parameters (i.e., *ci* and *φi*) of fill–rock interfaces are considered equal to those of backfill by assuming rough rock walls.


**Table 3.** Parameters of the Mohr–Coulomb and Soft Soil models applied for overlying uncemented backfill with unit weight *γ* = 18 kN/m3.

In the numerical model, the displacement along the third direction (*y*-axis) is constrained to simulate a two-dimensional plane strain condition. The top boundary of the numerical model is set free to simulate the ground surface while normal displacement is restricted on the lateral boundaries. For the bottom boundary, the displacements are constrained in all directions. Numerical simulations are conducted at mine depths *D* of 200 and 1000 m respectively. The lateral earth pressure coefficient *Kr* = 2 is employed by considering the typical stress regime of the Canadian Shield [51]. In numerical simulations, the overlying stope is excavated after obtaining the initial equilibrium state. The displacement is then reset to zero and overlying stope is sequentially backfilled with 1 m per layer. This is followed by excavating the underlying stope in one step to expose the sill mat. Figure 11 shows the displacement distributions in the numerical model at each step as references with the Soft Soil model at a mine depth of 200 m.

**Figure 11.** Distributions of the displacement with the Soft Soil model at a mine depth of 200 m for different simulation steps of (**a**) excavating the overlying stope; (**b**) backfilling the mined-out overlying stope; and (**c**) extracting the underlying stope.

Figure 12 shows the iso-contours of bulk modulus in the overlying backfill after placement by applying different constitutive models. In the figure, one sees the bulk modulus of the Mohr–Coulomb model is 250 MPa and is a constant independent on the stope height. The bulk modulus of the Soft Soil model around the middle height of the stope is around 250 MPa, but its value moderately increases along the height of the backfilled stope, because the elastic modulus of the Soft Soil model is mean stress dependent, as given by Equation (17).

**Figure 12.** Distributions of the bulk modulus in overlying backfill after placement simulated with the (**a**) Mohr–Coulomb and (**b**) Soft Soil models.

Figure 13 illustrates the variation of the vertical and horizontal stresses along the vertical central line (VCL) of the overlying backfill before the excavation of underlying stope. Results shown in Figure 13 are independent on different mine depths because the backfill is placed after the rock wall displacement (i.e., delayed placement). In the figure, one sees that both the vertical and horizontal stresses increase smoothly along the stope height while the arching effect is evident by comparing with the overburden stresses. The stress distributions in the backfilled stope prior to the underlying excavation by applying the Mohr–Coulomb and Soft Soil models are almost identical. At the lower part of the stope, the stress state of the Soft Soil model is slightly larger than that of the Mohr–Coulomb model. The results shown in Figure 13 agree well with the results reported by Sobhi and Li [11].

**Figure 13.** Variation of the (**a**) vertical and (**b**) horizontal stresses along the VCL of overlying backfilled stope before the underlying extraction.

Figure 14 shows the variation of vertical and horizontal stresses along the VCL of overlying backfill after excavating the underlying stope at a mine depth *D* = 200 m. Under the influence of rock wall closure associated with the underlying extraction, the vertical and horizontal stresses in the overlying backfill increase compared to the results shown in Figure 13. However, the vertical and horizontal stresses of the Soft Soil model are smaller than those of the Mohr–Coulomb model below the stope height of 2 m. The value of vertical stress at the bottom of overlying backfill is 256 kPa for the Mohr–Coulomb model and is 169 kPa for the Soft Soil model. The different results of two constitutive models are explained as that the Soft Soil model simulates plastic volumetric strain of backfill under the compression from rock walls. The backfill needs to be compacted with certain amount of compressive strain before it can sustain large compressive stress. This feature is not captured by the Mohr–Coulomb model, which can thus overestimate the stress state in the overlying stope when the rock wall closure is small at a shallow mine depth. Since the stability of sill mat largely depends on the stresses within the overlying backfill, using the Mohr–Coulomb model may further cause an inaccurate estimation on the required strength of sill mat.

**Figure 14.** Variations of the (**a**) vertical and (**b**) horizontal stresses along the VCL of overlying backfill after excavating the underlying stope at a mine depth of 200 m.

Figure 15 shows the variations of vertical and horizontal stresses along the VCL of overlying backfill after the excavation of underlying stope at a mine depth *D* = 1000 m. As the mine depth increases from 200 to 1000 m, the rock wall closure associated with underlying excavation becomes larger which increases the vertical and horizontal stresses within the overlying backfill. The stress distribution predicted by the Soft Soil model is similar to that of the Mohr–Coulomb model above the stope height of 6 m as shown in Figure 15. However, the stresses of the Soft Soil model rapidly increase as the stope height increases. At the bottom of overlying backfill, the vertical and horizontal stresses of the Soft Soil model reach 2.2 and 7.5 MPa, which are significantly larger than the values of 0.9 and 2.5 MPa as predicted by the Mohr–Coulomb model. The different results of two constitutive models shown in Figure 15 are attributed to that the Soft Soil model accounts the volumetric hardening and pressure-dependent behavior of backfill. With a significant rock wall closure at a large mine depth, the mining backfill demonstrates large volumetric plastic strain, during which it becomes harder with a large elastic modulus, resulting in an increase in stresses generated by rock wall closure. The Mohr–Coulomb model does not simulate the volumetric hardening of backfill, which, thus, underestimates the stresses in backfilled stopes when the walls convergence is significant at a large mine depth.

**Figure 15.** Variation of the (**a**) vertical and (**b**) horizontal stresses along the VCL of overlying backfill after excavating the underlying stope at a mine depth of 1000 m.

Numerical simulations have shown that the Mohr–Coulomb and the Soft Soil models predict similar stress distributions in an isolated backfilled stope when the rock wall closure is absent. However, the stresses within backfill simulated with two constitutive models can be different when closure of surround rock walls is applied. In this condition, the commonly used Mohr–Coulomb model can under- or overestimate the stresses due to the neglect of volumetric yield and pressure-dependent behavior of backfill. The Soft Soil model is deemed more applicable to simulating uncemented or lightly cemented backfill when its compressibility or closure resistance is a dominant factor.

#### **5. Discussion**

By comparing the numerical results against the experimental results of one-dimensional consolidation and consolidated drained triaxial compression tests, the Soft Soil model is identified as capable of describing the compressibility of mining backfill. Nonetheless, the Soft Soil model should not be applied for backfill with a very large cohesion. This is partially because that its elliptical yield surface crosses the *p*-axis at a value of *c*·cot*φ* on the left side and has an apex on the critical state line as shown in Figure 3b. It implies that the lower bond of preconsolidation pressure in the Soft Soil model is determined by the value of cohesion. When a large cohesion is applied, the backfill modeled with the Soft Soil model is over-consolidated with a significant preconsolidation pressure, which can be unrealistic in some conditions. Another reason is that the slope of the critical state line (Equation (15)) in the Soft Soil model is calculated based on the flow rule of the modified Cam-Clay model, which considers a nil cohesion [42,44].Therefore, the Soft Soil model is deemed applicable to modeling the uncemented or lightly cemented backfill with a small (or nil) cohesion value.

For cemented backfill, cementation increases its primary stiffness at low stress condition by binding together fill particles. Experimental results show that the cement bond can yield as the applied stress increases after which the mechanical behavior of cemented backfill approaches an uncemented condition [26,52]. The results in Figure 6c show that the effect of cementation on fill stiffness at low stress levels can be pseudo-simulated using an over-consolidated state in the Soft Soil model. However, one should note that the mechanisms of cementation and overconsolidation are different. More effort is needed to investigate the effect of cement content and curing time on the compressibility of mining backfill and incorporate it in a constitutive model [52–54].

In numerical simulations, the Poisson's ratio of backfill relates to the friction angle as *ν* = (1 − sin*φ*)/(2 − sin*φ*), which is based on a unique value of at-rest earth pressure coefficient *K*<sup>0</sup> [49,50]. Such equation is practical in numerical modeling with an elastoplastic model. Previous studies have proposed several forms of empirical equation to define the relationship between *ν* and *φ* [55,56]. More works are needed to investigate this aspect based on experimental results. Meanwhile, since the stress–strain curve of soil-like material is highly nonlinear, how to determine the Poisson's ratio of backfill based on the experimental results for numerical modeling is a problem, and needs to be studied in future works.

The Soft Soil model postulates a perfectly plastic behavior when the stress state reaches the Mohr–Coulomb yield line. The post-peak behavior of backfill is affected by the cement content and confining pressure level [27,29]. Experimental results show that mining backfill with large cement content demonstrates strain softening under small confining pressures at the post-peak stage. The large confining pressure can also result in a strain hardening behavior of mining backfill. The post-peak behavior of backfill is not analyzed in this study.

This study focusses on modeling the compressibility of backfill in the long-term condition while the pore water pressure is not considered. More effort is thus necessary to evaluate the hydraulic conductivity and the effects of pore water pressure and drainage condition on the compressibility of backfill in the short-term condition [57–59].

The simulations of a backfilled stope overlying a sill mat indicate the applicability of Soft Soil model and the significance of modeling fill compressibility. The commonly used Mohr–Coulomb model tends to under- or overestimate the stress states in a backfilled stope when the walls closure is applied due to the poor description of the fill compressibility. Although the capability of the Soft Soil model has been tested against some consolidation and triaxial tests, field measurements are still needed when available to make further verifications.

#### **6. Conclusions**

The Mohr–Coulomb elasto-plastic, double-yield, and Soft Soil constitutive models are recalled and evaluated for the applicability to describing the compressibility of mining backfill. Numerical results with different constitutive models in FLAC*3D* are compared with one-dimensional consolidation and consolidated drained triaxial compression tests made on lowly cemented backfill available in the literature. Part of the experimental results is used to calibrate some model parameters and the calibrated models are applied to predict the other part of the test results. Based on the comparisons, the Soft Soil model shows the satisfactory description of the experimental results while its prediction is also quite good. The prevalently used Mohr–Coulomb model demonstrates poor correlations with the experimental results due to the neglect of volumetric yield and pressure-dependent behavior of backfill. The double-yield model accurately describes the experimental results based on the calibration, but its predictive capability is limited when the test results are insufficient.

The comparisons between Soft Soil and Mohr–Coulomb models in simulating a backfilled stope overlying a sill mat at different mine depths show similar stress distributions when rock wall closure is absent. However, when the rock wall closure associated with the underlying extraction is applied, application of the Soft Soil model shows that the Mohr–Coulomb model tends to overestimate the stresses in backfill at a shallow mine depth and underestimate the stresses at a large mine depth due to the poor description of the fill compressibility. The Soft Soil model is recommended to describe the compressibility of uncemented or lightly cemented backfill with small cohesions under external compressions associated with rock wall closure.

**Author Contributions:** R.W.: conceptualization, numerical modeling, analysis, literature, writing, and editing of the original draft. F.Z.: editing of the original draft. L.L.: project administration, methodology, supervision, editing of the original draft. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by the Natural Sciences, and Engineering Research Council of Canada (RGPIN-2018-06902), Mitacs Elevate Postdoctoral Fellowship (IT12569 and IT12570), China Scholarship Council (201706420059), and industrial partners of the Research Institute on Mines and the Environment (RIME UQAT-Polytechnique; http://rime-irme.ca/ accessed on 24 November 2021). The authors are grateful for their support.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within the article.

**Acknowledgments:** The authors acknowledge financial support from the Natural Sciences and Engineering Research Council of Canada, Mitacs Elevate Postdoctoral Fellowship, China Scholarship Council, and industrial partners of RIME UQAT-Polytechnique.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **List of Symbols**



#### **Appendix A Validation of FLAC***3D* **against Analytical Solutions of Stresses and Displacements around a Cylinder Hole in the Linearly Elastic Material**

FLAC*3D* can be validated against analytical solutions of stresses and displacements around an infinite cylinder hole in the infinite linearly elastic material. The problem can be analyzed in a plane strain condition. The physical model of this problem is shown in Figure A1. The origin locates at the central point of the model. In the figure, *a* is the radius of the hole. Domain size is the distance from the hole to the model boundary. *r* and *θ* are the cylindrical coordinates. The model is characterized by Young's modulus *E* and Poisson's ratio *ν*.

**Figure A1.** Physical plane strain model of a cylinder hole in an infinite linearly elastic material.

When the infinite hole is subject to a stress field composed of *σ*<sup>∞</sup> *<sup>x</sup>* , *σ*<sup>∞</sup> *<sup>y</sup>* , *σ*<sup>∞</sup> *<sup>z</sup>* , *τ*<sup>∞</sup> *xy*, *τ*<sup>∞</sup> *xz*, *τ*<sup>∞</sup> *yz*, the analytical solutions for stresses around the hole are given as [60,61]:

$$\sigma\_r = \frac{\sigma\_x^{\infty} + \sigma\_y^{\infty}}{2} \left( 1 - \frac{a^2}{r^2} \right) + \frac{\sigma\_x^{\infty} - \sigma\_y^{\infty}}{2} \left( 1 + 3\frac{a^4}{r^4} - 4\frac{a^2}{r^2} \right) \cos 2\theta + \tau\_{xy}^{\infty} \left( 1 + 3\frac{a^4}{r^4} - 4\frac{a^2}{r^2} \right) \sin 2\theta \tag{A1}$$

$$\sigma\_{\theta} = \frac{\sigma\_{\chi}^{\infty} + \sigma\_{y}^{\infty}}{2} \left( 1 + \frac{a^2}{r^2} \right) - \frac{\sigma\_{\chi}^{\infty} - \sigma\_{y}^{\infty}}{2} \left( 1 + 3 \frac{a^4}{r^4} \right) \cos 2\theta - \tau\_{xy}^{\infty} \left( 1 + 3 \frac{a^4}{r^4} \right) \sin 2\theta \tag{A2}$$

$$\sigma\_z = -2\nu \left( \sigma\_x^{\infty} - \sigma\_y^{\infty} \right) \frac{a^2}{r^2} \cos 2\theta - 4\nu \tau\_{xy}^{\infty} \frac{a^2}{r^2} \sin 2\theta + \sigma\_z^{\infty} \tag{A3}$$

$$\pi\_{r\theta} = -\frac{\sigma\_x^{\infty} - \sigma\_y^{\infty}}{2} \left( 1 - 3\frac{a^4}{r^4} + 2\frac{a^2}{r^2} \right) \sin 2\theta + \tau\_{xy}^{\infty} \left( 1 - 3\frac{a^4}{r^4} + 2\frac{a^2}{r^2} \right) \cos 2\theta \tag{A4}$$

$$\tau\_{\theta z} = \left(-\tau\_{xz}^{\infty}\sin\theta + \tau\_{yz}^{\infty}\cos\theta\right)\left(1 + \frac{a^2}{r^2}\right) \tag{A5}$$

$$\tau\_{\overline{z}\mathcal{I}} = \left(\tau\_{xz}^{\infty}\cos\theta + \tau\_{yz}^{\infty}\sin\theta\right)\left(1 - \frac{a^2}{r^2}\right) \tag{A6}$$

where *σ<sup>r</sup>* is the radial stress; *σθ* is the tangential stress; *σ<sup>z</sup>* is the normal stress along third direction (*z*-axis); *τrθ*, *τθz*, and *τzr* are the shear stresses around infinite cylinder hole based on cylindrical coordinates.

The analytical solutions for displacements around the infinite cylinder hole were given by Li [62], as follows:

$$\mathcal{U} = \frac{1}{E} \left\{ \frac{\sigma\_x^{\infty} + \sigma\_y^{\infty}}{2} (1 + \nu) + \left[ -(1 + \nu) \frac{a^2}{r^2} + 4 \left( 1 - \nu^2 \right) \right] \left( \frac{\sigma\_x^{\infty} - \sigma\_y^{\infty}}{2} \cos 2\theta + \tau\_{xy}^{\infty} \sin 2\theta \right) \right\} \frac{a^2}{r} \tag{A7}$$

$$V = -\frac{(1+\nu)}{E} \left[\frac{a^2}{r^2} + 2(1-2\nu)\right] \left(\frac{\sigma\_\text{x}^\infty - \sigma\_\text{y}^\infty}{2} \sin 2\theta - \tau\_\text{xy}^\infty \cos 2\theta\right) \frac{a^2}{r} \tag{A8}$$

$$\mathcal{W} = \frac{2(1+\nu)}{E} \left( \mathfrak{r}\_{xz}^{\infty} \cos \theta + \mathfrak{r}\_{yz}^{\infty} \sin \theta \right) \frac{a^2}{r} \tag{A9}$$

where *U*, *V,* and *W* are the components of displacement in the directions of *r*, *θ*, and *z* (third direction), respectively.

Figure A2 shows the corresponding plane strain numerical model built with FLAC*3D*. In the numerical model, hole radius *a* = 1 m, *E* = 10 GPa, *ν* = 0.25. The applied stress components include *σ*<sup>∞</sup> *<sup>x</sup>* = 15 MPa, *σ*<sup>∞</sup> *<sup>y</sup>* = 10 MPa, and *τ*<sup>∞</sup> *xy* = 3 MPa. The displacement along the third direction (*z*-axis) is restricted. In order to ensure stable numerical results, the domain size and mesh size of the numerical model need to be determined based on the sensitivity analyses.

**Figure A2.** Plane strain numerical model of a cylinder hole in an infinite linearly elastic material built with FLAC*3D*.

The numerical results of radial displacement and tangential stress are obtained at point M in Figure A2. Figure A3 shows the variations of radial displacement and tangential stress at point M as functions of the domain size. The variation of numerical results reduces as the domain size increases from 1 to 20 m, and become stable when the domain size is larger than 10 m. Therefore, the optimal domain size is determined as 12 m that is 6 times the size of the hole.

**Figure A3.** Variations of (**a**) radial displacement and (**b**) tangential stress at point M as functions of domain size.

Figure A4 shows the variations of radial displacement and tangential stress at point M as functions of mesh size. The mesh size ranges from 1 to 0.01 m. The numerical results become stable when the mesh size is smaller than 0.1 m. Further reduction of the mesh size will not greatly change the results. Therefore, the optimal mesh size is determined as 0.05 m to ensure stable results.

**Figure A4.** Variations of (**a**) radial displacement and (**b**) tangential stress at point M as functions of mesh size.

Numerical simulations are then conducted by using the optimal domain and mesh sizes. Figure A5 shows the comparisons between the numerical results of *σr*, *σθ*, *τrθ*, *U, V* long *x*-axis with the analytical solutions. The good correlations between the numerical and analytical results validate the FLAC*3D*.

**Figure A5.** Comparisons between the numerical results and the analytical results of (**a**) *σr*, (**b**) *σθ* , (**c**) *τr<sup>θ</sup>* , (**d**) *V*, (**e**) *U* long *x*-axis around a cylinder hole in the linearly elastic material.

#### **Appendix B Validation of FLAC***3D* **against Analytical Solutions of Stresses and Displacements around a Cylinder Hole in the Mohr–Coulomb Material**

FLAC*3D* can be further validated against analytical solutions of stresses and displacements around an infinite cylinder hole in the infinite Mohr–Coulomb (MC) material. Figure A6 shows the plane strain physical model of the problem subject to the isotropic in-situ stress *P*0. The origin is on the central point of the hole. Due to axial symmetry, only one-quarter of the model is considered. In Figure A6, *a* is the radius of the hole. The model is characterized by Young's modulus *E*, Poisson's ratio *ν*, cohesion *c*, friction angle *φ*, and dilation angle *ψ*.

**Figure A6.** Physical plane strain model of a cylinder hole in an infinite MC material.

According to Salencon's analytical solutions [63], the radius of yield zone *R*<sup>0</sup> around the cylinder hole can be calculated as:

$$R\_0 = a \left(\frac{2}{K\_p + 1} \frac{P\_0 + \frac{q\_{MC}}{K\_p - 1}}{P\_{in} + \frac{q\_{MC}}{K\_p - 1}}\right)^{1/(K\_p - 1)}\tag{A10}$$

where *Pin* is the internal pressure; *qMC* = 2*c*·tan (45◦ + *φ*/2); *Kp* = (1 + sin*φ*)/(1 − sin*φ*).

The radial stress *σr*, tangential stress *σθ*, and radial displacement *U* in the elastic zone are given as:

$$
\sigma\_r = P\_0 - (P\_0 - \sigma\_{re}) \left(\frac{R\_0}{r}\right)^2 \tag{A11}
$$

$$
\sigma\_{\theta} = P\_0 + (P\_0 - \sigma\_{\text{rc}}) \left( \frac{R\_0}{r} \right)^2 \tag{A12}
$$

$$
\Delta U = \frac{R\_0}{2G} \left( P\_0 - \frac{2P\_0 - q\_{MC}}{K\_p + 1} \right) \frac{1}{r} \tag{A13}
$$

where *σre* is the radial stress at the elastic-plastic interface and is given as:

$$
\sigma\_{rc} = \frac{1}{K\_p + 1} (2P\_0 - q\_{MC}) \tag{A14}
$$

The stresses and radial displacement in the plastic zone are given as:

$$
\sigma\_{\!\!\!T} = -\frac{q\_{\!\!M\!C}}{K\_p - 1} + \left( P\_{\!\!iner{I}} + \frac{q\_{\!\!M\!C}}{K\_p - 1} \right) \left( \frac{r}{a} \right)^{(K\_p - 1)} \tag{A15}
$$

$$
\sigma\_{\theta} = -\frac{q\_{M\mathbb{C}}}{K\_p - 1} + K\_p \left( P\_{inner} + \frac{q\_{M\mathbb{C}}}{K\_p - 1} \right) \left( \frac{r}{a} \right)^{(K\_p - 1)} \tag{A16}
$$

*<sup>U</sup>* <sup>=</sup> *<sup>r</sup>* 2*G* (2*υ* − 1) *<sup>P</sup>*<sup>0</sup> <sup>+</sup> *qMC Kp* − 1 <sup>+</sup> (<sup>1</sup> <sup>−</sup> *<sup>υ</sup>*) *Kp* <sup>2</sup> − <sup>1</sup> *Kp* + *K<sup>ψ</sup> Pin* <sup>+</sup> *qMC Kp* − 1 *R*<sup>0</sup> *a* (*Kp*−1) *R*<sup>0</sup> *a* (*Kψ*+1) + (<sup>1</sup> <sup>−</sup> *<sup>υ</sup>*) *KpK<sup>ψ</sup>* + 1 *Kp* + *K<sup>ψ</sup>* − *υ Pin* <sup>+</sup> *qMC Kp* − 1 *r a* (*Kp*−1) (A17)

where *K<sup>ψ</sup>* = (1 + sin*ψ*)/(1 − sin*ψ*).

Figure A7 shows the plane strain numerical model built with FLAC*3D*. In the numerical model, hole radius *a* = 1 m, *P*<sup>0</sup> = 10 MPa. The MC material is characterized by *E* = 10 GPa, *ν* = 0.25, *c* = 1 MPa, *φ* = 35◦, *ψ* = 0◦. The displacement along the third direction (*z*-axis) is restricted. Normal displacement on the left boundary and the bottom are prohibited to simulate symmetric planes. Sensitivity analyses are then performed to determine optimal domain and mesh sizes to ensure stable numerical results.

**Figure A7.** Plane strain numerical model of a cylinder hole in an infinite MC material built with FLAC*3D*.

The numerical results of radial displacement and tangential stress are obtained at point M in Figure A7. Figure A8 shows the variations of radial displacement and tangential stress at point M as functions of domain size ranging from 1 to 30 m. The numerical results become stable when the domain size is larger than 12 m. To be conservative, the optimal domain size is determined as 20 m.

**Figure A8.** Variations of (**a**) radial displacement and (**b**) tangential stress at point M as functions of domain size.

Figure A9 shows the variations of radial displacement and tangential stress at point M as functions of mesh size. As the mesh size reduces from 1 to 0.005 m, the variation of numerical results decreases. The numerical results become stable when the mesh size reduces to 0.03 mm. Further reduction of the mesh size will not greatly affect the results. Therefore, the optimal mesh size is determined as 0.02 m to ensure stable results.

**Figure A9.** Variations of (**a**) radial displacement and (**b**) tangential stress at point M as functions of mesh size.

The optimal domain and mesh sizes are then used to conduct numerical simulations of the problem. Figure A10 shows the comparisons between the numerical results of stresses and radial displacement long *x*-axis with the analytical solutions. The good agreements between the numerical and analytical results indicate that the FLAC*3D* is validated.

**Figure A10.** Comparisons between the numerical results and the analytical results of (a) *σ<sup>r</sup>* and *σθ* , (b) *U* along *x*-axis around a cylinder hole in the MC material.

#### **Appendix C Sensitivity Analyses of Domain and Mesh Sizes in the Numerical Simulations**

Sensitivity analyses are conducted to determine the optimal mesh size for the numerical models of consolidation and triaxial tests in this study. For the numerical model of a backfilled stope overlying a sill mat, both optimal domain and mesh sizes are determined based on the sensitivity analyses.

For the one-dimensional consolidation simulation, the physical and numerical models are illustrated in Figures 4 and 5. The material parameters are provided in Table 1. Figure A11 shows the variation of compressive strain under an applied stress of 500 kPa with different constitutive models as a function of the mesh size. The values of mesh sizes

vary between 20 and 2 mm. The numerical results are considered stable when the mesh size is smaller than 5 mm.

**Figure A11.** Variation of compressive strain under an applied stress of 500 kPa in one-dimensional consolidation simulations as functions of the mesh size.

For the numerical simulations of consolidated drained triaxial compression tests, the physical and numerical models are shown in Figures 7 and 8, while the material parameters are provided in Table 2. Figure A12 shows the variation of axial stress under an axial strain of 5% and a confining pressure of 200 kPa with different constitutive models as functions of the mesh size. The values of the mesh sizes range from 10 to 1 mm. The numerical results are considered stable when the mesh size reduces to 3 mm.

**Figure A12.** Variation of axial stress under an axial strain of 5% and a confining pressure of 200 kPa in triaxial compression simulations as functions of the mesh size.

For the numerical model of a backfilled stope overlying a sill mat shown in Figure 10b, Mohr–Coulomb model is applied for backfill to conduct the sensitivity analyses at a mine depth of 1000 m. The material parameters are provided in Table 3. In the sensitivity analyses, two indicators are analyzed for different domains and mesh sizes. One is the total displacement of surrounding rock mass at Point A in Figure 10b (one corner of sill mat) after extracting the overlying stope. The other indicator is the horizontal stress after excavating the underlying stope at Point B in Figure 10b, which is at the middle height on the VCL of backfill. Figure A13 shows the variation of total displacement at Point A and horizontal stress at Point B as functions of domain sizes ranging from 35 to 550 m. The numerical results become stable when the domain size is larger than 200 m.

**Figure A13.** Variation of (**a**) total displacement at Point A after extracting the overlying stope and (**b**) vertical stress at Point B after excavating the underlying stope as functions of domain size.

Figure A14 shows the variation of total displacement at Point A and horizontal stress at Point B as functions of the mesh size, which ranges from 5 to 0.1 m. The numerical results become stable when the mesh size reduces to 0.5 m.

**Figure A14.** Variation of (**a**) total displacement at Point A after extracting the overlying stope and (**b**) horizontal stress at Point B after excavating the underlying stope as functions of the mesh size.

#### **References**

