A Vertical Joint Spacing Calculation Method for UDEC Modeling of Large-Scale Strata and Its Inﬂuence on Mining-Induced Surface Subsidence

: While universal discrete element code (UDEC) is widely used for understanding the mechanism of large-scale strata movement and the effects of mining subsidence on the environment, the fundamental knowledge of how to set vertical joint spacing (VJS) in UDEC is still not fully understood. To address the knowledge gap, we ﬁrst present a novel VJS calculation method, then conduct UDEC experiments, and ﬁnally compare the predictions of UDEC models with ﬁeld subsidence observation. The results suggest the following: (1) when compared to the conventional VJS setting (1 × to 3 × bed thickness), the maximum surface subsidence (MSS) prediction via UDEC models based on our proposed VJS setting method is closer to ﬁeld observation; (2) a smaller but varying VJS setting can also have the effect of a larger VJS setting; and (3) with the increase in VJS, MSS ﬁrst drops, then rises, and reaches the minimum when VJS is set at approximately 7 × bed thickness. This paper provides an explanation of the VJS setting in UDEC and establishes a bridge between the KS theory and VJS, which is helpful for the sustainable development of such an UDEC modeling strategy and for a better understanding of the inﬂuences of mining subsidence on the environment in mining-affected areas.


Introduction
Underground coal mining causes the deformation and movement of overlying strata and eventually leads to ground surface subsidence or cracking, which affects our living environment and the sustainable development of coal industries [1,2]. The discrete element commercial software UDEC is extensively used for understanding the mechanism behind strata movement and to predict surface subsidence more directly and realistically. However, the knowledge of how to establish a sound UDEC model especially involving discontinuity parameters is not entirely understood and merits further study.
In the UDEC modeling process of large-scale strata, it is difficult to construct detailed discontinuities such as joints and cracks of different directions and scales, mainly due to the following two reasons: (1) under the constraints of UDEC software algorithms and computer computing power, simulating a large number of joints will lead to extremely long calculation times; (2) under the current technology level of geological exploration, the detailed spatial distribution of discontinuities cannot be given. In such a case, the modeling idea of equivalent joint (i.e., simplified or fictitious joint) is popular [3][4][5][6]. Specifically, researchers often use a horizontal joint to simulate the bedding plane and use evenly spaced vertical joints (ESVJs) that are perpendicular to the bedding plane and run through the entire rock bed to simulate the phenomenon of stratum fracture. It should be mentioned that, in the field, although vertical joints (also called cross joints) [7,8] may actually exist in some strata in local, they are not always perpendicular to bedding planes and are not always run through the entire rock bed [9][10][11]. More commonly, they have limited length and intersect bedding planes at different angles in terms of their sedimentary environment and geological history. In addition, the spacing between vertical joints varies greatly from one place to another; some researchers [12][13][14][15] suggest that the spacing is proportional to bed thickness and the ratio of joint spacing to bed thickness may range from 0.77 to 100 as reported in literatures [16,17].
To be clear, the ESVJs are used to comprehensively characterize the cross joints of different directions and scales in the field; it is not necessarily a real existence but a simulation strategy in the context of large-scale modeling to achieve a balance between calculation time and simulation effects. Therefore, vertical joint spacing (VJS) will play a decisive role in determining the length of the fractured rock mass. If the VJS is too large, the actual broken stratum cannot be reflected in simulations in time; if the VJS setting is too small, the simulated stratum may break prematurely. Apparently, the setting of VJS determines the shape of the overburden structure to a certain extent and ultimately affects the degree of surface subsidence, whereas such effects have not been understood.
In practical, it is very common to use ESVJs to simulate large-scale strata movement. For example, Zhang et al. [18] built a UDEC model with dimensions of 800 m × 280 m for revealing the overburden failure characteristics in the case of a special mining and geological condition. Li et al. [19] constructed a UDEC model with dimensions of 850 m × 100 m for analyzing the influence of key stratum on mine pressure in the case of top-coal caving face with super great mining height. Zhang et al. [20] built a UDEC model with dimensions of 2200 m × 550 m for evaluating the fracture development around river valleys. The value of the VJS in the above and many other studies [18][19][20][21][22][23][24][25][26][27][28][29][30] is between 1× to 3× bed thickness, but how the VJS was determined is still not clear and there is no related literature that clearly explains the determination process. To the best of the authors' knowledge, practical experiences and subjective tendencies often play a leading role, which makes it impossible to compare between different UDEC models involving different engineering cases and limits the sustainable development of such UDEC modeling strategy.
In this study, we first proposed a method for calculating the VJS on the basis of the Key Stratum theory and the "Three-Zones" theory, which provides a direct explanation on how to set up the VJS in UDEC when modeling large-scale strata movement and therefore is beneficial for the continuous progress of such modeling strategy. The proposed method is then applied to a specific UDEC model based on the Ying-Pan-Hao coal mine, while other UDEC models with the VJS varying from 0.5× to 1000× bed thickness were also built (see Section 3) as a control group and for fundamental understandings of how the VJS influences surface subsidence prediction. Finally, the effectiveness of the proposed method was evaluated by comparing the difference between field observations and predictions in Section 4.

General Description
With the ongoing extraction of underground coal resources, the stress field of overburden will also be dynamically redistributed. During this process, the overlying rock bed failure starts at the weakest part where stress concentration exceeds its strength and then propagates in the least resistance path or, in other words, along the discontinuities (either a predominate one or a set of mircojoints). Recall that, in the field, the discontinuities between bedding planes have different lengths, directions, and persistence, and that ESVJs as a kind of fictitious joints were used to comprehensively characterize these discontinuities for large-scale modeling. Then, it seems that the major function of the vertical joints (VJs) are to mark the potential fracture locations of the overlying strata so as to simulate the fracture phenomenon of rock beds because, if not, they lose their physical meaning. More specifically, if they are not marking the potential fraction locations, they must indicate the real joints or need another theory to explain the meaning of setting up such fictious joints. Since there is no related theory in this regard for now and the real existence of ESVJs in the field is very rare as we have stated in the Introduction, we believe that assuming the vertical joint spacing (VJS) is essentially the same as the length of broken block in the Key Stratum theory [31] should be reasonable. In the following, we elaborate the method for calculating VJS by adding an additional hypothesis to the KS theory, and the effect of "Three-Zones" theory on VJS setting is also discussed.
Moreover, we must stress that the above analyses do provide an explanation for the setting of ESVJs. Although it may not be perfect, it is helpful for the sustainable development of such UDEC modeling strategy.

Determining VJS Based on the Key Stratum Theory
The Key Stratum (KS) theory is proposed by Qian [31] and is widely accepted by researchers, which emphasizes that hard rock layers play a key role in the control of movement and the deformation of overburden. While KS theory is a big leap for understanding the possible overburden structure after coal mining, the load acting on the stratum between two adjacent key strata is not specified, which results in only the VJS of key stratum being available. Rather, the VJS of the remaining rock layers is still not clear. Here, we present a direct calculation method for solving this knowledge gap as elaborated below.
As shown in Figure 1, we presume that there are M layers strata above a coal seam, and denote their thickness, volume force, elastic modulus, and tensile strength as h i , γ i , E i , and T i , respectively, where 1 i M. According to the composite beam used in KS theory, the load (q n ) 1 borne by the first layer of an n layers composite beam is given as: The first step of our proposed method is consistent with the step of identifying potential key stratum in KS theory, i.e., looking for the hard rock layers in the overburden and specifically to determine if Equation (2) is true. If so, the rock layer n + 1 can be considered as a potential key stratum, and then by making the layer n + 1 the new first layer, we can reuse Equation (2) to find all potential key strata: Supposing that we have obtained N 1 potential key strata, and one of them is numbered n i with the corresponding load q n i , where n i represents the layer n i and 1 i N 1 . Equation (3) can then be used for calculating the breaking distance l n i (i.e., VJS), where Equation (3) is derived under the condition of a fixed beam at both ends: In KS theory, the potential key strata also need to meet the condition that the breaking distance of the lower one is smaller than that of the upper one (see Equation (4)). Otherwise, all the rock loads controlled by the upper one need to be applied to the lower one, and then the breaking distance of the lower one is recalculated until all key strata are  In KS theory, the potential key strata also need to meet the condition that the breaking distance of the lower one is smaller than that of the upper one (see Equation (4)). Otherwise, all the rock loads controlled by the upper one need to be applied to the lower one, and then the breaking distance of the lower one is recalculated until all key strata are determined, i.e., Equation (4) holds for all potential key strata. Apparently, KS theory still does not describe the calculation method of the breaking distance of the rock layers between two adjacent key strata. Here, we add an additional hypothesis to solve this problem. Supposing that we have obtained N 2 key strata, and one of them is numbered k i with the corresponding load q k i , where k i represents the layer k i of the overburden and 1 i N 2 . The number of the rock layers controlled by layer k i can therefore denote as N 3 = k i+1 − k i − 1. When the layer k i is broken, the N 3 + 1 layers composite beam with layer k i as its first layer can then be regarded as an N 3 layers composite beam with layer k i + 1 as the first layer, as shown in Figure 2. In KS theory, the potential key strata also need to meet the condition that the breaking distance of the lower one is smaller than that of the upper one (see Equation (4)). Otherwise, all the rock loads controlled by the upper one need to be applied to the lower one, and then the breaking distance of the lower one is recalculated until all key strata are determined, i.e., Equation (4) holds for all potential key strata. Apparently, KS theory still does not describe the calculation method of the breaking distance of the rock layers between two adjacent key strata. Here, we add an additional hypothesis to solve this problem. Supposing that we have obtained N2 key strata, and one of them is numbered with the corresponding load , where represents the layer of the overburden and 1 ≪ ≪ . The number of the rock layers controlled by layer can therefore denote as = − − 1. When the layer is broken, the + 1 layers composite beam with layer as its first layer can then be regarded as an N3 layers composite beam with layer + 1 as the first layer, as shown in Figure 2. Hence, the load acting on layer + 1 is given as: By analogy, we can obtain the load acting on the layer + ( < ) as: Collectively, the load of every rock layer can be determined, and by introducing it into Equation (3), their VJS can be directly and clearly calculated. Hence, the load acting on layer k i + 1 is given as: By analogy, we can obtain the load acting on the layer k i + j (j < N 3 ) as: Collectively, the load of every rock layer can be determined, and by introducing it into Equation (3), their VJS can be directly and clearly calculated.

VJS Setting Based on the "Three-Zones" Theory
The "Three-Zones" theory [32] is also widely accepted by researchers, which divides the overburden into three zones: caved zone, fracture zone, and continuous deformation zone. Basically, the caved zone is a space where roof rocks collapsed and piled together after coal mining; the fracture zone is a space where strata have been separated or fractured but still can form a certain balanced structure without collapse; and the continuous deformation zone is a space where only bending deformation occurred. This suggests that it is not necessary to set up vertical joints in the continuous deformation zone since there is no stratum fracture, which is helpful to reduce modeling complexity and calculation time, especially for deep mining, where large amounts of joints are involved. However, the impact of such a setting on the prediction of surface subsidence remains unclear. In this paper, we first set up vertical joints in full modeling space to determine the characteristic of three zones, then build a same model with vertical joints in continuous deformation zone removed and re-execute the simulation process to observe the impacts on surface subsidence.

Modeling Prototype and Field Subsidence Observation
Our UDEC experiments are based on the Ying-Pan-Hao (YPH) coal mine located in the Wu-Shen banner of China's Inner Mongolia province. As shown in Figure 3, there are several boreholes around the 2201 panel, of which K82 and K83 are obviously the most representative, and there are 163 observation stations along and vertical the mining direction for surface subsidence measurement. The 2201 panel is the first mining panel of the YPH coal mine with no buildings or structures on the ground surface, which ensures the reliability of surface subsidence measurements. Its mining length, working face length, average coal thickness, and average depth are 2500 m, 300 m, 6.5 m, and 730 m, respectively. We list the strata information revealed by boreholes K82 and K83 in Tables 1 and 2, and more details of the mining and geological conditions of the YPH coal mine can be found in Gong and Guo [33]. Maximum surface subsidence was recorded to 411 mm after mining the 2201 panel, and additional comments will be made when discussing the differences between field subsidence observation and results predicted by UDEC models (Section 4).
it is not necessary to set up vertical joints in the continuous deformation zone since there is no stratum fracture, which is helpful to reduce modeling complexity and calculation time, especially for deep mining, where large amounts of joints are involved. However, the impact of such a setting on the prediction of surface subsidence remains unclear. In this paper, we first set up vertical joints in full modeling space to determine the characteristic of three zones, then build a same model with vertical joints in continuous deformation zone removed and re-execute the simulation process to observe the impacts on surface subsidence.

Modeling Prototype and Field Subsidence Observation
Our UDEC experiments are based on the Ying-Pan-Hao (YPH) coal mine located in the Wu-Shen banner of China's Inner Mongolia province. As shown in Figure 3, there are several boreholes around the 2201 panel, of which K82 and K83 are obviously the most representative, and there are 163 observation stations along and vertical the mining direction for surface subsidence measurement. The 2201 panel is the first mining panel of the YPH coal mine with no buildings or structures on the ground surface, which ensures the reliability of surface subsidence measurements. Its mining length, working face length, average coal thickness, and average depth are 2500 m, 300 m, 6.5 m, and 730 m, respectively. We list the strata information revealed by boreholes K82 and K83 in Tables 1 and  2, and more details of the mining and geological conditions of the YPH coal mine can be found in Gong and Guo [33]. Maximum surface subsidence was recorded to 411 mm after mining the 2201 panel, and additional comments will be made when discussing the differences between field subsidence observation and results predicted by UDEC models (Section 4).

UDEC Models and Setup
We built a total of 48 UDEC models, half of which are K82-based models (scenario 1), and the other half are K83-based models (scenario 2). We take scenario 1 as an example to explain the experimental process and setup. According to Table 1, we first built 11 UDEC models with the VJS set at 0.5×, 1×, 2×, 3×, 4×, 5×, 6×, 7×, 8×, 9×, and 1000× bed thickness, respectively, in order to study the influence of the VJS on surface subsidence. We also built a special UDEC model with the VJS calculated via the equations in Section 2.1 for comparison with the aforementioned models. We then conducted UDEC simulations to obtain the ranges of the "Three-Zones" and built another 12 UDEC models with vertical joints removed in the continuous deformation zone (e.g., Figure 4) to study the influence of the "Three-Zones" theory on surface subsidence. The VJS in coal seams is constant to ensure that the excavation range in all models is exactly the same. All models have the same dimension of 1902.8 m × 724.4 m, the same boundary condition with left, right, and bottom sides fixed (see Figure 5), and the same constitutive model (Mohr-Coulomb model for rock mass and Coulomb slip model for joints). For evenly spaced vertical joints (ESVJs), the position of the first vertical joint from the left boundary determines the position of the following vertical joints. In all models, we set ESVJs with the left boundary as the starting position to ensure that there is no manual manipulation of the joint position. The rock mass properties used here were calibrated in the previous study [33] and are shown in The rock mass properties used here were calibrated in the previous study [33] and are shown in Figures 6 and 7.

Influence of VJS on Maximum Surface Subsidence
Based on the methodology and engineering background elaborated above, we simulate the mining process of the 2201 panel, record the maximum surface subsidence (MSS) Figure 6. Calibrated rock mass properties used for strata of borehole K82. Strata numbers are consistent with that in Table 1.   Table 2.

Influence of VJS on Maximum Surface Subsidence
Based on the methodology and engineering background elaborated above, we simu-  Table 2.

Influence of VJS on Maximum Surface Subsidence
Based on the methodology and engineering background elaborated above, we simulate the mining process of the 2201 panel, record the maximum surface subsidence (MSS) of different scenarios (K83-based models and K82-based models), and present them in Figure 8. We found that, in general, the MSS predicted by the UDEC models based on borehole K83 is larger than that based on borehole K82, which is mainly because that the thickness and lithology of the stratum were revealed by different boreholes, and the mechanical properties used in different models vary. This suggests that a 3DEC (3D Discrete Element Code) model is needed to incorporate all strata data into one model so as to improve the accuracy of predictions, which, as a consequence, will also lead to an extremely complex modeling process and long calculation time and merits further study.
We also found that, with the increase ofin the VJS, the MSS first drops (stage 1), then rises (stage 2), and reaches the minimum at 7× bed thickness for both scenarios, which is significant because it is different from the intuition that MSS will continue to decrease as the VJS increasingincreases. Here, we explain this phenomenon as illustrated in Figure 9, where we can see that the length of the cantilever beam [34] in Figure 9b is much longer than that in Figure 9a, which well reduces the space for overburden to further overhang and collapse. We also note that the length of the cantilever beams in Figure 9d,e are almost the same, but due to the different VJS settings, the patterns of masonry beam [29,35,36] formed in Layer 31 vary. Specifically, the denser vertical joints in Figure 9d lead to greater subsidence of Layer 31, which in turn provides the overlying strata more space for deformation. Collectively, if we consider this from the perspective of a compressive arch [34], then we see that, with the increase ofin VJS setting, the earlier formation of compressive arch structure reduces deformation space, which explains the mechanism of stage 1. For stage 2, we consider an extreme case where the VJS is equal to 1000 times bed thickness (see Figure 9c,f). Then, we imagine that Figure 9c,f are transformed from Figure 9b,e, respectively. Or inIn other words, the rocks that have collapsed and fractured return to their original positions with bending and separation, and the only difference before and after this transformation is that the compressive arch structure in Figure 9c,f bear more rock loads that have been released in Figure 9b,e, which results in the greater subsidence and explains the mechanism of stage 2. We found that, in general, the MSS predicted by the UDEC models based on borehole K83 is larger than that based on borehole K82, which is mainly because the thickness and lithology of the stratum were revealed by different boreholes, and the mechanical properties used in different models vary. This suggests that a 3DEC (3D Discrete Element Code) model is needed to incorporate all strata data into one model so as to improve the accuracy of predictions, which, as a consequence, will also lead to an extremely complex modeling process and long calculation time and merits further study.
We also found that, with the increase in the VJS, the MSS first drops (stage 1), then rises (stage 2), and reaches the minimum at 7× bed thickness for both scenarios, which is significant because it is different from the intuition that MSS will continue to decrease as the VJS increases. Here, we explain this phenomenon as illustrated in Figure 9, where we can see that the length of the cantilever beam [34] in Figure 9b is much longer than that in Figure 9a, which well reduces the space for overburden to further overhang and collapse. We also note that the length of the cantilever beams in Figure 9d,e are almost the same, but due to the different VJS settings, the patterns of masonry beam [29,35,36] formed in Layer 31 vary. Specifically, the denser vertical joints in Figure 9d lead to greater subsidence of Layer 31, which in turn provides the overlying strata more space for deformation. Collectively, if we consider this from the perspective of a compressive arch [34], then we see that, with the increase in VJS setting, the earlier formation of compressive arch structure reduces deformation space, which explains the mechanism of stage 1. For stage 2, we consider an extreme case where the VJS is equal to 1000 times bed thickness (see Figure 9c,f). Then, we imagine that Figure 9c,f are transformed from Figure 9b,e, respectively. In other words, the rocks that have collapsed and fractured return to their original positions with bending and separation, and the only difference before and after this transformation is that the compressive arch structure in Figure 9c,f bear more rock loads that have been released in Figure 9b,e, which results in the greater subsidence and explains the mechanism of stage 2.  Tables 1 and 2. More importantly, the UDEC models corresponding to Figure 9c,f are very similar to FLAC2D or FLAC3D models where no vertical joints are set in abundant literatures [37][38][39][40][41][42][43][44][45], which suggests that when predicting the surface subsidence using FLAC2D or FLAC3D, the surface subsidence could be overestimated due to the fact that the collapse of roof strata above a coal seam cannot be effectively simulated [46].

Influence of the Proposed Method on MSS
Regarding the influence of our proposed VJS calculation method on MSS, we can see ( Figure 10) that, in K83-based models, the MSS (714 mm) predicted by our proposed model is between the values predicted by models with VJS setting at 3× bed thickness and 4× bed thickness, while in K82-based models, that (567 mm) is between the values predicted by models with the VJS setting at 4× and 5× bed thickness. If we consider that in abundant literatures the VJS is set between 1× and 3× bed thickness and the observed MSS in the YPH coal mine is 411 mm, then we can conclude that the prediction by our proposed model is better because it is closer to the field observation. More importantly, we can see (Figure 11) that the VJS based on the proposed method varies in different strata and is no more than 3× bed thickness in most thick beds, suggesting that a smaller but varying VJS setting can also have the effect of a larger VJS setting in large-scale UDEC simulations. We also note that K82-based and K83-based models with the VJS set at 4× to 9× bed thickness achieved smaller MSS predictions, which are closer to field observation, but such a large VJS setting has no precedent and runs counter to the breaking distance of rock roof in the YPH coal mine.  Tables 1 and 2. More importantly, the UDEC models corresponding to Figure 9c,f are very similar to FLAC2D or FLAC3D models where no vertical joints are set in abundant literatures [37][38][39][40][41][42][43][44][45], which suggests that when predicting the surface subsidence using FLAC2D or FLAC3D, the surface subsidence could be overestimated due to the fact that the collapse of roof strata above a coal seam cannot be effectively simulated [46].

Influence of the Proposed Method on MSS
Regarding the influence of our proposed VJS calculation method on MSS, we can see (Figure 10) that, in K83-based models, the MSS (714 mm) predicted by our proposed model is between the values predicted by models with VJS setting at 3× bed thickness and 4× bed thickness, while in K82-based models, that (567 mm) is between the values predicted by models with the VJS setting at 4× and 5× bed thickness. If we consider that in abundant literatures the VJS is set between 1× and 3× bed thickness and the observed MSS in the YPH coal mine is 411 mm, then we can conclude that the prediction by our proposed model is better because it is closer to the field observation. More importantly, we can see (Figure 11) that the VJS based on the proposed method varies in different strata and is no more than 3× bed thickness in most thick beds, suggesting that a smaller but varying VJS setting can also have the effect of a larger VJS setting in large-scale UDEC simulations. We also note that K82-based and K83-based models with the VJS set at 4× to 9× bed thickness achieved smaller MSS predictions, which are closer to field observation, but such a large VJS setting has no precedent and runs counter to the breaking distance of rock roof in the YPH coal mine. Sustainability 2021, 13, x FOR PEER REVIEW 11 of 15   Tables 1 and 2. In addition, it should be mentioned that the observed MSS in the YPH coal mine was recorded right after mining the 2201 panel, and further surface subsidence will occur afterward and account for 15% of the total subsidence according to He et al. [47]. Therefore, the final MSS may reach 483 mm, which suggests that the K82-based UDEC model is more representative for the 2201 panel. We note that the subsidence values (567 mm and 714 mm) obtained via the UDEC models still have a certain gap with the field observation (483 mm), suggesting that mining subsidence is a complex phenomenon that is unlikely to be accurately explained by UDEC models based on a single borehole and that a data-intensive UDEC model that considers many other factors, such as modeling based on multiple boreholes, is needed and merits further study.

Influence of "Three-Zones" Theory on MSS
Regarding the influence of the "Three-Zones" theory on MSS, we can see (Figure 12) that in K83-based models considering the "Three-Zones" theory, the MSS is slightly smaller when the VJS is set at 0.5× to 1× bed thickness, is larger when the VJS is set at 2× to 9× bed thickness, and is the same when the VJS is set at 1000× bed thickness. Similar relationships can be found in K82-based predictions. For models based on our proposed method, a larger MSS prediction can also be observed. Since the VJS is set larger than 1×  In addition, it should be mentioned that the observed MSS in the YPH coal mine was recorded right after mining the 2201 panel, and further surface subsidence will occur afterward and account for 15% of the total subsidence according to He et al. [47]. Therefore, the final MSS may reach 483 mm, which suggests that the K82-based UDEC model is more representative for the 2201 panel. We note that the subsidence values (567 mm and 714 mm) obtained via the UDEC models still have a certain gap with the field observation (483 mm), suggesting that mining subsidence is a complex phenomenon that is unlikely to be accurately explained by UDEC models based on a single borehole and that a data-intensive UDEC model that considers many other factors, such as modeling based on multiple boreholes, is needed and merits further study.

Influence of "Three-Zones" Theory on MSS
Regarding the influence of the "Three-Zones" theory on MSS, we can see (Figure 12) that in K83-based models considering the "Three-Zones" theory, the MSS is slightly smaller when the VJS is set at 0.5× to 1× bed thickness, is larger when the VJS is set at 2× to 9× bed thickness, and is the same when the VJS is set at 1000× bed thickness. Similar relationships can be found in K82-based predictions. For models based on our proposed method, a larger MSS prediction can also be observed. Since the VJS is set larger than 1× In addition, it should be mentioned that the observed MSS in the YPH coal mine was recorded right after mining the 2201 panel, and further surface subsidence will occur afterward and account for 15% of the total subsidence according to He et al. [47]. Therefore, the final MSS may reach 483 mm, which suggests that the K82-based UDEC model is more representative for the 2201 panel. We note that the subsidence values (567 mm and 714 mm) obtained via the UDEC models still have a certain gap with the field observation (483 mm), suggesting that mining subsidence is a complex phenomenon that is unlikely to be accurately explained by UDEC models based on a single borehole and that a dataintensive UDEC model that considers many other factors, such as modeling based on multiple boreholes, is needed and merits further study.

Influence of "Three-Zones" Theory on MSS
Regarding the influence of the "Three-Zones" theory on MSS, we can see (Figure 12) that in K83-based models considering the "Three-Zones" theory, the MSS is slightly smaller when the VJS is set at 0.5× to 1× bed thickness, is larger when the VJS is set at 2× to 9× bed thickness, and is the same when the VJS is set at 1000× bed thickness. Similar relationships can be found in K82-based predictions. For models based on our proposed method, a larger MSS prediction can also be observed. Since the VJS is set larger than 1× bed thickness in abundant literatures, we can conclude that removing vertical joints in the continuous deformation zone will result in the MSS increasing. We also note that, with the increase in the VJS, the MSS difference before and after removing vertical joints in the continuous deformation zone basically first increases and then decreases in both scenarios. We further calculate this difference and compare it with the error between measured and predicted MSS to get the percentage of its influence on the total error, as shown in Figure 13. We can see that, when the VJS is set to less than 1× bed thickness, the influence of the "Three-Zones" theory is considered to be very small, less than 5%, which we believe is acceptable because removing vertical joints in a continuous deformation zone will greatly reduce joints number, thereby reducing the complexity of modeling and saving calculation time, especially considering that a small VJS setting will produce a large number of joints in the case of deep coal mining. We also note that, for our proposed VJS setting method, the influence of the "Three-Zones" theory can reach 12.1% to 15.5%, which we believe is not acceptable because when the VJS is large (e.g., more than 3× bed thickness), the influence may be great, as shown in Figure 13, especially considering that the VJS in each stratum often differs.

Summary and Conclusions
Recall that our objectives of this study were to gain knowledge of a clear and direct VJS calculation method for UDEC modeling of large-scale strata movement, and to better We also note that, with the increase in the VJS, the MSS difference before and after removing vertical joints in the continuous deformation zone basically first increases and then decreases in both scenarios. We further calculate this difference and compare it with the error between measured and predicted MSS to get the percentage of its influence on the total error, as shown in Figure 13. We can see that, when the VJS is set to less than 1× bed thickness, the influence of the "Three-Zones" theory is considered to be very small, less than 5%, which we believe is acceptable because removing vertical joints in a continuous deformation zone will greatly reduce joints number, thereby reducing the complexity of modeling and saving calculation time, especially considering that a small VJS setting will produce a large number of joints in the case of deep coal mining. We also note that, for our proposed VJS setting method, the influence of the "Three-Zones" theory can reach 12.1% to 15.5%, which we believe is not acceptable because when the VJS is large (e.g., more than 3× bed thickness), the influence may be great, as shown in Figure 13, especially considering that the VJS in each stratum often differs. We also note that, with the increase in the VJS, the MSS difference before and after removing vertical joints in the continuous deformation zone basically first increases and then decreases in both scenarios. We further calculate this difference and compare it with the error between measured and predicted MSS to get the percentage of its influence on the total error, as shown in Figure 13. We can see that, when the VJS is set to less than 1× bed thickness, the influence of the "Three-Zones" theory is considered to be very small, less than 5%, which we believe is acceptable because removing vertical joints in a continuous deformation zone will greatly reduce joints number, thereby reducing the complexity of modeling and saving calculation time, especially considering that a small VJS setting will produce a large number of joints in the case of deep coal mining. We also note that, for our proposed VJS setting method, the influence of the "Three-Zones" theory can reach 12.1% to 15.5%, which we believe is not acceptable because when the VJS is large (e.g., more than 3× bed thickness), the influence may be great, as shown in Figure 13, especially considering that the VJS in each stratum often differs.

Summary and Conclusions
Recall that our objectives of this study were to gain knowledge of a clear and direct VJS calculation method for UDEC modeling of large-scale strata movement, and to better

Summary and Conclusions
Recall that our objectives of this study were to gain knowledge of a clear and direct VJS calculation method for UDEC modeling of large-scale strata movement, and to better understand the influence of VJS on surface subsidence. Specifically, we proposed a novel VJS calculation method and conducted UDEC experiments. We find that with the increase in the VJS, the MSS first drops, then rises, and reaches the minimum approximately at 7× bed thickness. We find that when the VJS setting is extremely large, UDEC models are similar to FLAC models where no vertical joints are set. We thus suggest that surface subsidence predicted by FLAC2D or FLAC3D models could be overestimated. We also find that when compared to conventional VJS setting (1× to 3× bed thickness), the MSS predictions via UDEC models based on our proposed VJS setting method is closer to measurement and that a smaller but varying VJS setting can also have the effect of a larger VJS setting.
Regarding the influence of "Three-Zones" theory on MSS, we find that removing vertical joints in a continuous deformation zone will result in MSS increasing and that with the increase in the VJS, the MSS difference before and after removing vertical joints basically first increases and then decreases. We suggest not removing vertical joints in the continuous deformation zone when building UDEC models based on our proposed VJS calculation method and suggest removing them when the VJS setting is small (e.g., 0.5× to 1× bed thickness in this case).
Given the importance and extensive application of UDEC and also in the context of precision coal mining [48], a clear and direct VJS calculation method and a better understanding of influences of the VJS on surface subsidence, as discussed in this study, are necessary and important, which helps for the sustainable development of the UDEC modeling strategy of large-scale strata and for a better understanding of the influences of mining subsidence on the environment in mining-affected areas.