Application of Roadway Backﬁll Mining in Water-Conservation Coal Mining: A Case Study in Northern Shaanxi, China

: The mining induced subsidence and strata deformation are likely to a ﬀ ect the stability of the aquiclude, resulting in loss of water resources in the mining area. In order to reduce the disturbance of coal mining to the overlying strata and to preserve the water resources in the coal mining area, the roadway backﬁll mining (RBM) method was trialed in Yuyang coal mine in Northern Shaanxi, China. Based on pressure arch theory and ultimate strength theory, a mechanical model was developed to analyze the stability of coal pillars. Then the maximum number of vacant roadways between the mining face and the backﬁlling face was determined according to the stability of coal pillar and ﬁlling body. The method to calculate aquiclude subsidence and deformation was also proposed. Furthermore, as indicated by FLAC3D numerical simulations, the maximum tensile stress subjected by the aquiclude was 0.14 MPa, which is smaller than its tensile strength; the horizontal deformation was 0.24 mm / m, which is also smaller than the critical deformation of failure. Field monitoring data demonstrated a maximum of 2.76 m groundwater level drop in the mining area after mining. The groundwater level was determined to be 4.45~10.83 m below surface, ensuring the normal growth of surface vegetation and realizing the water-conservation coal mining (WCCM).


Introduction
Coal is the primary energy source in China. It is estimated that by 2050, coal will still account for 40% of China's energy consumption [1]. The underground mining method is used for about 92% of the coal resources extraction in China [2]. However, high-intensity underground coal mining has brought about a series of social and environmental problems [3,4]. The overlying strata of coal seam will collapse and subside after coal extraction. This is likely to damage the aquiclude of surface water and ground water, changing the recharge, runoff and discharge conditions and causing water resources pollution and loss. In regions with arid climate and lack of rainfall, the water loss caused by mining-induced damage to the aquiclude can bring about severe ecological consequences [5,6].
To achieve water-conservation coal mining (WCCM), it is necessary to reduce the mining-induced disturbance to the overlying strata, so that the mining-induced stress and deformation of the aquiclude are smaller than the tensile strength and critical deformation of failure [7]. At present, WCCM is achieved by partial mining and backfill mining [8,9]. The former includes strip mining, and room and pillar mining. In partial mining, only part of the coal is mined, while the remaining coal is left to support the roof so as to reduce the movement and deformation of the overlying strata (including Yuyang coal mine is located in North of Shaanxi, China, on the southern margin of Mu Us Desert, where the ground surface is dominated by sand and semi-fixed dunes (Figure 1) [29]. This region has an arid climate and lacks in precipitation, with an annual precipitation of only about 400 mm. Groundwater is vital for industrial and agricultural development as well as a healthy ecosystem [30]. The RBM scheme designed by Tiandi Science & Technology Co. Ltd. was implemented in Yuyang coal mine, where the gravity driven flow system was used as the backfilling approach. The paste backfilling material mainly consisted of aeolian sand. The surface aeolian sand in the mining area is mainly composed of quartz, feldspar and mica, with particle sizes vary between 0.075 mm and 0.600 mm (0.249 mm on average), which belongs to superfine sand. Therefore, the filling material could be easily transported using pipelines with little chance of blockage [31].
The No.3 coal seam, with an average thickness of 4 m, dip angle of 0.3~3° and burial depth of 100 m~120 m, was extracted in Yuyang coal mine. There is a loose sandy aquifer covering the surface. The underlying weathered strata, with an average thickness of 30 m, acts as the main aquiclude. The stratigraphic column of the mining area is shown in Figure 2.  Yuyang coal mine is located in North of Shaanxi, China, on the southern margin of Mu Us Desert, where the ground surface is dominated by sand and semi-fixed dunes (Figure 1) [29]. This region has an arid climate and lacks in precipitation, with an annual precipitation of only about 400 mm. Groundwater is vital for industrial and agricultural development as well as a healthy ecosystem [30]. The RBM scheme designed by Tiandi Science & Technology Co. Ltd. was implemented in Yuyang coal mine, where the gravity driven flow system was used as the backfilling approach. The paste backfilling material mainly consisted of aeolian sand. The surface aeolian sand in the mining area is mainly composed of quartz, feldspar and mica, with particle sizes vary between 0.075 mm and 0.600 mm (0.249 mm on average), which belongs to superfine sand. Therefore, the filling material could be easily transported using pipelines with little chance of blockage [31].
The No.3 coal seam, with an average thickness of 4 m, dip angle of 0.3~3° and burial depth of 100 m~120 m, was extracted in Yuyang coal mine. There is a loose sandy aquifer covering the surface. The underlying weathered strata, with an average thickness of 30 m, acts as the main aquiclude. The stratigraphic column of the mining area is shown in Figure 2.  In RBM, the transport roadway was driven along the strike in the middle of the mining panel. The mining roadways (MRs), with an average length of 85 m and width of 5.5 m, were arranged on each side of the transport roadway. The MRs were successively numbered, mined and backfilled in two separate stages. In the first stage, skip mining was performed for the odd-numbered MRs; and in the second stage, skip mining was performed for the remaining even-numbered MRs. Each MR was sealed with brick wall built at the MR ends and backfilled with paste backfilling material as soon as the coal was extracted, then the mining of the next roadway began. In this process, the backfilling face was always closely behind the mining face ( Figure 3). While during the actual underground production, the mining rate and backfilling rate were not always the same due to specific geological conditions of the coal seam and technical constraints, thus a number of VRs would be generated between the mining face and backfilling face ( Figure 4). The coal pillars between VRs would therefore bear the overlying load of VRs.
Protodyakonov's equilibrium arch (PEA) theory is a classical theory for the evaluation of roadways stability [32]. According to the PEA theory, during the excavation of a roadway, an arch-like zone is developed in the strata above the roadway [33]. Strata within this zone will subside and deform, while beyond this zone, the rocks remain stable and a self-bearing arch is formed. When there are several adjacent VRs, the arch structures above each roadway may join together forming a new arch, which can be referred as extended pressure arch structure ( Figure 5), and pillars between roadways will only subject to the overburden loads below the arch [34].
If there are few VRs, the extended pressure arch is small, as well as the load to coal pillars and subsidence. With the increase of VRs number, the arch extends further upward, and the load to coal pillars and subsidence also increase. When the arch reaches the aquiclude, the aquiclude will deform and subside, which may result in failure. In addition, the width of coal pillars between VRs are generally small (width-height ratio is less than 2) for the convenience of mining. If one of the coal pillars fails under heavy load, the load is then transferred to adjacent coal pillars, increasing the chance of chain damage of the coal pillar group. Therefore, it is necessary to determine the maximum number of VRs between mining face and backfilling face while ensuring that the coal pillars will not fail. In RBM, the transport roadway was driven along the strike in the middle of the mining panel. The mining roadways (MRs), with an average length of 85 m and width of 5.5 m, were arranged on each side of the transport roadway. The MRs were successively numbered, mined and backfilled in two separate stages. In the first stage, skip mining was performed for the odd-numbered MRs; and in the second stage, skip mining was performed for the remaining even-numbered MRs. Each MR was sealed with brick wall built at the MR ends and backfilled with paste backfilling material as soon as the coal was extracted, then the mining of the next roadway began. In this process, the backfilling face was always closely behind the mining face ( Figure 3). While during the actual underground production, the mining rate and backfilling rate were not always the same due to specific geological conditions of the coal seam and technical constraints, thus a number of VRs would be generated between the mining face and backfilling face ( Figure 4). The coal pillars between VRs would therefore bear the overlying load of VRs.
Protodyakonov's equilibrium arch (PEA) theory is a classical theory for the evaluation of roadways stability [32]. According to the PEA theory, during the excavation of a roadway, an arch-like zone is developed in the strata above the roadway [33]. Strata within this zone will subside and deform, while beyond this zone, the rocks remain stable and a self-bearing arch is formed. When there are several adjacent VRs, the arch structures above each roadway may join together forming a new arch, which can be referred as extended pressure arch structure ( Figure 5), and pillars between roadways will only subject to the overburden loads below the arch [34].
If there are few VRs, the extended pressure arch is small, as well as the load to coal pillars and subsidence. With the increase of VRs number, the arch extends further upward, and the load to coal pillars and subsidence also increase. When the arch reaches the aquiclude, the aquiclude will deform and subside, which may result in failure. In addition, the width of coal pillars between VRs are generally small (width-height ratio is less than 2) for the convenience of mining. If one of the coal pillars fails under heavy load, the load is then transferred to adjacent coal pillars, increasing the chance of chain damage of the coal pillar group. Therefore, it is necessary to determine the maximum number of VRs between mining face and backfilling face while ensuring that the coal pillars will not fail.  Figure 5 shows the influence of the VRs on overlying strata. For a stratum at a certain height, if there is only one VR (VR 1), a pressure arch is formed (arch a). As the stratum is above the top of arch a, no subsidence and deformation would occur on the stratum. If there are two VRs (VR 1 and VR 2), an extended pressure arch will be formed (arch b). Part of the stratum (CD section) falls into the range of arch b, then slight deformation and subsidence would occur, e.g., the straight line CD gradually becomes the arc CD. If there are three VRs (VR 1, VR 2 and VR 3), a larger pressure arch   Figure 5 shows the influence of the VRs on overlying strata. For a stratum at a certain height, if there is only one VR (VR 1), a pressure arch is formed (arch a). As the stratum is above the top of arch a, no subsidence and deformation would occur on the stratum. If there are two VRs (VR 1 and VR 2), an extended pressure arch will be formed (arch b). Part of the stratum (CD section) falls into the range of arch b, then slight deformation and subsidence would occur, e.g., the straight line CD gradually becomes the arc CD. If there are three VRs (VR 1, VR 2 and VR 3), a larger pressure arch will be formed (arch c). If the majority of the stratum falls in the range of arch c, significant deformation and subsidence would then occur, e.g., the straight line BE turns into the arc BE.
Sustainability 2017, 9, x FOR PEER REVIEW 5 of 23 will be formed (arch c). If the majority of the stratum falls in the range of arch c, significant deformation and subsidence would then occur, e.g., the straight line BE turns into the arc BE.

Coal Pillar Load
3.1.1. Protodyakonov's Equilibrium Arch Theory Figure 6 illustrates the general structure of the PEA. h is the mining height and b is the width of the MR (the dimension of MR is identical to VR and coal pillar between VRs, so h and b can also respectively represent the height and width of VR and coal pillar). H and 2W represent the height and width of the arch, respectively. The angle between the slip planes and the roadway ribs is 45°-φ/2, where φ is the internal friction angle of the overlying strata [33]. Based on the PEA theory, roadway roof pressure q can be calculated as follows: where γ is the average bulk density of the overlying strata within the pressure arch (N/m 3 ). The horizontal pressures e1 and e2 are evaluated as  will be formed (arch c). If the majority of the stratum falls in the range of arch c, significant deformation and subsidence would then occur, e.g., the straight line BE turns into the arc BE.

Coal Pillar Load
3.1.1. Protodyakonov's Equilibrium Arch Theory Figure 6 illustrates the general structure of the PEA. h is the mining height and b is the width of the MR (the dimension of MR is identical to VR and coal pillar between VRs, so h and b can also respectively represent the height and width of VR and coal pillar). H and 2W represent the height and width of the arch, respectively. The angle between the slip planes and the roadway ribs is 45°-φ/2, where φ is the internal friction angle of the overlying strata [33]. Based on the PEA theory, roadway roof pressure q can be calculated as follows: where γ is the average bulk density of the overlying strata within the pressure arch (N/m 3 ). The horizontal pressures e1 and e2 are evaluated as

Coal Pillar Load
3.1.1. Protodyakonov's Equilibrium Arch Theory Figure 6 illustrates the general structure of the PEA. h is the mining height and b is the width of the MR (the dimension of MR is identical to VR and coal pillar between VRs, so h and b can also respectively represent the height and width of VR and coal pillar). H and 2W represent the height and width of the arch, respectively. The angle between the slip planes and the roadway ribs is 45 • − ϕ/2, where ϕ is the internal friction angle of the overlying strata [33]. Based on the PEA theory, roadway roof pressure q can be calculated as follows: where γ is the average bulk density of the overlying strata within the pressure arch (N/m 3 ). The horizontal pressures e 1 and e 2 are evaluated as The span of the PEA can be calculated as follows: Given the symmetry of the PEA, only the left half of the arch is analyzed (Figure 7). Taken the arch top as the origin, a rectangular coordinate system is built. If the axial force at the origin is T, according to the PEA theory, the torque balance equation at any point M(x,y) on the arch is [34]: The value of the axial force T acting on the arch top equals to the horizontal thrust T at the arch feet, but they are opposite in direction. In order to ensure the structural stability of the PEA, the horizontal thrust T should meet the following condition: That is, the horizontal thrust T should be less than or equal to the friction generated by the vertical pressure of the overlying strata. The Protodyaknove's number f can be approximately calculated as one tenth of the uniaxial compressive strength of the surrounding rock or coal [33]. For safety reasons, T usually takes the value as half of the friction force, [34] i.e.,: According to Formulas (5) and (7), the arch line equation for the PEA can be written as: and the maximum height of the PEA is: Sustainability 2017, 9, x FOR PEER REVIEW 6 of 23 The span of the PEA can be calculated as follows: ( ) Given the symmetry of the PEA, only the left half of the arch is analyzed (Figure 7). Taken the arch top as the origin, a rectangular coordinate system is built. If the axial force at the origin is T, according to the PEA theory, the torque balance equation at any point M(x,y) on the arch is [34]: The value of the axial force T acting on the arch top equals to the horizontal thrust T′ at the arch feet, but they are opposite in direction. In order to ensure the structural stability of the PEA, the horizontal thrust T′ should meet the following condition: That is, the horizontal thrust T′ should be less than or equal to the friction generated by the vertical pressure of the overlying strata. The Protodyaknove's number f can be approximately calculated as one tenth of the uniaxial compressive strength of the surrounding rock or coal [33]. For safety reasons, T′ usually takes the value as half of the friction force, [34] i.e.:

= 2
Wfq T  (7) According to Formulas (5) and (7), the arch line equation for the PEA can be written as: and the maximum height of the PEA is:

Extended Pressure Arch Structure
As shown in Figure 8, if the number of VRs between the mining face and backfilling face is too large, a weak coal pillar, referred as pillar 0, may fail due to primary fissures in coal pillar or stress concentration induced by mining, then the two arches on both sides of the failed pillar would join together to form a new arch, which is referred as arch 2. Thus the coal pillar adjacent to the failed pillar, named pillar 1, will subject to two different horizontal forces from arch 1 and arch 2, and the resultant force is ( ) This resultant force will cause increased deformation of pillar 1, leading to the stress transfer from pillar 1 to pillar 2. Thus a bigger arch (i.e., arch 3) is formed. For the ease of expression, pillars without sufficient supporting force, like pillar 0 and pillar 1, are referred as "failed pillars". Then the horizontal resultant force on a pillar is: where m is the number of failed pillars. According to Formula (11), the more the failed pillars, the greater the horizontal resultant force. A greater horizontal resultant force will cause the failure of the next pillar, resulting in continued damage of the coal pillar group until the horizontal force is balanced by the horizontal thrust from unmined coal seam [34,35]. This will eventually result in the extended pressure arch structure, as shown in Figure 9 [36,37].

Extended Pressure Arch Structure
As shown in Figure 8, if the number of VRs between the mining face and backfilling face is too large, a weak coal pillar, referred as pillar 0, may fail due to primary fissures in coal pillar or stress concentration induced by mining, then the two arches on both sides of the failed pillar would join together to form a new arch, which is referred as arch 2. Thus the coal pillar adjacent to the failed pillar, named pillar 1, will subject to two different horizontal forces from arch 1 and arch 2, and the resultant force is This resultant force will cause increased deformation of pillar 1, leading to the stress transfer from pillar 1 to pillar 2. Thus a bigger arch (i.e., arch 3) is formed. For the ease of expression, pillars without sufficient supporting force, like pillar 0 and pillar 1, are referred as "failed pillars". Then the horizontal resultant force on a pillar is: where m is the number of failed pillars. According to Formula (11), the more the failed pillars, the greater the horizontal resultant force. A greater horizontal resultant force will cause the failure of the next pillar, resulting in continued damage of the coal pillar group until the horizontal force is balanced by the horizontal thrust from unmined coal seam [34,35]. This will eventually result in the extended pressure arch structure, as shown in Figure 9 [36,37].

Extended Pressure Arch Structure
As shown in Figure 8, if the number of VRs between the mining face and backfilling face is too large, a weak coal pillar, referred as pillar 0, may fail due to primary fissures in coal pillar or stress concentration induced by mining, then the two arches on both sides of the failed pillar would join together to form a new arch, which is referred as arch 2. Thus the coal pillar adjacent to the failed pillar, named pillar 1, will subject to two different horizontal forces from arch 1 and arch 2, and the resultant force is ( ) This resultant force will cause increased deformation of pillar 1, leading to the stress transfer from pillar 1 to pillar 2. Thus a bigger arch (i.e., arch 3) is formed. For the ease of expression, pillars without sufficient supporting force, like pillar 0 and pillar 1, are referred as "failed pillars". Then the horizontal resultant force on a pillar is: where m is the number of failed pillars. According to Formula (11), the more the failed pillars, the greater the horizontal resultant force. A greater horizontal resultant force will cause the failure of the next pillar, resulting in continued damage of the coal pillar group until the horizontal force is balanced by the horizontal thrust from unmined coal seam [34,35]. This will eventually result in the extended pressure arch structure, as shown in Figure 9 [36,37]. According to Formulas (4), (8) and (9), the span and height of the extended pressure arch are respectively: where N is the number of VRs. As known from the effective area theory [38], the support area of the k-th coal pillar is Ak ( Figure  9). Therefore, the load subjected by the k-th coal pillar on the left (right) side of the Y axis is: The pillar bearing the maximum load is the first coal pillar on the left (right) side of the Y axis, i.e.:

Coal Pillar Strength
The coal pillar strength depends on the shape and size of the coal pillar, and the strength of coal sample as well. The commonly used empirical formula for coal pillar strength determination is given in Table 1.  According to Formulas (4), (8) and (9), the span and height of the extended pressure arch are respectively: where N is the number of VRs. As known from the effective area theory [38], the support area of the k-th coal pillar is A k (Figure 9). Therefore, the load subjected by the k-th coal pillar on the left (right) side of the Y axis is: The pillar bearing the maximum load is the first coal pillar on the left (right) side of the Y axis, i.e.,:

Coal Pillar Strength
The coal pillar strength depends on the shape and size of the coal pillar, and the strength of coal sample as well. The commonly used empirical formula for coal pillar strength determination is given in Table 1.  According to the ultimate strength theory, when the load subjected by the coal pillar exceeds its compressive strength, the coal pillar will fail. Since coal pillar with b/h < 2 is more likely to fail, a safety factor is introduced to ensure the stability of the coal pillar [39,40].
where F is the safety factor (1.5~2); P max is the maximum load borne by the coal pillar (MPa).

Maximum Number of VRs
Combining Formulas (16) and (17), together with the pillar strength σ p obtained through empirical formulas in Table 1, the critical safety distance W e between mining face and backfilling face can be determined. Substituting W e into Formula (12), then the integral part of the solution N will be the maximum number of VRs which could ensure the stability of coal pillars.

Uncertainties of the Analysis
In Yuyang coal mine, the mechanical properties of rocks were obtained from laboratory tests with standard samples (Φ50 mm × 100 mm). Mining parameters such as width of the MR and mining height were analyzed as influencing factors of the maximum number of VRs by orthogonal experiment. The mining height is set as 3~4.5 m. In order to improve mining efficiency and reduce the roof subsidence before backfilling, the width of the MR is set as 5~6 m [12]. Each factor was applied to several levels, i.e., 3 m, 3.5 m, 4 m and 4.5 m for the mining height, and 5 m, 5.5 m and 6 m for the MR width. All the values of parameters are shown in Table 2.     It can be observed from Figure 10 that the maximum load subjected by coal pillar showed an approximately linear increase with the increase of VRs number, while the pillar strength is irrelevant to the number of VRs. When the number of VRs reaches a certain value, the maximum load will be greater than the pillar strength indicating the possibility of failure of coal pillar. With the increase of MR width, the strength of pillar and the maximum load also increase. It is noted that the former (i.e., pillar strength increase) is conducive to maintain the stability of coal pillars while the latter (i.e., load increase) has negative influence on pillars stability. As the increasing rate of the  It can be observed from Figure 10 that the maximum load subjected by coal pillar showed an approximately linear increase with the increase of VRs number, while the pillar strength is irrelevant to the number of VRs. When the number of VRs reaches a certain value, the maximum load will be greater than the pillar strength indicating the possibility of failure of coal pillar. With the increase of MR width, the strength of pillar and the maximum load also increase. It is noted that the former (i.e., pillar strength increase) is conducive to maintain the stability of coal pillars while the latter (i.e., load increase) has negative influence on pillars stability. As the increasing rate of the maximum pillar load is higher than that of pillar strength, the maximum numbers of VRs decrease with the increase of MR width. Mining height has a significance influence on pillar strength while its influence on pillar load is minor.
With the increase of mining height, the pillar strength decreases, which will cause the decrease of VRs numbers. In summary, although using large size of MR could reduce the times of movements for mining face and backfilling face and improve mining efficiency, its influence on stability of coal pillar is negative. Therefore, the parameters need to be optimized further according to the actual situation.
It can be concluded from Table 3 that, with the width of the MR vary from 5 m to 6 m and mining height vary from 3 m to 4.5 m, the maximum number of VRs in the first mining stage is 8~10, while in the second stage it is 3~4. Considering mining safety, mining efficiency and backfilling effects, the width of MR is determined as 5.5 m and the maximum number of VRs is determined as 8 in the first stage and 3 in the second stage.

Immediate Roof Subsidence
During the two stages of RBM, the roof subsidence U mainly consists of three parts, the roof subsidence (U 1 ) after mining in the first stage, the roof subsidence (U 2 ) after mining in the second stage, and the roof subsidence (U 3 ) after final stabilization of overlying strata, i.e.,: After the extraction of roadways in the first stage, the load above the roadway is transferred to the coal pillars on each side of the roadway. The immediate roof subsides under extra loading, and the subsidence is U 1 . So the height of the filling body in the first stage is h-U 1 . During solidification of the paste backfill material which is mainly composed of aeolian sand, the size of the filling body reduces slightly. With small roof subsidence after backfilling, the support provided by the filling body is insufficient and the load is mainly borne by the coal pillar.
During the second stage of mining, the load previously borne by the coal pillars is transferred to the filling body as the recovery of coal pillars, leading to compression of the filling body and the subsidence of immediate roof with an amount of U 2 . Thus the height of the filling body in the second stage is h-U 1 -U 2 .
As the expansion of mining area, the extent of mining-induced overlying strata subsidence will eventually develop upwards to the ground surface. At this stage, the load of the coal seam is the weight of all strata ranging from the roof to the ground surface, and the immediate roof subsidence is U 3 .

Roof Subsidence U 1
After mining in the first stage, the two sides of the coal pillars are subjected to bidirectional compression. Under high load, the coal pillars tend to yield characterized with low load bearing capacity and considerable horizontal deformation. The middle section of the coal pillars is subjected to tri-directional compression characterized with high load bearing capacity and small horizontal deformation, which is beneficial to maintain its supporting capacity. From the sides to the inside of the coal pillars, the plastic zone and elastic zone are successively developed. The overall load bearing pattern is an isosceles trapezoid [41].
As shown in Figure 11, the load acting on elastic zone of coal pillar is: Due to the short existence time of the coal pillars, the creep characteristics are ignored. Then the width of the plastic zone bp can be calculated as follows [42]: (20) where λ is the lateral pressure coefficient which can be calculated using Poisson's ratio ν (i.e., λ=ν/(1-ν)); C is the cohesion (Pa). According to Formulas (16), (19) and (20), the load σ can be calculated, and the maximum roof subsidence after mining in the first stage is given by: where Ec is the Young's modulus of the coal pillar (Pa).

Roof Subsidence U2
During mining in the second stage, the overlying roof is supported by the filling body constructed in the first stage. According to Formulas (16), (19) and (20), the load σ′ can be calculated using mechanical parameters of the filling material. Thus the maximum roof subsidence after mining in the second stage can be calculated as: where Ef is the Young's modulus of the filling body (Pa); σ′ is the load borne by filling body in the second stage (Pa).

Roof Subsidence U3
After the stabilization of overlying strata movement, the weight of all the overlying strata will be acted on the backfilling bodies constructed during the two mining stages. According to the condition of static equilibrium, the following formula can be obtained: where γ′ is the average bulk density of all the overlying strata (N/m 3 ); H′ is the burial depth of coal seam (m). According to Formula (23), the roof subsidence U3 can be determined by: Uniaxial compressive tests of the standard samples (Φ50 mm×100 mm) of the coal and filling body were carried out in laboratory. Then the slope of the axial stress-axial strain curves in elastic deformation stage of the standard samples was estimated for the Young's modulus and the slope of Due to the short existence time of the coal pillars, the creep characteristics are ignored. Then the width of the plastic zone b p can be calculated as follows [42]: where λ is the lateral pressure coefficient which can be calculated using Poisson's ratio ν (i.e., λ = ν/(1 − ν)); C is the cohesion (Pa). According to Formulas (16), (19) and (20), the load σ can be calculated, and the maximum roof subsidence after mining in the first stage is given by: where E c is the Young's modulus of the coal pillar (Pa).

Roof Subsidence U 2
During mining in the second stage, the overlying roof is supported by the filling body constructed in the first stage. According to Formulas (16), (19) and (20), the load σ can be calculated using mechanical parameters of the filling material. Thus the maximum roof subsidence after mining in the second stage can be calculated as: where E f is the Young's modulus of the filling body (Pa); σ is the load borne by filling body in the second stage (Pa).

Roof Subsidence U 3
After the stabilization of overlying strata movement, the weight of all the overlying strata will be acted on the backfilling bodies constructed during the two mining stages. According to the condition of static equilibrium, the following formula can be obtained: where γ is the average bulk density of all the overlying strata (N/m 3 ); H is the burial depth of coal seam (m). According to Formula (23), the roof subsidence U 3 can be determined by: Uniaxial compressive tests of the standard samples (Φ50 mm × 100 mm) of the coal and filling body were carried out in laboratory. Then the slope of the axial stress-axial strain curves in elastic deformation stage of the standard samples was estimated for the Young's modulus and the slope of the lateral strain-axial strain curves was estimated for the Poisson's ratio [29,31], and the corresponding values are E c = 3.26 GPa, ν c = 0.25, E f = 0.13 GPa, ν f = 0.29, respectively. According to Formulas (18), (21), (22) and (24), the maximum roof subsidence U 1 , U 2 , U 3 and U can be calculated ( Figure 12). U 1 and U 2 are caused by the weight of rocks within the arch in the first and second stage, while U 3 is caused by the weight of strata ranging from arch top to the ground surface. So with the increase of VRs number, the height of the arch increases while the height of strata from arch top to the ground surface decreases, leading to the increase of U 1 and U 2 and the decrease of U 3 . Regarding U, with the increase of U 1 and U 2 , the height of filling bodies in the first and second stage (i.e., h-U 1 and h-U 1 -U 2 .) will decrease, resulting in the decrease of U.
Sustainability 2017, 9, x FOR PEER REVIEW 13 of 23 the lateral strain-axial strain curves was estimated for the Poisson's ratio [29,31], and the corresponding values are Ec=3.26 GPa, νc=0.25, Ef=0.13 GPa, νf=0.29, respectively. According to Formulas (18), (21), (22) and (24), the maximum roof subsidence U1, U2, U3 and U can be calculated ( Figure 12). U1 and U2 are caused by the weight of rocks within the arch in the first and second stage, while U3 is caused by the weight of strata ranging from arch top to the ground surface. So with the increase of VRs number, the height of the arch increases while the height of strata from arch top to the ground surface decreases, leading to the increase of U1 and U2 and the decrease of U3. Regarding U, with the increase of U1 and U2, the height of filling bodies in the first and second stage (i.e., h-U1 and h-U1-U2.) will decrease, resulting in the decrease of U.

Subsidence of Aquiclude and Surface
A planar mechanical model as shown in Figure 13 is developed to investigate the overlying strata subsidence in RBM. Assuming the subsidence area of the i-th stratum on this plane to be Si (from Figure 2，i≤10) where n is the total number of MRs; Hi is the distance from the upper surface of the i-th stratum to the upper surface of the coal seam (m); α is the influence angle of overlying strata subsidence (°); wi(x) is the subsidence of the i-th stratum (m). After backfilling, bending deformation with a small deflection occurs to the strata above the main roof (i≥2). It is reasonable to treat the subsidence curve as straight lines [43], therefore, the subsidence area Si of the i-th stratum can be expressed as the area of an inverted isosceles trapezoid, i.e., where β is the angle of full subsidence (°). Assuming the maximum subsidence of the main roof is equal to that of the immediate roof, i.e., w2(0)=U, then S2 can be written as: It is assumed that each stratum above the main roof has the same subsidence area, i.e., Si=S2 (i≥2), then the maximum subsidence at the top of the i-th stratum is expressed as:

Subsidence of Aquiclude and Surface
A planar mechanical model as shown in Figure 13 is developed to investigate the overlying strata subsidence in RBM. Assuming the subsidence area of the i-th stratum on this plane to be S i (from Figure 2, i ≤ 10), then (25) where n is the total number of MRs; H i is the distance from the upper surface of the i-th stratum to the upper surface of the coal seam (m); α is the influence angle of overlying strata subsidence ( • ); w i (x) is the subsidence of the i-th stratum (m). After backfilling, bending deformation with a small deflection occurs to the strata above the main roof (i ≥ 2). It is reasonable to treat the subsidence curve as straight lines [43], therefore, the subsidence area S i of the i-th stratum can be expressed as the area of an inverted isosceles trapezoid, i.e., where β is the angle of full subsidence ( • ). Assuming the maximum subsidence of the main roof is equal to that of the immediate roof, i.e., w 2 (0) = U, then S 2 can be written as: It is assumed that each stratum above the main roof has the same subsidence area, i.e., S i = S 2 (i ≥ 2), then the maximum subsidence at the top of the i-th stratum is expressed as: The equation for the subsidence curve at the top of the i-th stratum is written as: When calculating the subsidence of each stratum by Formula (28), the load bearing capacity of each stratum is not considered, so verification is carried out on the calculation. Assuming each overburden stratum is Euler-Bernoulli beam with unit width, then the overburden load Q i and length l i of the i-th stratum are respectively: where γ j is the bulk density of the j-th stratum (kN/m 3 ); h j is the thickness of the j-th stratum (m). If the stiffness and thickness of the i-th stratum are large enough to bear the overlying load, bed separation may occur in its underneath. In this case, the maximum subsidence of the i-th stratum can be calculated as the deflection of a clamped-clamped beam subjected to uniform load, i.e., where E i is the Young's modulus of the i-th stratum (GPa); h i is the thickness of the i-th stratum (m). If w i (0) ≤ w i (0), the i-th stratum subsides with the underlying stratum together and the subsidence is w i (0); If w i (0) > w i (0), the i-th stratum is separated from the underlying stratum and the subsidence of the i-th stratum is w i (0). The geological condition of Yuyang coal mine is simple. There is little change in the thickness and burial depth of coal seam. Considering the small research area which was only limited to the RBM panel, the thickness and Young's modulus were assumed to be constant. The thickness of each stratum was determined according to the drilling data. As for the Young's modulus, the uniaxial compression tests were performed on five rock samples from each stratum, then the slopes of the axial stress-axial strain curves in elastic deformation stage were averaged to obtain the Young's modulus of the stratum. The parameters used for the calculation of overlying strata subsidence are shown in Table 4. From Table 4, it can be observed that w i (0) < w i (0) for each stratum, indicating the subsidence of the i-th stratum is w i (0). The study showed that all the overlying strata subside together and there was no separation between two adjacent strata using RBM. The maximum subsidence of the aquiclude and surface are 23 mm and 22 mm respectively. 15  From Table 4, it can be observed that wi(0)<wi ) for each stratum, indicating the subsidence of the i-th stratum is wi(0). The study showed that all the overlying strata subside together and there was no separation between two adjacent strata using RBM. The maximum subsidence of the aquiclude and surface are 23 mm and 22 mm respectively. Figure 13. The planar mechanical model for the overlying strata subsidence and deformation.

Numerical Model
FLAC 3D (Fast Lagrangian Analysis of Continua in 3 Dimensions), an explicit finite-difference program developed by Itasca Consulting Group Inc., is one of the most widely used numerical simulation software in mining and geo-engineering [22]. The mohr-coulomb model is a classical constitutive model used to describe yield and failure behaviors of soil and rock [24]. In this study, FLAC 3D simulations were conducted to investigate the stress and deformation of the aquiclude using mohr-coulomb model. The transport roadway was developed along the x-axis of the numerical model, and a total of 100 MRs were designed on the two sides of the transport roadway. Considering the boundary effect, the model was developed with an extra 200 m wide boundary coal pillar surrounding the mining area as Figure 14b shows. The model size was 680 m×585 m×190 m, and the finite difference meshes generated by FLAC 3D was shown in Figure 14a. Roller constraint was assigned to the four sides of the model, fixed constraint were used at the bottom of the model, and no constraint was applied to the top. The mechanical properties of rocks and backfilling materials used for the simulation were shown in Table 5 [29,31].

Numerical Model
FLAC 3D (Fast Lagrangian Analysis of Continua in 3 Dimensions), an explicit finite-difference program developed by Itasca Consulting Group Inc., is one of the most widely used numerical simulation software in mining and geo-engineering [22]. The mohr-coulomb model is a classical constitutive model used to describe yield and failure behaviors of soil and rock [24]. In this study, FLAC 3D simulations were conducted to investigate the stress and deformation of the aquiclude using mohr-coulomb model. The transport roadway was developed along the x-axis of the numerical model, and a total of 100 MRs were designed on the two sides of the transport roadway. Considering the boundary effect, the model was developed with an extra 200 m wide boundary coal pillar surrounding the mining area as Figure 14b shows. The model size was 680 m × 585 m × 190 m, and the finite difference meshes generated by FLAC 3D was shown in Figure 14a. Roller constraint was assigned to the four sides of the model, fixed constraint were used at the bottom of the model, and no constraint was applied to the top. The mechanical properties of rocks and backfilling materials used for the simulation were shown in Table 5 [29,31].

Modeling Scheme
As can be seen from Section 4.2, the subsidence of overlying strata increases with the increase of VRs number, so the numerical simulation was carried out with 8 VRs in the first stage and 3 VRs in the second stage which was taken as the worst condition. During the numerical modeling, mining was simulated by applying null model to the MR which was Mohr-Coulomb model, backfilling of a VR was simulated by applying Mohr-Coulomb model to the VR which was null model.

Aquiclude Stability
The subsidence curve of aquiclude and surface obtained from numerical simulation are compared with that determined by Formula (29) (Figure 15). Tensile stress and deformation of the aquiclude was calculated to evaluate the stability of aquiclude (Figures 16-19). The correlation coefficient between two subsidence curves of the aquiclude is 0.94. And for the two surface subsidence curves, the correlation coefficient is 0.96. The maximum subsidence of aquiclude calculated by Formula (29) is 23 mm, which is 11.5% smaller than the numerical simulation results (26 mm); the maximum subsidence of surface calculated by formula (29) is 22 mm, which is 4.3% smaller than the numerical simulation results (23 mm). It demonstrates the good agreement between the theoretical model and the numerical model.          After backfilling, the maximum tensile stress of the aquiclude was calculated to be 0.14 MPa, which was smaller than its tensile strength; the maximum horizontal stretching deformation was 0.24 mm/m, which was smaller than the critical deformation for failure (Table 6). Therefore, it can be concluded from the model results that the aquiclude will remain stable during the mining process, through which both the ground water and surface water resources could be effectively preserved [29,44].

Analysis of An Engineering Application
In Yuyang mine, all of the MRs were backfilled after mining, and brick walls with a thickness of 700 mm were built at both ends of the MRs to prevent filling material from flowing out of the MRs. To ensure the stability of coal pillar, the number of VRs between mining face and backfilling face was kept at 4 or lower. In the first stage, a number of 120 MRs were mined and backfilled. One After backfilling, the maximum tensile stress of the aquiclude was calculated to be 0.14 MPa, which was smaller than its tensile strength; the maximum horizontal stretching deformation was 0.24 mm/m, which was smaller than the critical deformation for failure (Table 6). Therefore, it can be concluded from the model results that the aquiclude will remain stable during the mining process, through which both the ground water and surface water resources could be effectively preserved [29,44].

Analysis of An Engineering Application
In Yuyang mine, all of the MRs were backfilled after mining, and brick walls with a thickness of 700 mm were built at both ends of the MRs to prevent filling material from flowing out of the MRs. To ensure the stability of coal pillar, the number of VRs between mining face and backfilling face was kept at 4 or lower. In the first stage, a number of 120 MRs were mined and backfilled. One month after the first stage was over, field tests were conducted onsite to ensure the strength of filling body and its roof supporting rate meet the requirements, mining and backfilling in the second stage could then begin. A total of 278 000 tons of coal had been extracted using RBM, and more than 212 000 tons of filling material had been used to control the subsidence of overlying strata.
Both subsidence and water level were monitored during the application of RBM at Yuyang mine to evaluate its impact on water conservation. The maximum subsidence above the backfill mining area was measured to be 21 mm, which was significantly lower than that of fully mechanized longwall caving faces (86 mm) of the mine [45]. It is noted that, two boreholes (ZK0905 and ZK0906 as shown in Figure 20) were originally employed for water level observation before and after mining. However, due to some construction and maintenance issues, water level data after mining were not available at these two boreholes. Therefore, three more boreholes (YYT1, YYT2 and YYT3 in Figure 20) were drilled after mining to investigate the water levels. Table 7 shows the water levels monitored at each observation hole. Considering the relative positions of the boreholes, water altitudes at the area near ZK906 and YYT1 were compared, and a decrease of 2.76 m was found; water altitudes near ZK905 and YYT2 were also compared showing a 0.45 m drop. as shown in Figure 20) were originally employed for water level observation before and after mining. However, due to some construction and maintenance issues, water level data after mining were not available at these two boreholes. Therefore, three more boreholes (YYT1, YYT2 and YYT3 in Figure 20) were drilled after mining to investigate the water levels. Table 7 shows the water levels monitored at each observation hole. Considering the relative positions of the boreholes, water altitudes at the area near ZK906 and YYT1 were compared, and a decrease of 2.76 m was found; water altitudes near ZK905 and YYT2 were also compared showing a 0.45 m drop.  When the groundwater level is within 8~15 m below the surface, the surface vegetation can grow normally [46]. After coal mining in the studied area, the water level was 4.45~10.83 m (7.33 m on average) below surface, and it was less than 8 m for most of the region; therefore, the normal growth of surface vegetation was ensured. Field survey also demonstrated the generally good surface vegetation in the mining area, indicating the achievement of WCCM.

Discussion
Compared with the longwall backfill mining, it is more flexible to arrange the working face of RBM. The mining system and backfilling system is able to operate independently in different roadways, which not only ensures efficiency, but also avoids mutual interference. The exposed roof area is smaller after coal extraction from the roadway and the two sides of the roadways are always supported by the coal rib or the filling body. Therefore, the roof deformation and subsidence are small before backfilling, which is beneficial to the effective control of overlying strata and aquiclude. Considering these advantages, RBM has great application potential in WCCM.  When the groundwater level is within 8~15 m below the surface, the surface vegetation can grow normally [46]. After coal mining in the studied area, the water level was 4.45~10.83 m (7.33 m on average) below surface, and it was less than 8 m for most of the region; therefore, the normal growth of surface vegetation was ensured. Field survey also demonstrated the generally good surface vegetation in the mining area, indicating the achievement of WCCM.

Discussion
Compared with the longwall backfill mining, it is more flexible to arrange the working face of RBM. The mining system and backfilling system is able to operate independently in different roadways, which not only ensures efficiency, but also avoids mutual interference. The exposed roof area is smaller after coal extraction from the roadway and the two sides of the roadways are always supported by the coal rib or the filling body. Therefore, the roof deformation and subsidence are small before backfilling, which is beneficial to the effective control of overlying strata and aquiclude. Considering these advantages, RBM has great application potential in WCCM.
In order to reduce the subsidence of the roadway roof and facilitate mining, the width of roadway or coal pillar is generally less than 6 m, which can be categorized into slender pillar (b/h < 2). The slender pillar is likely to fail under additional loads. Once failed, it will transfer pressure to adjacent pillars, resulting in the formation of an extended pressure arch structure. Therefore, with the increase of pillar width, the load subjected by coal pillar increases, the maximum number of VRs thus decreases. However, when the pillar width reaches a certain degree, there is an original stress zone in the middle of the wide coal pillar, then the pressure arch on the two sides of the wide coal pillar will not be able to form a new extended pressure arch. In this case, there is no limit on the number of VRs. If the coal pillar does not need to be recovered in the second stage, and the MRs can be mined without backfilling in the first stage, the RBM then becomes the traditional strip mining.
The paste backfill material is mainly composed of aeolian sand, and it takes time for the filling body to reach its designed strength. Within a short time after backfilling, the roof is still not effectively supported. So the number of roadways that have been backfilled but without effective roof support should be calculated according to mining rate, backfilling rate, and the properties change of filling body over time. These roadways should also be considered as VRs when determining the maximum number of VRs.
The subsidence curves of aquiclude and surface calculated by theoretical analysis were compared with that predicted by numerical simulation, and the correlation coefficients were 0.96 and 0.94 respectively, indicating good agreement of the two methods. Surface subsidence in RBM area was studied by theoretical analysis, numerical simulation and field measurement, and the results were 22 mm, 23 mm and 28 mm respectively. The results obtained from theoretical analysis and numerical simulation were similar, both of which were smaller than the field measurement. The reason for this difference could be attributed to the slightly reduced volume of filling material during its solidification, resulting in decreased roof supporting rate.
RBM was trialed at the studied mining panel. However, the adjacent mining panel was mined without backfilling and the roof strata caved into the gob after coal extraction. Although the overlying aquiclude above the studied mining panel is effectively preserved, the groundwater can still flow towards adjacent zones where backfilling was not used. Therefore, water may run off laterally along the aquifer or even leak further downwards. Consequently, the decline of groundwater level after mining cannot realistically reflect the WCCM effect of RBM. If RBM were applied to the entire mining area, the WCCM effect would be further enhanced.

Conclusions
Based on pressure arch theory and ultimate strength theory, a mechanical model was developed to analyze the stability of coal pillar group. The influence of mining height and MR width on the maximum number of VRs between mining face and backfilling face was analyzed. With a mining height of 3~4.5 m and MR width of 5.5 m, the maximum number of VRs in the two stages were 8 and 3, respectively. Then the stability of coal pillar could be ensured.
Based on the pressure arch theory, the stratum subsidence equation was derived, which was then verified by the subsidence equation developed by clamped beam theory considering the bearing capacity of each stratum. It was found that the stratum subsidence decreases with the decrease of stratum depth. The maximum subsidence of aquiclude and surface were calculated to be 23 mm and 22 mm, respectively.
The theoretical model was in good agreement with the numerical model as indicated by the high correlation coefficients (0.96 for the aquiclude subsidence curves and 0.94 for surface subsidence curves). Through numerical simulation, the maximum horizontal tensile stress and horizontal deformation of the aquiclude was 0.14 MPa and 0.24 mm/m respectively, which were lower than the tensile strength and the critical deformation for horizontal tensile failure, indicating the aquiclude would remain intact after mining.
After RBM, the surface subsidence of the mining panel decreased by 21 mm. The groundwater level in the mining area dropped by a maximum of 2.76 m, and the ground water level was 4.45~10.83 m below surface, which is sufficient to ensure the normal growth of surface vegetation. Therefore, it can be concluded that WCCM is achieved at Yuyang Mine by RBM.
Author Contributions: Y.Y. conceived the research and wrote the paper. L.M. revised and reviewed the manuscript. All authors have read and approved the final manuscript.