Abstract

The caved zone during longwall mining has high permeability, resulting in a mass of groundwater storage which causes a threat of groundwater inrush hazard to the safe mining. To investigate the hazard mechanism of granular sandstone and mudstone mixture (SMM) in caved zone, this paper presents an experimental study on the effect of sandstone particle (SP) and mudstone particle (MP) weight ratio on the non-Darcy hydraulic properties evolution. A self-designed granular rock seepage experimental equipment has been applied to conduct the experiments. The variation of particle size distribution was induced by loading and water seepage during the test, which indicated that the particle crushing and erosion properties of mudstone were higher than those of sandstone. Porosity evolution of SMM was strongly influenced by loading (sample height) and SP/MP weight ratio. The sample with higher sample height and higher weight ratio of SP achieved higher porosity value. In particular, a non-Darcy equation, for hydraulic properties (permeability and non-Darcy coefficient ζ) calculation, was sufficient to fit the relation between the hydraulic gradient and seepage velocity. The test results indicated that, due to the absence and narrowing of fracture and void during loading, the permeability κ decreases and the non-Darcy coefficient ζ increases. The variation of the hydraulic properties of the sample within the same particle size and SP/MP weight ratio indicated that groundwater inrush hazard showed a higher probability of occurrence in sandstone strata and crushed zone (e.g., faults). Moreover, isolated fractures and voids were able to achieve the changeover from self-extension to interconnection at the last loading stage, which caused the fluctuation tendency of κ and ζ. Fluctuation ability in mudstone was higher than that in sandstone. The performance of an empirical model was also investigated for the non-Darcy hydraulic properties evolution prediction of crushing and seepage processes. The predictive results indicated that particle crushing and water erosion caused the increase of hydraulic properties, being the main reason that the experimental values are typically higher than those obtained from the predictive model. The empirical model has a high degree of predictive accuracy; however, κ has a higher predictive accuracy than ζ. Furthermore, the predictive accuracy of κ increases and ζ decreases with increasing weight ratio of SP.

1. Introduction

As shown in Figure 1, Xiaojihan Coal Mine is located in Yulin City, Northern Shaanxi Province, China. The interbedded deposit of sandstone and mudstone in Xiaojihan Coal Mine is mainly formed in the Lower Cretaceous and the Upper Jurassic periods above the coal seam [1]. The interbedded deposit formed in Lower Cretaceous period, that is, Luohe formation, is about 60–90 m in thickness. The interbedded deposit formed in Upper Jurassic period includes the Anding formation, the Zhiluo formation, and the Yan’an formation. The thickness of the three formations is from 90 m to 200 m, 200 m to 320 m, and 320 m to 560 m, respectively.

Longwall mining is an efficient production technology of coal mining with several advantages, including high recovery and safe working environment; therefore, it is adopted in Xiaojihan Coal Mine. In longwall mining, as shown in Figure 2, there are four disturbance zones that are formed in the overlying strata of a coal seam [2]. Equilibrium will be progressed in the overburden strata, with the consolidation of the caved zone (gob) after mining extraction [3]. The granular rocks are composed of an immediate roof and overlying strata, which were filling in the caved zone termed as “gob.” The caved zone (gob) has high porosity and permeability resulting in a mass of water storage, which is primarily derived from weak and fractured adjacent coal beds and formations [4]. The vertical thickness of a gob (caved zone) with a highly irregular shape and different granular rocks can reach 4–11 times the mining height [5]. Therefore, the strata of sandstone and mudstone in Xiaojihan Coal Mine can be naturally filled in the caved zone of mixture state.

Recent laboratory and field investigations show that groundwater inrush hazard is mostly localized to the gob [6]. In Xiaojihan Coal Mine, the coal mining activity under the Luohe formation aquifer is threatened by the rich source of groundwater, which is responsible for the increase in the hydraulic conductivity of gob and roof strata [7]. Therefore, large amount volumes of groundwater will inrush into exposed space, for example, working face or tunnel [810]. Considering the effect of some more factors such as the hydraulic pressure [11], geological setting, and extraction percentage in underground mines [1214], the hydraulic properties of the granular rocks in gob (caved zones) have drawn attention [15].

Significant permeability variation of the granular medium under the effect of hydraulic pressure with compacted regions and the formation of pathways has been reported [16, 17]. The relationship between the pressure gradient and the velocity of the groundwater inrush in granular rock mass does not generally satisfy the Darcy equation, but a distinct nonlinear correlation [15, 18], that is, non-Darcy flow Forchheimer equation [19]. It is of great theoretical and practical significance to establish a nonlinear flow model for the identification of seepage mechanisms and the reasonable prediction of groundwater inrush [20]. A nonlinear dynamic model of variable mass system for non-Darcy flow in granular rocks was established based on the mechanical plug in porous media [6]. Hou et al. 2018 [21] applied a non-Darcy flow model (Forchheimer equation) to describe the water flow in karst collapse (a kind of granular rock); the dynamic process [22] of mining-induced water inrush mechanism in karst collapse has also been investigated by numerical simulation. Therefore, it is also a key issue for experimental study of hydraulic properties of granular sandstone-mudstone mixtures SMM samples, especially the non-Darcy hydraulic properties prediction using the limited laboratory data.

This paper aims to investigate the groundwater inrush mechanism of mining-induced granular SMM in cave zone. To carry out the study, self-designed granular rock seepage experimental equipment has been invented, which is able to quantify the effect of sandstone and mudstone weight ratio and sample height (porosity) on the hydraulic properties of granular SMM. An empirical model is also applied to analyze the non-Darcy hydraulic properties evolution for the crushing and seepage processes.

2. Experimental Methods

2.1. Granular Rock Materials

The sandstone and mudstone blocks, which were in Zhiluo and Yan’an formation formed in Upper Jurassic period and excavated from drills in Xiaojihan Coal Mine, were used in this study. The dry densities of the sandstone and mudstone, which were determined in the laboratory, are 2.561 × 103 kg/m3 and 2.423 × 103 kg/m3, respectively.

In this test, a vertical vibrating-shaker with a pan was used for dry-sieving. There was a series of selected standard sieves with decreasing size installed in the experimental system. By setting nearly 50 Hz of the frequency and loading rate and using 40% of the relative upper sieve volume lasting for 10 min, the fractions were prepared by sieving with different sizes, as 20 mm, 15 mm, 12 mm, 10 mm, 8 mm, 5 mm, and 2.5 mm, respectively. Firstly, to prepare the tested materials, rock blocks were hammered into particles with variable sizes less than 20 mm diameters. Secondly, the granular rock particles were then separated into seven size ranges using sieves: <2.5 mm, 2.5–5 mm, 5–8 mm, 8–10 mm, 10–12 mm, 12–15 mm, and 15–20 mm. Thirdly, particle sizes above were divided into groups according to their quantity; SMM samples were prepared by mixing particles according to the particle size distribution (PSD) curve shown in Figure 3.

Tested granular rock materials, granular sandstone-mudstone mixtures (SMM), were mixed together according to the weight ratio of sandstone particle (SP) mixture to mudstone particle (MP) mixture 1 : 0, 2 : 1, 1 : 1, 1 : 2, and 0 : 1, respectively. The weight of each granular sample was 1800 g. The PSD curves of the five tested granular rock materials were consistent with the curve shown in Figure 3. The weight of SP and MP of five samples is, respectively, mixed according to the PSD listed in Table 1.

2.2. Experimental Systems

As shown in Figure 4, the experimental system for testing water seepage in granular rocks comprises four different parts: (a) sample height control equipment, (b) hydraulic pressure control equipment, (c) automatic data collection equipment, and (d) granular rock seepage equipment.

Sample height control equipment shown in Figure 4(a) is made up of an oil pump, a hydraulic cylinder, a relief valve, and a linear variable differential transformer (LVDT). To achieve different height levels, a steady-loading pressure is conducted using the hydraulic cylinder and the oil pump, and then the LVDT is applied to measure the sample height.

A water pump and a relief valve are the main components of the hydraulic pressure control equipment (Figure 4(b)), which supplies a stable and constant hydraulic pressure to the granular rock sample. The water pump is conducted to fill the cylindrical tube with water by a constant flow pressure.

Automatic data collection equipment (Figure 4(c)) includes a hydraulic pressure sensor, a flow rate transducer, a data recorder, and a personal computer, which can collect the data of the hydraulic pressure and flow velocity during tests.

Granular rock seepage equipment (Figure 4(d)) is the key component of this experimental system, which is mainly made up of a cylindrical tube, a loading platen, two porous plates and a pedestal. Figure 5(a) shows the connection of the designed equipment. The loading plate and loading platen are conducted to offer the steady-loading pressure. The porous plate (Figure 5(b)) can ensure a steady water seepage. Two epoxy resin rings are applied to seal the aperture between the cylindrical tube and porous plates, guaranteeing the axial flow of water. The axial gap between the pedestal and the cylinder tube is enclosed by O-shaped rubber seal rings. This self-designed testing system for water flow in granular rock can supply flow quantity for a long time so that it can support the steady state flow test.

2.3. Testing Methods and Procedure

In the following part, detailed experimental steps for investigation on hydraulic properties of the SMM are explained.

2.3.1. Sample Saturation

The sandstone particle and mudstone particle should be mixed uniformly according to a designed weight ratio with a total mass of 1800 g (Table 1). Then, place the granular particle into the cylinder tube. The SMM sample should be saturated before testing with low hydraulic pressure.

2.3.2. Sample Height Control

In order to achieve the designed height level, a steady axial pressure was conducted by the hydraulic cylinder and oil pump, and sample height in the cylindrical tube was measured by the LVDT. The porosity of the SMM samples at each height level could be calculated as follows:where and are the density and mass of the sandstone particle, respectively; and are the density and mass of the mudstone particle, respectively; is the radius of cross section of the cylinder tube.

2.3.3. Forming Steady Seepage

During the experimental investigation of seepage process, water was provided by a constant hydraulic pressure or flow velocity. When hydraulic pressure gradually became steady value along with time, the gradient of the curve of the seepage with time was a constant value, and then steady flow had formed in the SMM sample. Furthermore, to obtain stable readings, the axial loading was maintained at the setting level firstly. Then start the water seepage measurement, namely, the loading change and seepage flow, which belong to two different test stages.

2.3.4. Experimental Data Recorded

Until a steady seepage was acquired, the water flow velocity and hydraulic pressure with time at each testing step were recorded.

2.3.5. Calculation of the Non-Darcy Hydraulic Properties

The Forchheimer equation [19] has been verified to model the water seepage in granular rocks with good performance [23]. It was applied here to model the relation between the seepage velocity and the hydraulic pressure of SMM sample aswhere is the hydraulic gradient, is the hydraulic pressure, is the water viscosity, is the seepage velocity. is the water density, is the hydraulic permeability of the SMM sample, and is the non-Darcy coefficient. It must be noted that the seepage obeys Darcy law if the value of was close to zero obtained from the test.

In the test, the steady seepage method was adopted to measure the hydraulic properties of SMM sample. As shown in Figure 5(a), the hydraulic pressure at the upstream and that at the downstream ends were signed as and , respectively. The upstream end of the SMM sample was linked with air, that is, . can be automatically recorded with time , which was a steady value and assumed as . Therefore, the hydraulic pressure gradient was a constant. We obtainwhere is the length of SMM sample. If a steady seepage dominates, then

Hydraulic properties (permeability κ and non-Darcy coefficient ζ) were calculated by a least-squares regression method from the relationship between hydraulic gradient and seepage velocity curve (- curve); that is,where is total seepage times of each sample height, , in the test.

2.3.6. Particle Sieve of the Remainder Sample

The screening test sieve machine was applied to monitor PSD mixture alterations during tests; weights of the different particle size ranges in the samples were determined before and after the testing.

Figure 6 illustrates the complete testing procedure. Meanwhile, five groups of SMM samples were tested one by one. To keep the accuracy, each of the SMM samples was carried out for three times, and the average experimental value was determined as the final test results.

3. Test Results and Discussion

3.1. PSD Alterations during Tests

PSD weight of tested rock particles, before and after experiment (loading and seepage), is presented in Figure 7(a). It can be seen that compared with the initial PSD before test, the mean particle size of the sample decreases with the decrease of the weight ratio of SP, which means sample within larger MP weight ratio has a bigger variation rate of PSD. As can be seen from Figure 7(b), the weight of size (20–10 mm) decreases and size (8–0 mm) increases for each sample, respectively. This is because of the crush of larger particles during loading; the weight of smaller particles increases. Moreover, since larger particle sizes can be crushed by loading more easily, the weight of larger size has a more rapid decrease; for example, the largest size (20–15 mm) tends to the most rapid decrease in the size mixture. In detail, the decreased position of weight variation for the sample within less weight ratio of SP was greater than that of MP. This is mainly because MP is crushed more easily than SP during loading, which can also explain why the smaller particle gains more size within less weight ratio of SP. The tiny minority of size (10–8 mm) increased for samples within SP/MP weight ratio of 0 : 1 and 2 : 1 but reduced for the others, indicating that particle crushing can be occurred by fracture and loading. Furthermore, the variation of PSD in Figure 7(a) also shows the water erosion effect by seepage and the weight loss of some fine particles, which were also observed in Figure 7(b). Weight loss increases with decreasing weight ratio of SP. The weight loss for sample within SP/MP weight ratio of 0 : 1 is maximum, that is, 68.7 g, but that of 1 : 0 is minimum, that is, 17.6 g. The former value is approximately 3.9 times the latter one. This indicates that mudstone has higher variation under water erosion than sandstone. In summary, the PSD variation was induced by loading and water seepage. PSD variation of mudstone is higher than that of sandstone. This indicates that particle crushing and erosion properties of mudstone were higher than those of sandstone.

3.2. Porosity Variation during Loading

Although there is some weight loss in each sample because of water erosion, the largest value of weight loss displayed in sample within weight ratio 100% of MP, that is, 56.7 g, is only 3.8% of the original weight (1800 g). Therefore, the effect of weight loss on porosity calculation can be ignored. According to (1), the porosity of samples at different height levels and its variation rate are presented in Figures 8(a) and 8(b), respectively. The results in Figure 8(a) show that the porosity of SMM sample is intensely affected by loading (sample height) and weight ratio of SP. To one given sample, the fracture and void in SMM sample become less during loading; therefore, the porosity decreases with the decrease of sample height. To the same loading (sample height) level, the porosity decreases with the decreasing weight ratio of SP; the sample within SP/MP weight ratio of 1 : 0 is maximum, but that of 0 : 1 is minimum. The density of SP is larger than that of MP, which means the volume of fracture and void in SP mixtures is greater than that in MP mixtures when they have the same weight. As shown in Figure 8(b), the porosity variation rate in sample within less weight ratio of SP is greater, which indicates the rock weight ratio is significant in loading; that is, at the sample height 110 mm, the porosity variation rates of each sample with variable SP/MP weight ratio are as follows: sample 1 : 0 (0.568), 2 : 1 (0.538), 1 : 1 (0.523), 1 : 2 (0.51), and 0 : 1 (0.483), respectively; that means the porosity variation induced by loading to sample within less weight ratio of SP has a bigger sensibility than that of MP.

3.3. Non-Darcy Hydraulic Properties

Relations of hydraulic gradient with seepage velocity in height of 140 mm of each sample are shown in Figure 9. It is obvious that relations are highly nonlinear for the tested SMM samples. Hydraulic gradient decreases with the decrease of seepage velocity, and the decreasing trend becomes slower at a lower seepage velocity. In general, hydraulic gradient-seepage velocity curves display quadratic relation, which indicates Forchheimer equation (2); that is, the non-Darcy model could be applied to simulate the water seepage in SMM samples.

According to (5a) and (5b), the variation of non-Darcy hydraulic properties (permeability κ and non-Darcy coefficient ζ) of SMM samples during loading (variable porosity ϕ) is presented in Figure 10. It is indicated that the decreasing permeability κ and the increasing non-Darcy coefficient ζ are caused by axial loading, mainly due to the absence and narrowing of fracture and void during the loading process.

However, in the loading process, due to the readjustment of the particle structure and the crushing of corners and edges, water flow (fracture) paths become more inconstant; as a result, there is an exception of the abnormal climbing in permeability κ and declining in non-Darcy coefficient ζ; this means some fluctuations are obtained for both κ and ζ; it is mainly resulted from the accession of granular particles in expanded fractures and accelerated seepage. In addition, compared with larger porosity, the sample within smaller porosity seems to have a higher probability of the fluctuation. A rational explanation is that samples at earlier loading stage own much more effective fractures and voids among the rock particles that decrease κ and increase ζ rapidly. However, at the last loading stage, isolated fractures and voids in some samples are able to achieve the changeover from self-extension to interconnection, which causes the tendency fluctuation of κ and ζ, indicating the increase of permeability and the decrease of non−Darcy coefficient. Furthermore, the fluctuation in sample within more weight ratio of MP is greater; for example, there exists fluctuation of κ in sample within SP/MP weight ratio of 1 : 2 and 0 : 1 (Figures 10(d) and 10(e)) at the last loading stage; the fluctuation of ζ in sample within sandstone weight ratio 0 is more obvious than that in others. This indicates that the rock particle crushing and fractures and voids from self-expansion to interconnection in mudstone have higher ability than that in sandstone during loading.

For the same sample height, due to the higher porosity of the sample within more weight ratio of SP, the permeability is larger and non−Darcy coefficient is smaller than other samples. For example, when the sample height is 140 mm, as shown in Figure 10(a), the permeability κ of sample within sandstone weight ratio 100% is 3.7 times greater than that within sandstone weight ratio 0 in Figure 10(e), yet the non-Darcy coefficient ζ is one order of smaller magnitude. Considering the variation of hydraulic properties of the sample within the same particle size and SP/MP weight ratio, it can be concluded that groundwater inrush hazard shows a higher probability of occurrence in sandstone strata and crushed zone (e.g., faults). To summary, hydraulic properties of SMM samples can be affected by sample height (porosity), SP/MP weight ratio, and fractures and voids among the rock particles, and the water flow (fracture) paths cannot be ignored.

3.4. Prediction of Hydraulic Properties Evolution

Due to the inaccessibility to investigate hydraulic properties of granular rocks underground, hydraulic properties are often predicted using limited experimental data. Wang et al. 2017 [24] analyzed the relationship between hydraulic conductivity and PSD. The grain diameter at 10% passing () was regarded in porosity effect. However, cannot be measured in our test because of particle crushing by loading and fine particle loss by water erosion. Therefore, it is better to predict the permeability tendency considering porosity only. To predict the permeability evolution in sands and sandstones, Chilingar [25] investigated a mathematical model to simulate the relationship between permeability and porosity:A good accuracy of the model has been proved in granular mudstone [15]; it is therefore used to model the relation of permeability and porosity in SMM sample in this study. According to (6), we obtainwhere is the effective diameter of particle sizes; and are the permeability and porosity when sample height is less than 140 mm, respectively; and and () are the original permeability and porosity of the sample with a height of 140 mm, respectively.

Furthermore, the non-Darcy coefficient ζ in (2) is an inherent property of the porous media. ζ was also called velocity coefficient [26], inertial resistance coefficient [27, 28], turbulence factor [29, 30], inertial coefficient [31], and so on. Although ζ can be determined by (5a) and (5b), it is important to hold that ζ is a function of the porous medium such as κ and ϕ [32]. Chen et al. [33] and Zhou et al. [34, 35] confirmed that the relation between κ and ζ follows a power law in the form ofwhere λ is the nondimensional coefficient and β is the index. Considering the principle of dimensional consistency and convenience for calculation of engineering applications, an empirical correlation between the κ and ζ is used [36, 37] as where is the Forchheimer coefficient, which can be determined according to the hydraulic experiments in porous media with variable PSD and pore connectivity [38, 39]. Considering (7), we can get a similar form as

Tested and predictive values of hydraulic properties (permeability κ and non-Darcy coefficient ζ) were presented in Figure 11. The model performance was analyzed by using a statistical efficiency criterion, which was on account of the determination coefficient :where is the total number (sample height) of experimental data, . , , and are the experimental value, model predictive value related to a test, and determination factor of permeability κ, respectively. , , and are the experimental value, model predictive value related to a test, and determination factor of non-Darcy coefficient ζ, respectively.

The results of values for variable test sample are shown in Figure 12. It can be seen from Figure 12 that values of all samples are larger than 0.91, which shows a high level of predictive accuracy. Therefore, the adoption of porosity calculation (voids and fracture) for the prediction of hydraulic properties variations presents adequately accurate results. However, as presented in Figure 11, experimental values of hydraulic properties (permeability and non-Darcy coefficient) are typically higher than those obtained from the predictive model by the empirical equation. This could indicate that the increase of hydraulic properties during loading was caused by the particle crushing and water erosion. In (2), the power of non-Darcy coefficient is quadratic, but that of permeability is linear. Therefore, as shown in Figure 12, the average value of permeability (0.977) is larger than that of non-Darcy coefficient (0.942). Moreover, the value of permeability decreases with decreasing weight ratio of SP; for example, value of permeability within weight ratio 100% of SP is the largest (0.991), while that of MP is the lowest (0.936). This suggests that a larger weight ratio of mudstone particles led to a higher percentage of crushed solids, giving greater volume reduction of voids and fracture and altered hydraulic properties. This indicates the particle crushing of mudstone is higher than that of sandstone. Therefore, the predictive accuracy of permeability increases with the increasing weight ratio of SP. Compared with the predictive trend of permeability, non-Darcy coefficient seems to have an opposite trend except the sample within weight ratio 100% of MP. This is mainly resulting from two parameters which are shown in multiplication form in (9).

4. Conclusions

The strata of sandstone and mudstone in Xiaojihan Coal Mine can be naturally filled in the caved zone of mixture state during longwall mining. The caved zone (gob) has a high porosity and permeability resulting in a mass of groundwater storage, which causes a threat of groundwater inrush hazard to the safety of coal mines. To investigate the effect of SP and MP weight ratio on hazard mechanism, a self-designed granular rock experimental system has been conducted to monitor the non-Darcy hydraulic properties of SMM.

The PSD variation before and after test of each sample indicates the PSD variation was induced by loading and water seepage. The weight variation of variable particle size reveals that the sample within more weight ratio of MP has a bigger variation rate of PSD. The weight loss observed in Figure 7(b) indicates mudstone has a higher variation of water erosion than sandstone. In summary, the particle crushing and erosion properties of mudstone were higher than those of sandstone. Furthermore, the effect of weight loss on porosity calculation can be ignored, owing to the large amount of original weight and relatively tiny variation in weight during the test. Porosity evolution of SMM is strongly influenced by loading (sample height) and SP and MP weight ratio. Since the void in SMM sample becomes less during loading, porosity decreases with the decreasing height of the sample. Porosity increases with the increasing weight ratio of SP at the same loading (sample height) level; this is because the density of SP is larger than that of MP. Porosity variation induced by loading to sample within less weight ratio of SP has a bigger sensibility than that of MP.

According to experimental observations, hydraulic gradient-seepage velocity curves show Forchheimer equation (2); that is, the non-Darcy model can be applied to simulate the water seepage in SMM samples. Non-Darcy hydraulic properties (permeability κ and non-Darcy coefficient ζ) of SMM samples are strongly influenced by the SP and MP weight ratio during loading (decrease of porosity ϕ). The decrease of κ and increase of ζ are caused by the absence and narrowing of fracture and void during loading. Moreover, because of the accession of granular particles in expanded fractures and accelerated seepage, κ-ϕ and ζ-ϕ curves show some local fluctuation. The sample within smaller porosity or within more weight ratio of MP seems to have greater fluctuation. Isolated fractures and voids are able to achieve the changeover from self-extension to interconnection at the last loading stage, which cause the fluctuation tendency of κ and ζ. Fluctuation ability in mudstone is higher than that in sandstone during loading. Considering the variation of the hydraulic properties for the rock particle mixtures within the same particle size and SP/MP weight ratio, it can be concluded that groundwater inrush hazard seems to have a higher probability of occurrence in sandstone strata and crushed zone (e.g., faults).

To predict hydraulic properties (permeability κ and non-Darcy coefficient ζ) using limited experimental data, an empirical model was used to analyze the hydraulic properties evolution for the crushing and seepage courses. Particle crushing and water erosion that cause increasing hydraulic properties during loading are the main reason that experimental values of hydraulic properties are typically higher than those obtained from the predictive model. values indicate a high degree of predictive accuracy of hydraulic properties. As the power of non-Darcy coefficient is quadratic but that of permeability is linear, permeability κ has a higher predictive accuracy than non-Darcy coefficient ζ. Furthermore, the predictive accuracy of permeability increases with the increasing weight ratio of SP; however, non-Darcy coefficient has an opposite trend.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work is supported by the National Basic Research Program of China (2015CB060200), the National Natural Science Foundation of China (41772313, 51322403, and 51274254), the Key Research and Development Program of Hunan (2016SK2003), the Natural Science Foundation of Hunan (2015JJ4067), and the Philosophy and Social Sciences Foundation of Hunan (14YBA400). The first author would like to thank the financial support by the Open Project of Research Center of Coal Resources Safe Mining and Clean Utilization, Liaoning (LNTU17KF03).