Stiffness Determination of Backfill-Rock Interface to Numerically Investigate Backfill Stress Distributions in Mine Stopes

Numerical modeling is an effective and efficient method to investigate the stress distributions of backfill in stopes, which should be well understood in underground mining. Interface elements between backfill and rock in simulated stopes had been proved to be essential components, for which the stiffness parameters need to be assessed and assigned. However, few reports have revealed the effects of interface stiffness on backfill stress distributions, and there is not yet a clear solution to determine the interface stiffness to simulate stresses in backfilled stopes, except an empirical method for simply applying a high value suggested in FLAC manual. In this study, a new solution is first proposed to determine the normal stiffness and shear stiffness of interface elements, respectively, in numerical modeling of backfill stresses. *e applicability of the solution has been verified by investigating backfill stress distributions in mine stopes of two widely used mining methods with variable stiffness values. *e results show that the newly proposed method leads to totally the same backfill stress distributions with models applying the interface stiffness by the method in FLACmanual based on a “rule-of-thumb” but will save at least 20%–30% calculation time to improvemodeling efficiency under the same simulation conditions and will carry much clear physical meanings corresponding to the interaction between backfill and rock walls in mine stopes. In addition, the vertical and horizontal stresses show good agreements with the analytical stresses predicted by the Marston equation under the at-rest state, which validates the reliability of the proposed solution for interface stiffness. Moreover, the plotting methods of stress distributions and the coefficient of lateral earth pressure of backfill in simulated stopes with proposed interface stiffness were discussed to further clarify the reasonable methods to investigate the backfill stresses in mine stopes, especially after considering the effects of the convergence from rock walls, which is a very significant and common phenomenon in practical mining engineering.


Introduction
Mining with backfill is used to extract valuable minerals as much and safe as possible by filling underground voids with mine wastes (tailings, waste rocks) and/or binders. e backfill in stopes can provide steady working platform, increase ore recovery, and improve ground conditions in mining processes [1][2][3]. e consumption of the mine wastes for backfill is also beneficial for the environmental protection of mine communities, which has accelerated the development and utilization of the mining methods world widely [4][5][6].
Open stoping with subsequent backfill and upward cut and fill are the most widely used mining methods with backfill in metal mines. e first method is designed with large stopes for mass mining, leading to high efficiency and safe working environment. Backfill slurry is poured into the mined voids to build vertical pillars for the following extraction of adjacent stopes. It is the most efficient mining method with backfill, which is usually employed in excavating thick ore bodies (e.g. copper, iron). e second method is commonly used in mining of midthick to narrow and/or irregularly shaped ore bodies, which often host highvalue minerals (e.g. gold, silver). Normally, after cutting two layers of an upward height (about 4∼6m) in stopes, backfill slurry is poured about half height (2∼3 m) of the mined void to construct a platform for next upward cut. By looping through the process, the stopes are mined out with a high recovery ratio because it can make flexible decisions to excavate one area in every cut layer after identifying the minerals contents within the boundaries of the ore veins.
To successfully implement the desired functions of backfill in different mine stopes, backfill stress distributions should be well understood and evaluated, which will further construct the basis for required strength design and stability analysis of backfill in the above mining methods [7][8][9][10][11]. Arching theory [12,13] has been used to develop some analytical solutions to calculate backfill stresses in vertical or inclined stopes [14][15][16][17]. Some physical tests were carried out to validate the analytical results of backfill stresses [18,19]. In addition, numerical simulations typically with FLAC [20] and other programs are proved to be effective and efficient to investigate the stresses in backfilled stopes with sensitive analysis of influential factors, including stope dimensions and backfill strength properties. [21][22][23][24][25][26]. e convergence of rock walls in a stope can have variable effects on backfill stress distributions due to different interaction mechanism of backfill-rock contact surfaces in above two mining processes.
It is known that the backfill cured with tailings and/or binders is a soft material compared with the hard rock in mine stopes. e contacting surfaces should be simulated by interface elements to represent the particular behaviors of sliding, rotating, or separating the backfill from the rock walls. Recently, it had been proved that backfill stress distributions were dominated by physical (e.g. roughness) and mechanical (e.g. cohesion, internal friction angle) parameters of the contacting interfaces between backfill and rock walls rather than by shear strength properties of backfill [27,28], if the interface shear strength is lower than that of backfills. Direct shear tests further demonstrated that the shear strength parameters along the backfill-rock interfaces are proportional to those of the backfill [29][30][31]. erefore, it is very necessary to consider the interface elements between backfill and rock walls when numerically simulating stress distributions of backfill in mine stopes.
Apart from the shear strength parameters of interfaces, the numerical models built with FLAC (or other programs) also require stiffness parameters. In the FLAC manual [20], the elastic stiffness of slip and/or separation interfaces are not considered as important, and it is suggested to use an apparent stiffness recommended for glued interfaces based on a "rule of thumb." e physical meanings of the suggested empirical equation are unclear neither. A method of stiffness determination for backfill-rock interface with clear physical meanings is desired, while there are very few published reports in current literatures.
Kong et al. [32] analyzed interface stiffness between mudstone and concrete piles by the back-analyses method based on field measured shear stresses and displacements of the piles. A sufficiently high value for interface stiffness was suggested to match the field dilation effect of installed piles into the mudstone. A similar approach with a high value for interface stiffness had been adopted by Pirapakaran [33] to numerically analyze vertical stress and displacement of filling materials. In the study by Sivakugan et al. [34], interface stiffness were determined by following the recommended empirical equations from the FLAC manual to compare backfill stresses in different positions of stopes. In other numerical simulations with self-developed programs [22,35], the determination of interface stiffness was not presented in detail. ese works indicate that the stiffness along the contact interfaces between backfill and wall is significant in numerical models but all using very large values based on empirical methods. A definition of interface stiffness with clear physical meanings is needed to be proposed.
In this study, the empirical solution suggested in the FLAC manual to determine the apparent stiffness of interface elements will be firstly discussed. And then, new expressions with practical physical meanings for normal and shear stiffness of interface elements are proposed. To validate the reliability of the proposed method for interface stiffness, the vertical and horizontal stresses of the backfill in mine stopes of the two mining methods are simulated by FLAC3D after applying different interface stiffness values. In addition, the simulated backfill stresses with the proposed interface stiffness will then be verified against the analytical stresses with the Marston equation [12]. Furthermore, different plotting methods for backfill stresses and the coefficient of lateral earth pressure in simulated backfilled stopes of two mining methods with different interface stiffness values will be further discussed. e advantages of newly proposed solutions for interface stiffness in numerical modeling will be summarized at the end of the article.

Suggested Empirical Method in FLAC.
When interface elements are applied in numerical models with FLAC program [20], the normal stiffness (k n ) and shear stiffness (k s ) have to be assigned based on Mohr-Coulomb constitutive model, as shown in Figure 1.
ree types of interface elements are defined in the FLAC manual [20], including glued interface (artificial component to connect subgrids together), slip and/or separation interface (real interface that is stiff and can slip and/ or open under anticipated loading), and interface containing fill material (real interface that is soft to affect the system, such as a geological fault or joint with soft clay).
Without any doubt, the second type should be applied to the contacting interfaces between backfill and rock walls in mine stopes. For such type of interfaces, the FLAC manual [20] states that the strength properties of the interfaces are important, but the stiffness properties are not important. It further suggests using the stiffness solution recommended for glued interfaces based on a "rule of thumb,": where k neib.max is the equivalent stiffness of the stiffest neighboring zone, K (Pa) is the bulk modulus of the backfill, G (Pa) is the shear modulus of the backfill, and Δz min (m) is the size of the smallest backfill element at the contact zones of the interface.

Advances in Civil Engineering
Clearly, equation (1) shows that the normal (k n ) and shear (k s ) stiffness have a unit of "Pa/m." However, it is unclear whether the suggested utilization of an extremely large value for interface stiffness is just to ensure stable numerical results, and whether it can reflect the physical response of backfill-rock interaction in stopes or not. Besides, there is still no explanation why the manual argues that the stiffness properties of the interface elements are not important. So far there have been no reasonable explanations for these doubts.
To better understand this issue, one considers the classical definition of the stiffness k in the field of material mechanics. In the SI system (International System of Units) [36], stiffness is the resistance of an elastic body to deflection or deformation by an applied force, which can be detailed expressed as a ratio between the applied force and the corresponding displacement.
us, the stiffness k can be expressed as follows: where ΔF (N) is the variation of an applied force on an area A (m 2 ), which leads to a length variation ΔL (m) for the material, Δσ (Pa) is the variation of induced stress, and ε is elastic strain (relative deformation, dimensionless) of the material with an original length L (m). From equation (2), it is found that the stiffness k should have a unit of "N/m" or "Pa·m," which is different from the unit of "Pa/m" of the normal (k n ) and shear (k s ) stiffness in manual of FLAC [20]. e updated solutions of apparent stiffness should be consistent with methods in the manual about the stiffness unit to meet the requirements in following numerical simulations by the FLAC program.

Newly Proposed Solution.
Based on the comparison of definitions for interface stiffness, one can find the physical meanings of the normal stiffness k n characterizing the ratio of normal stress to displacement perpendicular to the interface elements and the shear stiffness k s characterizing the ratio of shear stress to displacement parallel to the interface elements, near the contacting surface between the backfill and rock walls in mine stopes, as shown in equations (3) and (4).
en, the normal stiffness k n and the shear stiffness k s of the stiffest neighboring zone along backfill-rock interface can be approximately defined as equations (5) and (6): where K (Pa) is the bulk modulus of the backfill, G (Pa) is the shear modulus of the backfill, and Δz min (m) is the size of the smallest backfill element at the contact of the interface. Equations (5) and (6) construct the newly proposed expression for determining the apparent normal and shear stiffness of the interface elements between backfill and rock walls in simulated backfilled stopes. e reliability of the proposed solutions and the induced influences on backfill stress distributions will then be compared with the cases using different interface stiffness values including those suggested parameters of FLAC manual [20] from the viewpoints of stably computational accuracy and efficiency. e simulated stresses will also be validated with the analytical results by the Marston equation [12].

Numerical Models and Cases.
To validate the proposed solution for apparent interface stiffness and investigate their influences on backfill stress distributions in stopes of two most widely used mining methods, numerical models were built with FLAC3D, which is a finite difference program using an explicit Lagrangian calculation scheme and a mixed discretization zoning technique. It is a well-adapted tool for handling geomechanical problems, including excavation and backfilling of mine stopes [3,8,9,15,19,[25][26][27][28]. Figure 2(a) shows a vertical backfilled stope with a width of B. e mined void of the stope is filled with two programs to a height of H to simulate the open stoping with subsequent backfill mining and the upward cut and fill mining methods.
For the first filling program, the whole height H of the stope is mined out with one excavating step, which simulates mass mining of open stoping with deep blasting holes. e energy release and closure deformation of surrounding rock walls caused by excavation are assumed to be fully completed before the backfilling, as it practically takes months to muck Advances in Civil Engineering out the blasted ores and then return to fill the mined void. In this program, an equivalent multistep loading for the whole height of the backfill was simulated with 10 layers [28,37] to obtain a stably static numerical solution based on a series of sensitivity analyses (not shown here) about the effects of filling layers on simulated backfill stresses, which revealed that the backfill stresses became numerically stable and insensitive to the further increase the filling layers after reaching 10. Besides, a space of 0.5 m in height was left at the top part of the backfilled stope, to simulate the practical difficulty of roof-contact filling due to drainage and settling of backfill slurry, or in some occasions, the space is designedly left to be a platform for continuous mining of upper stope.
For the second filling program, the ore rocks in the stope are excavated from bottom to top and layer by layer following the practical cut and fill process. In detail, after cutting the first and second layers (2 m per layer in height) starting from stope bottom, a backfill of 2 m in height will then be loaded to the first layer void. en the third layer will be excavated on the backfill platform of the first layer. By looping through the process, a stope of 40 m in height will be mined out with 20 layers' cutting and filling respectively, which also meets the requirements of using more than 10 filling layers for numerically stable results. e simulations of every cutting and filling were solved to reach the same equilibrium state in FLAC. In this program, the convergence of rock walls induced by one layer cut might be applied to the previous layers of backfill. e stress distributions of backfill will then be affected by its self-weight and the closing pressure of rock walls which is transferred from the interfaces between backfill and rock mass. Figure 2(b) shows a typical numerical model constructed with FLAC3D corresponding to the physical model. e bottom boundary is fixed in all directions, and the two lateral boundaries are only fixed in horizontal direction. All elements of the model are fixed in Y-direction (perpendicular to paper plane) to simulate plane strain conditions. After a series of domain and mesh sensitivity analyses, the external boundaries are determined at 120 m from stope boundaries, and the optimal mesh size for the backfill is set to 0.2 m. e rock mass is considered homogeneous, isotropic, and linearly elastic, characterized by c r � 27 kN/m 3 (unit weight), E r � 30 GPa (Young's modulus) and μ r � 0.25 (Poisson's ratio). e values were selected based on typical values given in the literature [38,39]. e backfill is modeled (1) Rock Mass -Elastic   [1,18,19,[40][41][42][43] and partly on authors' experience gained from mining projects with backfill. e mechanical properties of the interfaces between backfill and rock walls are characterized by the cohesion c i � 0 kPa and internal frictional angle δ � 0.6ϕ � 21°, which are based on the published shear tests of interfaces [29][30][31].
e normal (k n ) stiffness and shear (k s ) stiffness of interfaces are also key inputs for the numerical models. Table 1 shows the fixed parameters of rock mass, backfill, and interfaces in all models. It needs to be emphasized that the effects of backfill compaction and density-dependent stiffness of the backfill along the height of mine stopes were not considered in the numerical models here. e physical and mechanical parameters of the backfill were uniformly distributed and kept as constant in the simulated stopes shown in Table 1. Table 2 presents numerical programs to simulate and investigate the backfill stresses varying with key influential factors, which are emphasized on the normal (k n ) and shear (k s ) stiffness of interface elements. e backfill stresses in stopes of two typical mining methods and the stress plotting methods were also critical factors that should be discussed. e influence of a given parameter on the backfill stress distributions in mine stopes is analyzed by varying its value and keeping other parameters unchanged. Besides, what should be mentioned here is that the normal and shear stiffness for the interfaces between backfill and rock walls are kept unchangeable along the height of the backfilled stopes, which was inherited from the conventional hypothesis used in the FLAC manual [20] and relevant published articles [21,28,[32][33][34] without considering the cumulative effects of interface variability with filling depth in mine stopes.
Based on the above numerical models and backfill properties, one can have backfill bulk modulus K � 250 MPa and shear modulus G � 115 MPa calculated with backfill Young's modulus E � 300 MPa and Poisson's ratio μ � 0.3 and the smallest backfill element size Δz min � 0.2 m at the contacting surface to the interfaces.
With the empirical solution given in FLAC manual [20], the normal stiffness (k n ) and shear stiffness (k s ) of the interfaces are k n � k s � 20.2 GPa/m � 2.02 × 10 10 Pa/m after applying (1).

Effect of Interface Stiffness on Backfill Stress in Open
Stoping with Subsequent Backfill. e numerical results presented in Figure 3 shows that the normal stiffness (k n ) of the interface elements can significantly affect the backfill stresses along the VCL of the stope. e horizontal stress increases slightly with the increase in the interface normal stiffness, but the vertical stress decreases obviously when the interface normal stiffness increases from 2 to 200 MPa/m. When the normal stiffness is above 200 MPa/m, the vertical stress becomes almost insensitive to its further increase. What should be emphasized is that when k n � 2.02 × 10 10 Pa/m is set for interface normal stiffness after applying equation (1) of the empirical solution suggested in FLAC manual [20], the vertical and horizontal stresses follow the same distributions with the case using the interface normal stiffness of the newly proposed method here (k n � 1.25 × 10 9 Pa/m, k s � 5.8 × 10 8 Pa/m).
However, it is not computationally efficient to increase interface stiffness to higher values because the running time to reach the same equilibrium state will become undesirably higher under the same calculation conditions. In the numerical model with only one backfilled stope shown in Figure 2, the required running time to complete the calculation with the interface stiffness by the newly proposed method can be statistically 20%-30% less than the simulating backfill stress using interface stiffness of the solution suggested in FLAC manual. is computationally efficient effect will be exponentially highlighted if a lot of backfilled stopes are simulated in a model with interface elements between backfill and rock walls. erefore, it can be preliminarily concluded that the newly proposed equations to determine the interface stiffness carry more advantages in calculating efficiency of numerical simulation to obtain the same backfill stresses than the empirical solution suggested in FLAC manual; meanwhile, there have been clearer physical meanings corresponding to the interaction between backfill and rock walls in mine stopes. Figure 4 illustrates vertical (σ v ) and horizontal (σ h ) stresses along vertical central line (VCL) of stopes with different shear stiffness k s and fixed normal stiffness k n of interfaces in open stoping with subsequent mining method. Figure 4 shows that the vertical and horizontal stresses in the backfilled stopes are less sensitive to the variation of shear stiffness k s comparing with the normal stiffness k n shown in Figure 3. However, it presents the vertical stress decreases with the increase in shear stiffness of interfaces from the model of k s � 2.02 × 10 6 Pa/m to model of k s � 2.02 × 10 8 Pa/m, whereas the horizontal stress only presents a slight decrease. Besides, the same distributions about the vertical and horizontal stresses are shown for the models after applying the shear stiffness of k s � 2.02 × 10 10 Pa/m by equation (1) suggested by FLAC manual and the models applying the interface stiffness by the newly proposed method of equations (5)        Advances in Civil Engineering Figure 5 shows that the vertical stress reduces significantly with the normal stiffness of interface increasing from k n � 2.02 × 10 6 Pa/m to k n � 2.02 × 10 8 Pa/m, whereas the horizontal stress has an obvious trend of enlargement. However, the vertical and horizontal stresses become almost insensitive to the variation of the normal stiffness k n of interface when its value is above 2.02 × 10 8 Pa/m, even though the stresses in the model of k n � 2.02 × 10 8 Pa/m still have a bit difference with the stresses distributions in the model of k n � 2.02 × 10 10 Pa/m. It is worth mentioned that both vertical and horizontal stresses after applying interface stiffness by the newly proposed method (k n � 1.25 × 10 9 Pa/m, k s � 5.8 × 10 8 Pa/m) have the same numerically stable stress distributions with the model of k n � 2.02 × 10 10 Pa/m by the suggested empirical method in FLAC manual, but the newly proposed method took much less computation time (i.e., higher efficiency for simulation). Figure 6 illustrates the similar trends of vertical and horizontal stresses with the variations of shear stiffness k s compared with Figure 4. In details, the vertical stress along vertical central line (VCL) of backfilled stopes decreases gradually with the increase in the shear stiffness of interface from k s � 2.02 × 10 6 Pa/m to k s � 2.02 × 10 10 Pa/m, whereas the horizontal stress reduces slightly under these variations of shear stiffness. In addition, the stresses in the backfill model with interface shear stiffness calculated by suggested empirical method of FLAC manual show again the same and stable distributions with the backfill model of applying interface stiffness determined by the newly proposed method. e vertical and horizontal stresses varied with different normal and shear stiffness values shown in Figures 5 and6 demonstrate again the good reliability to determine reasonable interface stiffness parameters to achieve numerically stable results on stresses of backfilled stopes in upward cut and fill mining method. e numerical modeling about the stresses in backfilled stopes considering interfaces with stiffness values determined by the newly proposed solution has better calculating efficiency and clearer physical meanings to simulate the interaction between backfill and rock walls, compared with the numerical models of using interface stiffness by empirical solution suggested in the FLAC manual.

Comparison of Simulated Stresses with Analytical Results of the Marston Model.
e sensitive analyses about the simulated stress distributions in backfilled stopes with variable stiffness values of interface elements between backfill and rock walls shown in Figures 3-6 have demonstrated the strong reliability and stable accuracy of the newly proposed methods to determine the interface stiffness, together with much better calculating efficiency than the empirical solution suggested in FLAC manual. To further validate the newly proposed solution, the numerical results about backfill stresses in mine stopes of two typical mining methods will be compared with the analytical results obtained by the Marston solution [12,13], which is a classical theory to evaluate the backfill stresses under arching effect in mine stopes [14][15][16][17]. Figure 7 presents the numerically simulated vertical (σ v ) and horizontal (σ h ) stresses of the backfill along the vertical central line (VCL) of stopes in the open stoping with subsequent backfill mining method and the upward cut and fill mining method, respectively. In these numerical models, the normal (k n ) and shear (k s ) stiffness of interfaces were set up after applying the proposed solution (equations (5) and    Advances in Civil Engineering (6)), and the other parameters were referred to the values shown in Table 1. e analytical results with the Marston theory [12,13] by applying equations (7) and (8) to calculate the vertical and horizontal stresses are also plotted in Figure 7 for verifying the simulated stresses. In addition, the stress distributions based on the overburden solution are added in Figure 7 for comparison which has no consideration of arching effect for the backfill in stopes. e overburden backfill stresses can be assessed by σ v � c × h, σ h � K 0 × σ v , where h (m) is the depth from the top surface of the backfill, K 0 � 1-sinϕ is the at-rest coefficient of lateral earth pressure, and ϕ is the internal friction angle of backfill.
where c (� 18 kN/m 3 ) is unit weight of backfill, B � (10 m) is the width of backfilled stope, δ (� 0.6×ϕ � 21°) is internal friction angle along interfaces between backfill and rock walls in stopes, K earth is coefficient of lateral earth pressure of backfill that can be set to active state K earth � K a � (1-sinϕ)/ (1+sinϕ) or at-rest state K earth � K 0 � 1-sinϕ. It can be seen from Figure 7 that for the backfilled stopes in the open stoping with subsequent backfill mining method (red dots), the vertical and horizontal stresses along the vertical central line (VCL) of the backfill are much lower than the overburden stresses (black lines), which indicates the occurrence of arching effect, and these stress distributions show good agreements with the analytical stresses predicted by the Marston solution (blue lines) under the at-rest state (K earth � K 0 ). Based on the comparisons, it validates again the good reliability of the newly proposed methods to determine the interface stiffness values for numerically simulating backfill stresses in practical mining methods.
In addition, Figure 7 illustrates that for the backfilled stopes in the upward cut and fill mining method (purple triangles), the vertical stress is much lower, and the horizontal stress is much larger than the analytical stresses by the Marston solution (equations (7) and (8)). Furthermore, the simulated horizontal stress at the height of h � 0-20 m of the backfilled stope in the upward cut and fill mining method is even much larger than the overburden stress. ese differences mean a much stronger arching effect for the backfill stresses in the stopes of the upward cut and fill mining method. It can be preliminarily speculated that this phenomenon is mainly caused by the lateral convergent displacements of rock walls acting on the previous backfill layers during the later cutting and filling process. To further explain it, Figure 8 illustrates the horizontal convergent displacements along the rock walls in the stopes of the two mining methods.
From Figure 8, it can be seen that the convergent displacements along the right rock wall of the stope ( Figure 2, the same results can be obtained to plot along the left rock wall) in open stoping with subsequent backfill mining method are only slightly larger than zero (0-0.1 mm), which indicates a positive displacement going forward along the Xaxis. is is because the whole height of the stope was excavated at once step in the mining method, and the closure deformation of rock walls had finished before backfilling the stope void. e self-weight of the subsequent backfill will lead to a pushing pressure on the rock walls to make a slight increment of rock walls' extending displacement. e highest positive displacement is closing to the bottom part (height 30 m) of the stope and then gradually decreasing to zero at the top and bottom surfaces of the backfill.
However, Figure 8 shows that the convergent displacements along the right rock wall of the stope in the upward cut and fill mining method are much lower than zero (−0.4∼−1.3 mm), showing an approximately semicircular distribution with a diameter of stope height H (highest closure deformation near the midheight of stope and gradually decreasing towards the top and bottom surfaces). What is more interesting is that the convergent displacements are all negative values, indicating obvious convergence of rock walls leading to closure pressure on every backfill layer. is gives a good explanation of the much stronger arching effect (shown in Figure 7) in the backfilled stope of the upward cut and fill mining method and the nonagreement with the analytical stresses based on the Marston method.
Combined with the results in Figures 7 and 8, it not only further verifies the reliability of the newly proposed methods for interface stiffness determination to numerically simulate backfill stresses in mine stopes but also reveals that the backfill stress distributions are affected obviously by different levels of rock walls convergence based on the excavating and filling processes in different mining methods, which forms a practical factor to modify the arching theory originally from the field of soil mechanics [12,13].

Plotting and Monitoring Methods for Backfill Stresses.
A good understanding of the stress distributions in mine stopes is essential for underground mining design with backfill, but it is more important to reasonably evaluate the vertical and/or horizontal stresses in backfill to meet their practical usages in specific engineering conditions. For numerical modeling of stress distributions in backfilled stopes, there are currently two typical methods to plot and/ or monitor the backfill stresses (see details in Figure 2(a) and Table 2). e first method is to plot stresses along vertical central line (VCL) after filling the whole height of excavated stope void, which mainly illustrates the stress variations along with backfill height. Also, the second method is to monitor stresses at the midpoint of stope bottom (MPSB) after filling each layer in stopes, which usually presents the stress or pressure at stope bottom varying with the increasing process of backfill height. Similar discussions on the plotting and monitoring methods were reported by Sivakugan et al. [34] but only for the vertical stress analysis without considering horizontal stress. Based on its final conclusions [34], the second method was preferable to monitor the vertical loadings on the preburied conduits at the bottom of trenches with subsequent filling by granular materials. e preferable second method of monitoring the vertical stress in numerical models was also employed to compare with the laboratory physical tests of the vertical stress in dry granular fills [18,19,44], to investigate arching effects of sands behind retaining walls, fills in trenches, and hydraulic backfills in stopes.
However, in the practical engineering conditions of underground mining with backfill, not only the vertical stress but also the horizontal stress of backfill should be well evaluated. In most conditions, the horizontal stress distributions are relatively more important. For example, in order to assess the lateral pressure on adjacent structures, including the cemented backfill in adjacent primary stopes [3] or the man-made barricades in access drifts throughout the backfilling process [45], it is the horizontal stress that should be paid more attention, which will be transferred to horizontal pressure to affect the stability of these structures.
To illustrate the differences of the two methods to plot and monitor stresses, Figure 9 shows the contrastive results about the vertical and horizontal stress distributions plotted along vertical central line (VCL) and monitored at the idpoint of stope bottom (MPSB) with different normal stiffness k n of interfaces between backfill and rock walls in open stoping with subsequent backfill mining method.
It can be seen from Figure 9 that there are obvious differences about the backfill stress distributions between plotted along VCL and monitored at MPSB under the same numerical models. What is more interesting is that only at the top and the bottom surfaces of the backfilled stopes are the stress values that are same for the two plotting methods. Both the vertical stress and the horizontal stress plotted along VCL are lower than the stresses monitored at MPSB in the same backfilled stope condition. In addition, a peculiar phenomenon can be seen from Figure 9 (also in Figure 3-6) that there are some sudden increases near the stope bottom about the horizontal and vertical stresses of the backfill plotted along VCL. is is mainly caused by the restricted movement of the backfill layers near the stope bottom by the base rock to limit the mobilization and arching development of the backfill layers. However, this phenomenon has not happened for the monitored stresses at MPSB shown in Figure 9 because these stresses were collected at the fixed midpoint of the backfill bottom surface as shown in Figure 2a, which cannot illustrate the stress distributions along the height of the backfilled stope.
Besides, with the above two different stress plotting methods, both the vertical and horizontal stress distributions are not numerically stable at the interface normal stiffness of k n � 2.02 × 10 6 Pa/m, but the numerical results will no longer change when the normal stiffness k n reaches the values determined by the newly proposed method (equations (5) and (6)). It can be further explained that the stresses will not change with the variations of interface shear stiffness k s any more (not shown here) after using values by the proposed solution under the two plotting methods. ese comparisons further demonstrate the good reliability and adaptability of the proposed solution to determine the interface stiffness in different stress plotting methods within different mining processes with backfill.
However, in order to further reveal the different influences of the two plotting methods about backfill stresses to clarify their applicable conditions, the vertical (σ v ) and horizontal (σ h ) stresses monitored at MPSB of simulated stopes in the open stoping with subsequent backfill mining method and the upward cut and fill mining method have been shown in Figure 10. e analytical stresses with the Marston theory [12,13] by applying equation (7) and (8) have also added in Figure 10 for comparison. It can be seen from Figure 10, when using the stress monitoring method at MPSB, the horizontal stresses of the backfill in stopes of open stoping with subsequent backfill mining and upward cut and fill mining show almost the same distributions, which cannot present the significant differences in the convergent displacements by the closure deformation of stope rock walls, as shown in Figures 7 and8.
erefore, the effects of the mining-induced convergence of rock walls on backfill stresses can only be presented by plotting along VCL, whereas the stress monitoring at MPSB is not adaptable to rightly illustrate stress distributions of backfill in underground stopes after considering the mining closure, which is a very significant and common phenomenon in practical mining engineering [46][47][48][49].

Coefficient of Lateral Earth Pressure for Backfilled Stopes.
To further verify the reliability of the proposed solution for interface stiffness determination to simulate backfill stresses and to clearly distinguish the applicable conditions of two stress plotting methods, the coefficient of lateral earth pressure K earth under different interface normal stiffness k n in open stoping with subsequent backfill mining method has been shown in Figure 11. According to the definition in soil mechanics, the coefficient of lateral earth pressure K earth is a ratio of the horizontal stress over the vertical stress at different positions in soils, sands, fills, or other granular materials. It can be used to conveniently illustrate the stress state of filled materials. It is also an important input to do analytical calculations about backfill stresses [12][13][14][15][16][17] and strength requirements [3,7,8].
It can be seen from Figure 11(a) that with the increase of interface normal stiffness k n , the coefficient of lateral earth pressure K earth along with the whole height of backfilled stopes gradually decreases to numerically stable states and keeps constant after reaching the stiffness values by the newly proposed method, and the values of K earth generally are close to the at-rest state (K earth � K 0 ). However, Figure 11(b) shows that when the interface normal stiffness k n increases from around 2 MPa/m to 2000 MPa/m, the coefficients of lateral earth pressure K earth calculated with stresses by monitoring midpoint of stope bottom (MPSB) after filling each layer in stope are almost insensitive to the stiffness variations. erefore, if the investigations on the effects of variable interface stiffness values on backfill stresses in the two mining methods shown in Figure 3-6 were presented by the monitored stresses at MPSB, there would be no clear and right illustrations, which means the stress monitoring method at MPSB has the applicable limitations, and this is under the conditions of nonconsidering the convergence of rock walls acting on backfill in practical mine stopes. erefore, in the discussed conditions of considering different interface stiffness and the closure effects of stope walls, the stress plotting along the vertical central line (VCL) of the stope is a preferable method than monitoring the stresses at the midpoint of stope bottom (MPSB) after filling each layer in stopes.
It is worth further clarification that the presented earth pressure distributions with depth shown in Figure 11 were obtained in plane strain numerical models of mine stopes to simulate the backfill stresses under 2D arching effect. e variations of coefficient of lateral earth pressure in the context of 3D soil arching will be further investigated in the future by the authors and their collaborations.

Conclusions
To better understand stress distributions of backfill in mine stopes of two most widely used mining methods with backfill (i.e. open stoping with subsequent backfill, upward cut and fill mining), numerical simulations with FLAC3D were employed to do the calculations after inserting the necessary interface elements between backfill and rock walls. Analyses about the interface stiffness determination had been firstly carried out to propose a new reasonable solution to calculate the normal and shear stiffness, respectively. e effects of interface stiffness on stress distributions of backfilled stope under different excavating and filling processes in mining methods with backfill had been comprehensively investigated considering the accuracy and efficiency of the numerical simulations. Some key conclusions have been obtained as follows: (1) e proposed expressions to determine the apparent normal k n and shear k s stiffness of the interface elements between backfill and rock walls in simulated backfilled stopes can be suggested as k n � K/Δz min and k s � G/Δz min . e new validated solution will lead to totally the same backfill stress distributions with the results in numerical models applying the interface stiffness values by the empirical method in FLAC manual but will save at least 20%-30% calculation time to improve modeling efficiency under the same simulation conditions and will carry much clear physical meanings corresponding to the interaction between backfill and rock walls in mine stopes.
(2) With the proposed solution to determine interface stiffness in numerical models, the vertical and horizontal stresses along the vertical central line (VCL) of backfill in stopes of open stoping with subsequent backfill mining method show good agreements with the analytical stresses predicted by the Marston solution under the at-rest state (K earth � K 0 ). For the backfilled stopes in the upward cut and fill mining method, the vertical stress is much smaller and the horizontal stress is much larger than the analytical stresses by the Marston solution, which means a much stronger arching effect. e obvious differences of backfill stress distributions in mine stopes of two mining methods are mainly induced by different pressure of convergent displacements from rock walls.  stope bottom (MPSB) after filling each layer. Only at the top and the bottom surfaces of the backfilled stopes are the stress values the same for the two plotting methods. Both the vertical stress and the horizontal stress plotted along VCL are lower than the stresses monitored at MPSB in the same backfilled stope condition. (4) e effects of the mining-induced convergence of rock walls on backfill stress distributions and on the coefficient of lateral earth pressure K earth can only be presented by the plotted stresses along VCL, whereas the monitored stresses at MPSB is almost insensitive and not applicable to rightly illustrate stress distributions of backfill in underground stopes after considering the mining closure, which is a very significant and common phenomenon in practical mining engineering.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
ere are no conflicts of interest regarding the publication of this paper.