Numerical investigations to assess ground subsidence and fault reactivation during underground coal gasification

This paper aims to assess potential environmental impacts associated with commercial-scale application of in situ coal conversion in a target area in Hungary. The site is an environmentally protected forested area. Parametric numerical modelling techniques were employed using a discrete element method (DEM) to evaluate the surface subsidence and fault activation during mining processes. The Mohr–Coulomb elastic–plastic material model adopted to simulate rock formations. Zero-thickness interfaces with friction and cohesive characteristics were employed to simulate geological faults. A sensitivity study on rock formations and fault properties was conducted to address the importance of geological parameters that have an impact on surface subsidence as well as fault activation, which could in turn result in pollutant migration from the deep coal seams to the surface. The analysis of the results demonstrated that surface subsidence is affected by the average Young's modulus of the geological strata, whereby the activation of the faults is influenced by the friction angle between faults. Also, it was found that shallower seams are more likely to produce surface subsidence; i.e. as excavation depth increases, the surface subsidence decreases. Finally, computational outputs from this work were used to develop the web-based and interactive Environmental Hazards and Risk Management Toolkit (EHRM) for planning and decision-making processes during in situ coal conversion.

1 Introduction

Conventional and unconventional mining may trigger several environmental risks to human health and the environment in both the short and long term. Subsurface mining could lead to tunnel collapses and ground surface subsidence. Also, as mine water takes on harmful concentrations of minerals and heavy metals, it becomes a contaminant. This contaminated water can pollute the surroundings of the mine and beyond. Similarly, during in situ coal conversion, uncontrolled migration of synthesis gas to the surface can result in adverse impacts on adjacent aquifers and local ecosystems. As coal is progressively converted, the overlying strata become unsupported, resulting in surface subsidence. Also, subsurface mining, if not controlled, could trigger near-field seismicity due to fault reactivation and develop potential pathways for pollutant migration to the surface, exposing environmental risks. Therefore, understanding the nature and extent of potential contamination as well as quantification of the potential risks posed to human health and the environment are of paramount importance for sustainable mining operations. With the continuous and rapid development of computer and information technology, various numerical methods have been applied in the past to study surface subsidence and fault activation during mining. For example, Shoemaker et al. (1978) studied the effect of thermomechanical properties of overburden coal seams associated with underground coal gasification (UCG) on ground subsidence, while Advani et al. (1983) used numerical modelling to determine the ranges of UCG cavity shape, roof collapse, surface subsidence, pore pressure, and creep response at the overburden rock formations and coal. Vorobiev et al. (2008) presented results of geo-mechanical processes induced by UCG activities using the combination of continuum and discrete approaches for predicting ground subsidence. Tian et al. (2009) used numerical thermomechanical models to evaluate the temperature impact on ground subsidence, while a few years later, Tian (2013) studied the UCG-induced ground subsidence in controlled-retraction and injection point (CRIP) configuration by application of a thermomechanical model using Abaqus software. The same year, Sarhosis et al. (2013) developed a coupled hydrothermal analysis for UCG reactor evaluation, while a year after, Yang et al. (2014) developed a thermal mechanical coupled analysis using Abaqus software to study heat transfer, stress distributions around a UCG cavity, and the consequent surface subsidence. Also, Kiryukhin et al. (2014) investigated surface subsidence in geothermal fields through advanced finite-element models (FEMs), while Otto et al. (2016) developed a regional-scale thermomechanical 3D model to understand the fault reactivation process during UCG operations. More recently, Kempka et al. (2016) presented the flexible simulation framework to couple processes in complex 3D models for subsurface utilization assessment, while Kempka et al. (2022) investigated the probability of contaminant migration from in situ coal conversion reactors due to fault reactivation and ground subsidence. From the above, although a lot of work has been done to investigate the mechanical behaviour of rock formations during in situ coal conversion and UCG process, there are limited studies to investigate the influence of geological parameters on the geoenvironmental issues during unconventional mining operations. Thus, this study aims to present the implementation of a sensitivity study, using the commercial three-dimensional numerical modelling code 3DEC developed by Itasca (Itasca, 2021), to assess environmental risks such as ground surface subsidence and water contamination caused by fault reactivation during in situ coal conversion process.

2 Case study: geological setting

The coal deposit under investigation is in the Mecsek mountain region and located approximately 200 km south-west from Budapest, Hungary. The target site contains 20 mineable coal seams with 1.4 to 3.4 m thickness and total coal resources of 438 Mt, with 280 Mt as mineable reserves. One-third of the coal resources are located at depths of up to 600 m below ground level at a dip angle of ≤ 30^∘, whereby two-thirds of the coal is situated at depths up to 1100 m. The complex structured coal deposit consists of ash and sulfuric highly volatile coals with a methane content of 20 to 50 m³ t⁻¹ of coal. Three main types of fault systems (basin-margin normal-growth faults, reverse faults, and late-stage thrust faults) subdivide the target area into several coal blocks; see Fig. 1a and b (Fodor, 2006).

Figure 1 (a) SW–NE cross-section through the geological model of the target area; the considered blocks for in situ coal conversion are marked by numbers and the coal-bearing sequences are displayed in brown, grey, lavender, and green. (b) Sectional view of blocks and coal-bearing sequences considered for in situ coal conversion whereby red blocks were selected for the in situ coal conversion (Fodor, 2006).

In this study, the deeper and less accessible parts of the coal resources were considered for in situ coal conversion. The lower blocks in the south-western part of the license area were selected as favourable for mining. Figure 1a shows the SW–NE cross-section which runs through the target area. From Fig. 1a, it can be observed that the coal seams are increasingly deeper towards the western border of the deposit. The smaller tectonic blocks in the SW are situated in deeper and geologically unfavourable settings. Hence, a perpendicular cross-section trending SW–NE was considered to contain the main geological features, thrust faults, and back-thrust structures. Figure 1b shows the target blocks; the dashed red line indicates that block nos. 3, 5, and 6 were selected for the in situ coal conversion.

3 Numerical model set-up, mining sequences, and excavation process

3.1 Development of the 3D geometry and mesh generation

A numerical model has been developed to represent the 3D geology of the target area, including the coal layers and faults, which divides the coal deposit into a series of individual geological blocks. The geometry was created by merging the survey data and NE–SW cross-section. In total, eight coal seams were created. The coal layers were modelled in three blocks. The thickness of a coal seam ranged from 2.5 to 7.5 m. As shown in Fig. 2a, a significant number of zones and grip points, i.e. 995 517 zones and 174 075 grid points, were generated for the development of the model. Mesh refinement was used around the area of interest (excavation), while coarse meshes were implemented near the model boundary. Faults are represented as zero-thickness interfaces; see Fig. 2b.

Figure 2 (a) Generated 3D geometry of the model considering the mesh refinement around the excavation area; “block ID” refers to the zones generated in the model. (b) Faults in the model represented as zero-thickness interfaces.

3.2 In situ stresses and boundary restraints: setting to initial equilibrium state

Initially, the model's boundary was considered sufficiently large and extended to approximately 1000 m in the south-west direction to eliminate any boundary effects. Also, the boundary conditions of the model were fixed at the four lateral and bottom boundaries. The upper surface of the model was set free. For each case, the model was brought into equilibrium under gravity conditions. To ensure equilibrium conditions, the average force ratio was measured. When the average force ratio was close to zero, then equilibrium conditions were assumed.

3.3 Constitutive behaviour for rock, coal, and faults

Rock and coal were assumed to behave in an elastic–plastic manner according to the Mohr–Coulomb yield criterion. Faults were modelled as zero-thickness interfaces obeying Coulomb sliding constitutive law, which allows for friction, cohesion, dilation, stiffnesses (normal and shear), and bond strength (tensile and shear) characteristics. The normal and shear stiffnesses at the joints along the faults were estimated based on 10 times the equivalent stiffness of the stiffest neighbouring zone. The normal and shear stiffness was set to 10 GPa m⁻¹. The smallest width of an adjoining zone was taken equal to 1 m.

Figure 3 (a) Vertical displacement for fault friction equal to 10^∘ considering the downward to upward scenario (DU); (b) plastic behaviour of coal seams situated in blocks 3, 5, and 6; (c) displacement plot in terms of metres for ground subsidence when Young's modulus was minimum (a2 min); (d) comparative parametric vertical displacement representation when maximum ground surface subsidence was observed.

4 Evaluation of surface subsidence

Since geotechnical data on rock mass and coal were limited, a sensitivity study was undertaken to identify the most sensitive parameters that influence surface subsidence. In addition, the mining excavation sequence was also investigated. In total, 66 cases were evaluated. In the parametric evaluation, seven variables were investigated considering their minimum, average, and maximum values. The parameters varied were the unit weight (a1), Young's modulus (a2), Poisson's ratio (a3), the friction angle (a4), cohesion (a5), and fault friction (ff). In each computational experiment, one parameter was varied at a time, while the others were kept constant. In all cases, the horizontal to vertical stress ratio constant (k) was assumed to be equal to 1. In this way, all parameters were analysed and their sensitivity to the output results was evaluated. Surface subsidence was recorded using a total of 66 monitoring lines above the mining and excavation areas. Each monitoring line had 50 monitoring points. The distance between adjacent monitoring points was taken equal to 100 m. A comparative study was also undertaken to predict the subsidence displacement of multi-seam mining. In addition, two mining sequences were considered including (a) a downward to upward (DU) and (b) an upward to downward (UD) mining sequence. In this way, the effect of the mining sequence on surface subsidence and fault activation was identified. For the DU mining sequence, the excavation process started from the coal seam situated in block 6 and completed at the upper coal seam in block 3. The UD mining sequence started in the upper seam in block 3 and completed in block 6. Figure 3a shows the results of the surface subsidence contours for downward to upward (DU) excavation when the fault friction was set to 10^∘, which was the worst-case scenario. From Fig. 4a, it can be observed that the maximum vertical displacement is related to the minimum Young's modulus (a2min), while the minimum vertical displacement was observed when Poisson's ratio had a maximum value (a3max). Figure 3b presents the plastic deformation around the excavation area. In the target coal deposit, the coal seams were thin enough compared to the overburden material. So, the elastic parameters governed the material deformation and ground subsidence. As shown in Fig. 3b, there is a tensile plastic zone of rock mass around the excavation area connected with the plastic zone of coal pillars as the mining progressed. Moreover, at the ground level there is an elastic zone with a small range of coalescence to the plastic zone in the area above the mining and excavation. The depth and range of the tensile damage zone of the coal and host rock mass are shown in Fig. 3c. The coal extraction width to coal height ratio was taken equal to 1. The pillars that remained between excavations significantly contributed to the stability of the surrounding rock with a relatively smaller tensile damage zone. From Fig. 3c one can observe the increase in the maximum vertical subsidence when Young's modulus decreased. This could be explained by the effect of elastic deformation. Also, when the value of Young's modulus was relatively large, the overburden deformed primarily by concentrated elastic deformation at the overburden. Therefore, it was confirmed that the softer overburden rock undergoes much less plastic strain deformation than the stiffer one. Also, Fig. 3d presents comparative result of the sensitivity study. As can be seen from Fig. 3d, the minimum displacement was related to the minimum Poisson ratio of the studied material in all cases, while the maximum displacement was recorded when Young's modulus was maximum. From Fig. 3d, the excavation sequence does not significantly affect the surface subsidence, while the friction angle of the fault does. The latter is more prominent when the fault friction was considered to be equal to 10^∘, where significant shear displacement between blocks occurred, resulting in vertical displacement at the ground level. However, for higher fault friction values (e.g. 20 and 35^∘), the amount of surface subsidence was relatively small or even negligible. Finally, the maximum and minimum vertical displacements are related to elastic characteristics of the material (minimum Young's modulus and maximum Poisson's ratio) when the friction angle of the fault was 10^∘.

Figure 4 (a) Maximum normal displacement (m) when the friction at the faults was 10^∘, (b) maximum shear displacement (m) when the fault friction was 10^∘, (c) interface normal stress (Pa) when the up-to-down excavation case was considered and when Young's modulus was minimum with the value of friction along the faults of 10^∘, and (d) interface shear stress (Pa) for the up-to-down excavation case when the value of Young's modulus was minimum and the value of friction at the faults was 10^∘.

5 Fault reactivation analyses during the mining process

One of the crucial environmental issues associated with in situ coal conversion is the estimation of fault activation during the mining processes. The possibility for fault activation was evaluated by performing 66 computational simulations with varying mechanical properties of the rock and coal as well as frictional characteristics at the faults. In particular, the fault friction varied from 10 to 35^∘, and the normal and shear displacement stresses along the faults were estimated. Figure 4a shows the outputs of the maximum normal displacement for the fault with numbers 6, 7, 13, and 14 when the friction angle along the faults was equal to 10^∘. Also, Fig. 4b shows the maximum shear displacement at faults 7 and 8 when the fault friction was 10^∘. In Fig. 4c, the maximum normal displacement was observed around block 6, which is characterized by thick coal seams. However, one should bear in mind that the system of blocks is interconnected and there could be block interlocking not allowing them to displace any further irrespective of the frictional characteristics of the faults. Figure 4c and d show the variation of the normal and shear stress in surrounding rock mass during the in situ coal conversion process for fault friction of 10^∘. Also, as expected, during the process of mining near the fault, the change in fault stress near the mining area is mainly caused by the increase in shear stress. During the mining process, the working face and the shear stress on both sides of the excavation area were the same. From the overburden to the coal seam, the initial state of normal and shear stress of each measuring point increased, which conforms to the distribution law of in situ stress. Also, both normal stress and shear stress increased, but the increase in shear stress was larger than that of normal stress. This situation was found mainly around coal seams at blocks 3, 5, and 6. Also, the stress difference on the fault plane was found to relate to the stress evolution of the surrounding rock formation during coal mining. So, the further away from the coal seam, the earlier the stress change starts. As can be seen in Fig. 4c and d, considering the fault plane near the excavation area, the shear stress mainly increased, while the rate of change in normal stress was smother and less sensitive to fault activation and thus subsidence.

6 Conclusions

The main objective of this study was to present the implementation of a probabilistic analysis to assess and mitigate environmental risks while improving the productivity of the UCG process. One of the most important issues that should be considered in this process is the management of underground subsurface mining. The crucial environmental issues associated with the mining process are ground surface subsidence and contamination of water caused by fault reactivation, especially in bad geological conditions. In order to assess the environmental risks, parametric modelling techniques were employed for the coal deposits at Máza–Váralja (Hungary). Also, parameter uncertainty was evaluated to identify the optimum operational conditions concerning potential risks to human health and the environment as well as to increase coal conversion and productivity. In particular, the outputs of the parametric study focused on assessing mine stability, fault integrity, and surface subsidence. In the studied area, the parametric numerical method was employed using discrete element approaches to address the parameters with the most important impact on the mining process regarding enhanced risk, providing some information concerning the factors that affect subsidence and fault reactivation as an output. The simulation results demonstrate that surface subsidence is affected by the average Young's modulus of the geological strata and the fault activation to the friction angle of the faults. Also, shallower seams are more likely to produce surface subsidence, while as excavation depth increases, surface subsidence decreases. From the analysis of results, the following was found.

• Low values of fault friction significantly influence ground surface subsidence since blocks around the excavation can slide due to the concentration of the stress and plastic deformation around the excavation area.

• The excavation sequence and direction of the mining at the Máza site had a slight impact on the magnitude of surface subsidence.

• Small surface subsidence was observed when overburden rock formations were characterized by a high Poisson's ratio. On the other hand, maximum ground surface subsidence was found to occur when Young's modulus of the overburden rocks was low.

• Maximum shear displacement occurred at faults with a low friction angle.

Result outputs of the parametric study of ground subsidence and fault reactivation were used for the development of the web-based and interactive Environmental Hazards and Risk Management Toolkit (EHRM) of in situ coal conversion (Tranter et al., 2022).