Surface Subsidence Modelling Induced by Formation of Cavities in Underground Coal Gasification

Abstract

Underground coal gasification (UCG) is an efficient method for the conversion of deep coal resources into energy. The scope of this work is to model the subsidence of four gasification cavities with a size of 30 m × 30 m × 15 m, separated by 15 m wide pillars. Two scenarios of gasification sequence are modelled, one with the gasification of cavities 1 and 2 followed by 3 and 4, and the other one with the sequence of cavities 1 and 3, followed by 2 and 4. The results show that the final surface subsidence after gasification of four cavities is 9.8 mm and the gasification sequence has an impact only on the subsidence at the intermediate stage but has no impact on the final subsidence after all four cavities are formed, when only the elasticity regime is considered. Additionally, the maximum surface subsidence for the studied cavities of different sizes ranges from 0.016 mm to 7.14 mm, and the relationship between the subsidence and the cavity volume is approximately linear. Finally, a prediction model of surface subsidence deformation is built up using the elastic plate theory, and the formula of surface deformation at a random point is given. The maximum difference between measured and calculated deformation is 4.6%, demonstrating that the proposed method can be used to predict the ground subsidence induced by UCG.

1. Introduction

Underground coal gasification (UCG) is a process in which coal is converted into product gas without mining by artificially enhancing gas permeability in a section of a coal seam, igniting the coal, partially combusting and gasifying coal by means of injected oxidants, and producing the product gas for cleanup and processing for a variety of end uses. During the UCG process, the volume of the cavity increases progressively with coal consumption and thermo-mechanical spalling, if any, from the roof. The formation of an underground cavity will lead to the redistribution of stress field near the cavity and the disturbance will lead to surface subsidence. Hence, it is significant to investigate the effect of the formation of a cavity on surface subsidence.

Avasthi et al. [1] more accurately simulated surface subsidence during UCG by introducing new elastic parameters and empirical correction factors. Ma et al. [2] studied the issues of surface settlement due to the mining of a horizontal coal seam by using the elastic thin plate theory. Wang et al. [3] proposed a new thermal–mechanical–chemical displacement multifield coupled numerical simulation method, which was based on controllable pumping injection point (CRIP) technology to numerically simulate the coal seam gasification working face, and analyzed the evolution law of temperature and displacement fields after the gasification reaction. Li et al. [4] applied the complex function analysis method to solve the stress field of an elliptical formation cavity under a self-weight stress field and a vertically uniformly distributed load. By developing a thermomechanical coupling model, Yang et al. [5] simulated heat transfer, stress distribution, and surface subsidence during underground coal seam gasification using ABAQUS. The material properties of rocks and coal were obtained from the existing literature and geomechanical tests were carried out on samples derived from the demonstration site in Bulgaria. Three days of gasification ere simulated by assigning a moving heat flux on a cell of size 2 m × 2 m × 2 m at a velocity of 2 m/day. Results of temperature and stress distribution showed that the developed numerical model was able to simulate the heat propagation and the stress distribution around cavities under a thermal–mechanical coupled loading during the UCG process. Also, the surface subsidence was found to be 0.08 mm after three days of gasification for the case studied. Elahi et al. [6] studied the influence of different coal seam constitutive models on the thermal–mechanical response during the UCG process by coupling thermodynamic and geomechanical modules and used stress rebalancing technology to simulate the elastic–plastic behavior of coal seams. McMahon et al. [7] proposed an elliptical cave expansion model for estimating the bearing capacity and settlement of circular shallow foundations on clay. Yoo et al. [8] developed a soil flow protector (SFP) through field experiments and numerical analysis to reduce soft soil settlement caused by cave formation. Zha et al. [9] simulated the movement of a coal seam roof and surface subsidence during shallow underground coal gasification process using a finite element model coupled with a heat transfer module in COMSOL. Ekneligoda et al. [10] developed a coupled thermomechanical numerical model to simulate the underground coal seam gasification (UCG) process and simulated and predicted the ground settlement of the test site under different working conditions. Zhao et al. [11] discussed the risk assessment and protective measures for ground subsidence caused by underground mining. Lee et al. [12] evaluated the impact of cave types on ground stability through a finite element analysis. Made et al. [13] proposed an approximate solution for predicting the settlement of underground cylindrical caves based on empirical and analytical methods. Liu et al. [14] studied the rock movement characteristics and surface settlement prediction method during underground coal gasification (UCG) process. Shubham et al. [15] used surrogate models to predict settlement above foundation caves and conducted a reliability analysis. Jing et al. [16] studied the mechanical mechanism and prediction method of the settlement mechanism through numerical simulation and theoretical analysis. Rezaei et al. [17] used numerical analysis (NA), neural network (NN), fuzzy logic (FL), and statistical analysis (SA) models to estimate the vertical displacement of the middle and key points of the roof and floor of a power station cave under different conditions. Cui et al. [18] discussed the influence of voids in the strata on soil deformation during shield tunnel construction and established a three-dimensional symmetric calculation model to calculate soil deformation and surface settlement. Derbin et al. [19] used FLAC3D for numerical modeling and studied the mechanism of surface subsidence during underground coal gasification, taking into account the influence of thermal effects on soil rock mechanical properties. The impact of mining on surface streams was studied by Agioutantis et al. [20] based on self-developed SDP software. In addition, some people have studied the gasification process of underground coal, explored the development of underground cavities and their sustained effects on surrounding rock layers [21,22,23]. In order to more accurately predict the process of surface subsidence, some scholars continuously improve the dynamic prediction model of underground mining subsidence in order to achieve better prediction results [24,25,26,27,28,29,30].

The scope of this work is to model, under the elasticity regime, the surface subsidence when multiple gasification cavities are formed within a coal seam at a depth of 300 m. Four cavities (C1, C2, C3, and C4) of size 30 m × 30 m × 15 m (13,500 m³) separated by 15 m wide pillars are modelled. It is also intended in this work to investigate the influences of the gasification sequence on the surface subsidence based on two different scenarios: (I). C1 and C2 first, followed by C3 and C4; (II). C1 and C3 first, followed by C2 and C4. These cavities are within an 18-m-thick coal seam overlain by carbonaceous mudstone with occasional siderite bands.

2. Model Configuration

2.1. UCG Project

The coal mine at Leigh Creek, South Australia, is exploited by UCG technology. Similar to longwall coal mining, UCG also faces the risk of overburden caving. Because UCG does not have workers down the well and other processes, the cavity tunnel does not have any support. However, the collapse of overlying strata intensifies the land subsidence and bring negative effects to the surrounding environment. Hence, a surface subsidence evaluation is essential for UCG.

2.2. Model Setup

Due to the stress redistribution around a cavity, in the numerical models, the distance between model boundaries and cavities should be large enough to reduce the boundary effects. A sensitivity analysis was conducted in this study and this distance was found to be at least ten times the size of the disturbed area for the geological conditions of this project. All numerical models set up in this study had the same dimensions of x × y × z = 3465 × 1665 × 418 m³, which ensured the distance (1650 m) between the front, back, left, and right boundaries of the cavities was ten times the size of the final disturbed area (4 × 30 + 3 × 15 = 165 m), as shown in Figure 1. Note, due to symmetry, only 1/2 of the full model was used. For boundary conditions, no x displacements were assumed at the left and right boundaries (perpendicular to the x axis) and no y displacements were assumed at the front and back boundaries (perpendicular to the y axis). This was to simulate the actual condition that the horizontal disturbance caused by the cavity only had a limited zone of influence. No vertical (z) displacements were assumed at the bottom boundary of the model, and the top boundary was a free boundary. The geomechanical properties of the rocks considered in the model are listed in

Table 1. Geomechanical properties of rocks.

Property	Quaternary	Main Series (Mudstone)	Coal	Lower Series Overburden
Density (kg/m3)	2040	2040	1040	2040
Young’s modulus E (GPa)	3.4	3.4	0.785	3.4
Poisson’s ratio ν	0.22	0.22	0.25	0.22

, which are sourced from Lucas and TriLab [31,32]. The commercial code ANSYS was used to calculate stress field and deformation field based on this model setup. ANSYS code is characterized by strong multi-physical modeling, solution, and nonlinear analysis.

To investigate the effect of the gasification sequence, two scenarios were simulated, as listed in

Table 2. Scenarios of gasification sequence investigated in this study.

Scenario	Stage	Gasification
I	1	Cavity 1	Cavity 2
	2	Cavity 3	Cavity 4
II	1	Cavity 1	Cavity 3
	2	Cavity 2	Cavity 4

. In Scenario I, the gasification of cavities 1 and 2 happened first (stage 1), followed by the gasification of cavities 3 and 4 (stage 2). In Scenario II, the gasification of cavities 1 and 3 was completed first (stage 1) before that of cavities 2 and 4 (stage 2).

3. Simulation Results

3.1. Stress Distributions

Figure 2 shows the vertical stress profile generated by the gravitational force, indicating that the vertical stress was about 6.1 MPa at a depth of 300 m. Figure 3 and Figure 4 show the distributions of maximum and minimum principal stresses of the two stages for Scenario I. The corresponding distributions for Scenario II are shown in Figure 5 and Figure 6. As expected, there were significant tensile stresses developed at both the roof and the floor of the cavities, and the maximum compressive stresses appeared in the side walls. For Scenario I, the maximum tensile stresses were 1.32 and 1.42 MPa, and the maximum compressive stresses were 11.82 and 11.90 MPa, for stage 1 and stage 2, respectively. For Scenario II, these values were 1.25 and 1.42 MPa and 11.76 and 11.90 MPa for stage 1 and stage 2, respectively. Compressive and tensile stresses at stage 1 for Scenario I were slightly higher than those of Scenario II. However, the values after four cavities were completely formed for both scenarios were identical. In other words, the gasification sequence only affected the stresses at the intermediate stage but not the final stress distributions. Note the validity of this conclusion is clearly limited only to the elasticity regime without the consideration of plasticity and failure during the gasification process (cf. Section 3).

3.2. Subsidence

Figure 7 and Figure 8 shows the z-direction displacement of the area around the cavities for Scenario I and Scenario II, respectively. As expected, the maximum z-direction displacements occurred at the roof (negative value) and floor (positive value) of the cavity. The maximum z-direction displacements at the roof were 154.2 and 156.0 mm for stage 1 and stage 2 in Scenario I, and 147.4 and 156.0 mm for stage 1 and stage 2 in Scenario II.

Figure 9 shows the z-direction displacement profiles on horizontal cross sections at different depths from the surface. It can be seen that at the depth near the cavities, the z-direction displacement decreased quickly as it moved away from the center position. The z displacements on a cross section became more and more uniform as the section moved further away from the cavities. At the surface (z = 0 m), the z-direction displacement of the cross section was the modelled surface subsidence, which showed a profile with small variations (Figure 9), suggesting the impact of cavities on the subsidence profile was significantly reduced for areas far away from the cavities.

Figure 10 show the maximum z-direction displacement at different depths from the surface. This position is at the center line through subsidence. Figure 11 shows the z-direction displacement contour plots on the horizontal cross section at z = 0 m (note the different color scales). After four cavities were formed, the same z-direction displacement profiles were observed for Scenario I and Scenario II and they followed approximately an upside-down elongated bell shape, with the maximum z-direction displacement in its center (shown in Figure 11). The maximum z-direction displacement was 9.8 mm, occurring at the center position (x = 82.5, y = 0). However, the subsidence profiles at stage 1 for the two scenarios were different. The maximum subsidence was 5.0 mm for Scenario I and 4.7 mm for Scenario II at the completion of the stage 1 gasification. The surface subsidence curves along the central horizontal line at y = 0 m are displayed in Figure 12. Clearly, the surface subsidence profiles were asymmetrical in stage 1 (largest on top of stage 1 gasification cavities) but became symmetrical in stage 2 for both scenarios (note the center point is at x = 82.5 m). The variability of the subsidence of stage 1 for Scenario I was higher compared with that for Scenario II.

Additionally, the effects of the cavity size on surface subsidence were also investigated. One-cavity models with 5.3 × 5.3 × 5.3 m, 15 × 15 × 15 m, 20 × 20 × 20 m, 22 × 22 × 15 m, 26 × 26 × 15 m and 30 × 30 × 15 m were built and modeled, respectively. Figure 13 shows the maximum z-direction displacement on different horizontal cross sections at the center point (x = 0, y = 0) as a function of depth from the surface. As can be seen, major z-direction displacements occurred between 400 and 500 m. This suggests that the major impact on subsidence caused by the cavity was only up to about 100 m above the cavity, which can be considered as the critical range of influence in this study. Clearly, the greater the cavity size, the greater the critical range of influence, as demonstrated in Figure 13. For the actual value of surface subsidence (x = 0, y = 0 and z = 0 m), as expected, it increased as the volume of cavity increased, as shown in Figure 14. For the smallest cavity (152 m³), the surface subsidence was relatively small at only 0.016 mm. For the largest cavity (13,500 m³), the surface subsidence reached 7.14 mm, which was around a 450-fold increase compared with that of the smallest cavity, though the increase in volume was only around 90 folds. The relationship was approximately linear as demonstrated in Figure 14, but it is expected that it will become nonlinear as the increasing rate accelerates as the cavity size becomes significantly larger.

4. Discussions

4.1. Numerical Simulation Discussion

In this study, the modelling was carried out under the elasticity regime, and no rock failure was considered. In reality, due to stress concentrations around cavities and significant cavity roof deformations as modelled, as shown in the simulation results given above, parts of the cavity roof or walls would fail. In other words, the cavities are not expected to stay stable after they are created. Under the action of the vertical stress, caving of the cavity will develop upwards. However, this caving development is limited due to two reasons: (1). from the stress re-distribution perspectives, the caving eventually stops after it forms a stable arch which is strong enough to support the redistributed stresses; (2). there is a volume swelling factor that needs to be considered after the roof caves. The cavity voids are eventually filled by caved rocks, which help prevent further caving development. From the subsidence perspectives, cavity caving, if it happens, may actually reduce the surface subsidence. This is due to the caving possibly causing a more favorite stress redistribution around the cavity. Since no failure parameter data of rock and coal were obtained, the subsidence analysis was not conducted with an elastoplastic model. However, the subsidence calculation would be more adequate using a rock failure model, such as Drucker–Prager or Mohr–Coulomb models.

This study also shows that the final stress and subsidence distribution profiles are identical between two gasification scenarios after all four cavities are formed. This is only true under the elasticity regime. If rock failure and cavity caving are considered, due to significant stress redistributions and alteration of cavity geometries after stage 1 gasification, very different initial conditions are expected for stage 2 operation and therefore, the final stress and subsidence distribution profiles of the two scenarios are also expected to be very different. A detailed modelling considering these issues should be conducted in next work. Also, the order of the gasification of cavities can affect intermediate subsidence and stress distribution, which provides guideline for the order of gasification.

It is noted that massive gas (such as carbon monoxide, methane) would be produced and form cavities by internal pressure during gasification process. This pressure keeps the cavities in the tension condition, which affects the stability of cavities. This need to be considered in a 3D model. However, temperature, internal pressure distribution, and value need to be obtained based on field monitoring. This work will be conducted in the future.

4.2. Surface Subsidence Prediction Model

To determine the subsidence at any point in the surface moving zone due to the UCG process, it is necessary to establish the subsidence equation of the surface subsidence. The established coordinate is shown in Figure 15.

At depth H, the coal seam per unit volume is mined, and a stable subsidence is formed [33,34，35],

$$W_e\left(x,y\right)=W_0\left(\frac{x^2}{2p}+\frac{y^2}{2q}-1\right)$$

(1)

where W₀ is the maximum subsidence value. This value can be determined by

$$W_0=mE\cos \alpha$$

(2)

where m is the thickness of coal seam; E is the subsidence coefficient; and α is the coal seam inclination. p and q are the corresponding strike and dip half-length of the subsidence curve, which can be determined by a geometry relationship or a measured method [1,9].

When mining is carried out at the underground location (x, y, z), the expression of the surface subsidence curve caused by this is as follows:

$$\mathrm{d}W\left(x,y\right)=W_e\left(x,y\right)\mathrm{d}s$$

(3)

where ds is the element area of the upper surface S at (x, y, z) of the underground excavation space Ω.

The equation of S at the upper part of the Ω surface of the underground excavation space is

$$z=z\left(x,y\right)$$

(4)

Hence, the corresponding area ds is

$$\mathrm{d}s=\sqrt{1+z_x^2\left(x,y\right)+z_y^2\left(x,y\right)}\mathrm{dxdy}$$

(5)

When the upper surface of the excavation space is flat, the upper-part surface equation is,

$$z=c c \text{is} \text{constant}$$

(6)

Under this condition, ds = 1. Hence, the surface subsidence prediction model is,

$$W\left(x,y\right)={\displaystyle \underset{\mathsf{\Omega }}{\iint }{W}_e\left(x,y\right)}\mathrm{dxdy}$$

(7)

It is noted that the surface coordinates’ integral extends to the point on the ground where the coordinates are (x, y). When the upper surface of the excavation space is an inclined plane, the upper-part surface equation is,

$$z=ax+by+c$$

(8)

where a, b, and c are plane parameters. The surface subsidence prediction model is

$$W\left(x,y\right)={\displaystyle \underset{D}{\iint }{W}_e\left(x,y\right)\sqrt{1+a^2+b^2}}dxdy$$

(9)

In order to validate the subsidence prediction accuracy of the proposed model, the numerical simulation results of four cavities were compared to that of the proposed model (Equation (9)). Figure 14 shows the comparison of the subsidence of the numerical simulation with that of the proposed model.

The calculated maximum subsidence value was 10.2 mm, which was a small difference from the numerical value. The maximum difference between numerical and calculated deformations was 4.6%, which is acceptable for engineering analysis. Additionally, the numerical surface subsidence was greater than the calculated value, plotted in Figure 16. The reason was that the interaction between cavities was not considered in the proposed model. Actually, the interaction between cavities would affect the surrounding ground deformation. Hence, it is indicated that the proposed method can be used to analyze the ground deformation induced by UCG.

In order to illustrate the proposed model for actual mine gasification, an example was given with three gasification cavities, taken from Avasthi and Harloff, (1982) [1]. Data for this example were normalized. Figure 17 shows the comparison of the proposed model and field data. The curve tendency of the proposed model was the same as that of the field data. It is noted that the two profiles were closest over the first cavity, whereas large discrepancies were observed over the other cavities. The whole error between the two curves was less than 17.8%. Hence, the proposed model can be used to calculate subsidence profiles for multiple cavities.

5. Conclusions

In this study the subsidence corresponding to the gasification of four cavities (C1, C2, C3, and C4) of size of 30 m × 30 m × 15 m, separated by 15 m wide pillars at a depth of 300 m was modeled under the elasticity regime. Two gasification scenarios were investigated, including Scenario I: gasification of C1 and C2 (stage 1), followed by C3 and C4 (stage 2); and Scenario II: gasification of C1 and C3 (stage 1), followed by C2 and C2 (stage 2). For both scenarios, the final maximum surface subsidence after the completion of all four cavities was 9.8 mm and the maximum compressive stress and tensile stress around cavities were 11.9 MPa and 1.42 MPa. The surface subsidence profile was an approximate upside-down elongated bell shape with the maximum subsidence occurring at the center position of the four cavities (x = 82.5 and y = 0). The maximum z-direction displacement at the cavity roof was 156.0 mm and the z-direction displacement profile at the horizontal section directly above the cavities was an upside-down bell shape. The variations in z-direction displacements on horizontal cross sections became more and more uniform as the section moved further away from the cavity and got closer to the ground surface.

The gasification scenarios (sequences) only affected the stress and subsidence distributions at the intermediate stage (stage 1), not the final results after all four cavities were formed, which was due to the elasticity regime considered. For Scenario I, at stage 1, the maximum surface subsidence was 5.0 mm, the maximum z-direction displacement at the cavity roof was 154.2 mm, the maximum compressive and tensile stresses around the cavities were 11.82 MPa and 1.32 MPa, respectively. For Scenario II at stage 1, the maximum surface subsidence was 4.7 mm, the maximum z-direction displacement at cavity roof was 147.4 mm, the maximum compressive and tensile stresses around the cavities were 11.76 MPa and 1.25 MPa, respectively.

For the smallest cavity (152 m³) modelled, the maximum surface subsidence was 0.016 mm, while for the largest cavity (30 × 30 × 15 m³), the maximum surface subsidence reached 7.14 mm. The relationship between surface subsidence and cavity volume in this case was approximately linear, though a nonlinear relationship is expected at larger cavity volumes.

A prediction model of surface subsidence deformation was built up using the elastic plate theory, and the formula of the surface deformation at a random point was given. The comparison between numerical results and the proposed model demonstrated that the proposed method can be used to predict the ground subsidence induced by UCG.