Effects of CO2 injection factors on preventing spontaneous coal combustion in coal mine gobs: a simulation approach

The spontaneous combustion of residual coal in the gob seriously threatens the safety of coal mining. Injecting CO2 into the gob not only prevents the residual coal from spontaneous combustion but also realizes CO2 storage in the mined areas. Injection flux and burial depth of the port are crucial for CO2-preventing fire in coal mine gobs. In this study, the distribution of the oxidation zone in the Tanyaoping coal mine was field-measured, and the coal oxidation kinetic model was built by the adiabatic test. Then, a 3-D mathematical model was constructed based on the conditions of the 5011 working face by COMSOL Multiphysics. Furthermore, the coupled effects of the two factors on the distribution of the oxidation zone were investigated. Increases in both injection flux and burial depth result in a decrease in the oxidation zone volume. The reasonable ranges of the injection flux and burial depth are 540–720 m3 h-1 and 30–40 m, respectively. These results provide some guidelines on how to prevent the spontaneous combustion of residual coal in mine gobs.


Introduction
China is the largest producer of coal in the world. Coal mining normally faces several types of accidents that can cause environmental disruption, fatalities, and equipment loss Li et al., 2022). The reaction between coal and O 2 at low temperatures usually results in spontaneous combustion and uncontrolled coal fires (Zhang et al., 2023). During underground coal mining, spontaneous combustion of residual coal in the mine gob often leads to a wide range of secondary disasters, such as methane explosions, carbon monoxide generation, and ventilation disorders. The oxidation process of residual coal in the gob is extremely complex. High-mechanized technology, especially coal caving mining, results in considerable coal resources left in the gob. However, it is necessary to keep enough fresh air through the airways and working face for human survival and methane dilution. Therefore, the residual coal in the gob is exposed to continuous O 2 and has the possibility of spontaneous combustion. The hazardous areas of self-ignition in the gob depend on many factors, such as air leakage, thickness of residual coal, and the advanced rate of the working face (Zhu and Wen, 2023). Normally, the "three zones" (the cooling zone, oxidation zone, and suffocation zone) in the gob are determined based on the distribution of O 2 concentration in the gob area (Ma et al., 2020).
Spontaneous combustion of coal is influenced by coal oxidation reactions, heat accumulation, and the characteristics of porous media (Si et al., 2021). Therefore, many technical solutions were presented to prevent coal from spontaneous combustion and inhibit coal fire, which normally focus on two aspects: isolating coal active sites from O 2 and cooling the temperature (Dou et al., 2022). The present strategies can be broadly classified as uniform pressure ventilation (Lu et al., 2017), grouting (Zhang et al., 2018), spray inhibitors (Lu et al., 2020;Sun et al., 2022), pumping inert gas (Zhu and Wen, 2023), and injecting gel foam or gel (Lu et al., 2018;Lu et al., 2021;Han et al., 2022;Huang et al., 2022). CO 2 is an inert gas and can dilute the O 2 concentration and protect the coal from heat release. Furthermore, broken coal can absorb CO 2 and thus saved from oxidation (Guan et al., 2018). With the advancement of the working face, the volume of the mine gob increases. The resistance of internal airflow is small, and the porosity is large. Thus, the coal mine gob is an ideal location for storing CO 2 . Injecting CO 2 into the gob can achieve the dual benefits of fire prevention and CO 2 storage.
This study is based on the actual conditions of the 5011 working face of the Tanyaoping coal mine. To achieve the fire prevention of residual coal in the gob by injecting CO 2 , the oxidation zone was calculated by mathematical simulation using COMSOL Multiphysics. This simulated result without CO 2 injection was verified to be reasonable by comparing it with the field measurement. By studying the effect of the injection flux and burial depth of the injection port on the oxidation zone, the superior injection parameters were determined.

Overview of the 5011 working face
The Tanyaoping coal mine is located in Lvliang city, China. The recoverable coal seams consist of 5# and 10# coal seams. The 5011 working face is in coal seam 5# and adopts the "U" ventilation method. Coal seam 5# is the main production seam and is around 4.4 m thick. The width and design length are 180 and 2,000 m, respectively. The width and height of the intake and return airways are 5 and 4.5 m, respectively. The thickness of the residual coal that remains in the mined area, that is, the gob, is about 0.4 m. The mined coal has a spontaneous combustion tendency.
As shown in Figure 1, two sensors were fixed in the intake and return airways to collect the gas. The monitoring system was based on the measurement of bundle tubes. With the advancement of the working face, O 2 concentrations at different buried depths were drawn out through the bundle tubes using a vacuum pump. A gas chromatography analysis system was used to analyze the O 2 concentrations of the collected bags. Figure 2 shows the measured O 2 concentration values in the intake and return airways. The O 2 concentrations at the observation points decrease with increasing buried depth. According to the actual  condition and other studies, the "three zones" were divided by using an O 2 concentration of 10%-18% (Ma et al., 2020). The O 2 concentration at the point in the intake airway decreased to 18% and 10% at 90 and 147 m, respectively. For the measured point in the return airway, the buried depths of O 2 concentrations of 18% and 10% are 30 and 73 m, respectively. Normally, there are three parameters to determine the "three zones" in the gob: the air leakage velocity, O 2 concentration, and heating rate (Deng et al., 2018). However, the O 2 concentration is the most accessible and was, therefore, applied in this study to judge the oxidation zone.

Coal oxidation kinetic model
Coal exothermic processes under low temperatures include three forms: physical adsorption, chemical adsorption, and chemical reactions (Liang et al., 2021). Chemical adsorption and chemical reactions are recognized as the main reasons for promoting the self-heating of coal under low temperatures (Zhang et al., 2019;Qu et al., 2023). Although partially active functional groups are selfreacting, a significant amount of heat will be released from the coal-O 2 reaction. Therefore, a linear relationship between the O 2 consumption rate and heat generation rate is generally considered to exist, expressed as follows: where q is the heat generation rate of the coal reaction (W·m −3 ); ΔH is the coal oxidation heat (J·mol −1 ); and r is the O 2 consumption rate (mol·m −3 ·s −1 ). Based on the coal oxidation dynamics theory, another major influence on coal oxidation is temperature, which increases, causing the activation of functional groups, i.e., the exothermic intensity increases with increasing temperature. Most of the models on O 2 consumption in the previous studies were based on the Arrhenius equation, expressed as follows: where A is the pre-exponential factor (s − ); E is the activation energy (J·mol −1 ); R is the gas constant (8.314 J mol −1 ·K −1 ); c o is the O 2 concentration (%); and T is the temperature (K). The kinetic parameters in Equation 2 are normally obtained from the adiabatic oxidation test. As shown in Figure 3, the coal sample was mined from the 5011 working face of the Tanyaoping coal mine. Moisture, ash, volatile matter, and fixed carbon of coal are 7.64, 2.21, 26.24, and    (Wang et al., 2010). Once the temperature was stable, the mixed gases with 21% O 2 were injected into the reaction vessel at 10 mL min −1 . It should be ensured that the oven temperature is synchronized with the coal temperature using the temperature controller. The change in coal temperature over time was recorded using the recorder for later analysis. Figure 4 shows the temperature profile versus time of the adiabatic oxidation test. Under the adiabatic conditions, the coal sample was subjected to accelerated self-oxidation so that the heat released in the last step was completely converted to internal energy, which can be described by Equation 3 as follows: where ρ c is the density of coal (kg·m −3 ) and C pc is the specific heat capacity of coal (J·kg −1 ·K −1 ). Equation 4 can be rearranged by taking the natural logarithm, which is expressed as follows: The relationship between lndT/dt and −1/T shown in Figure 5 is a nonlinear positive correlation, which does not follow the standard Arrhenius equation, and this phenomenon also appears in other studies (Zhang et al., 2020;Yoruk and Arisoy, 2022). It is difficult to fit the curve with one or two straight lines. The curve was divided into four stages with three break points of 48.1, 56, and 68.7°C. In each stage, lndT/dt and −1/T exhibit good linearity, as shown in Equation 1. The slope of the fitting line in each stage presents the value of -E/R, and the intercept is the value of lnΔHAc 0 /ρ c C pc . The apparent activation energies of the four stages are 143. 25, 57.45, 32.89, and 19.48 kJ mol −1 , respectively. Accordingly, the preexponential factors of the coal sample are 2.53×10 20 , 2.83×10 6 , 358.31, and 3.2 s −1 , respectively. In the numerical calculation, different E and A are invoked based on the coal temperature for the heat generation model.

Numerical simulation
The process of spontaneous combustion of residual coal in the gob area will involve heat transport and generation, O 2 consumption and transport, CO 2 emission and migration, and gas airflow. As shown in Figure 6, all the involved equations are interrelated based on the fluid-solid-thermal coupling (Xia et al., 2016;Zhang et al., 2019;Qi et al., 2021). The continuity equation and the momentum conservation equation are as follows: where ρ g is the density of air (kg·m −3 ); ε is the porosity of gob; t is the time (s); U is the seepage air velocity vector (m·s 1 ); p is the static pressure (Pa); τ is the stress tensor (Pa); μ is the air viscosity (Pa s); K is the permeability (m 2 ); and F is the gravitational body force (N·m −3 ). The relationship between porosity and permeability is expressed as where d is the particle size (mm). With the working face advanced, the stress status of overlying strata is damaged, and the roof-bed formation collapses. The overlying strata are vertically divided into a caving zone, fissure zone, and bent deformation zone based on the differences in the balance status and re-compaction degree of the damaged strata (Qin et al., 2016). Therefore, it is a significant challenge to build mathematical models with reasonable permeability and porosity in the gob area. According to the observation of the mine pressure, the compaction bulking factor, k, of the gob follows the law (Xia et al., 2015): where k p,min and k p,max are the initial bulking factors before and after compacting in the gob area, respectively, (1.15 and 1.5). a 0 and a 1 are the decay ratios of the bulking factor in the dip and strike directions of the working face, respectively, (0.268 and 0.0368 m −1 ). d 0 and d 1 are the distances between the point (x,y) in the gob and the working face and coal pillar, respectively. ξ 1 is the adjustment factor controlling the distribution pattern of the model, 0.233. The local thermal non-equilibrium is considered in the energy equation of conservation, which is given by For gas, ερ g C pg zT g zt + ρ g C pg U∇T g ε∇ λ g ∇T g + q gc , where C pg is the specific heat capacity of air (J·kg −1 ·K −1 ); λ g and λ c are the thermal conductivities of air and coal, respectively (W·m −1 ·K −1 ); T g and T c are the temperatures of air and coal, respectively (K); q o is the heat generation rate of coal (W·m −3 ) obtained using Equation 3; q gc is the solid-gas heat exchange around coal particles (W·m −3 ), which has q gc h gc a T c − T g , where h gc is the convective heat transfer coefficient (W·m −2 ·K −1 ), indicating the heat exchange between coal and pore air; a is the specific surface area, that is, the exposed area of coal particles per unit volume (m -1 ). O 2 conservation is a unique consideration without CO 2 injection. However, we must consider the conservation of O 2 and CO 2 for fire prevention by injecting CO 2 into the gob, which is expressed as follows: ε zc m zt where c m is the concentration of CO 2 (mol·m -3 ); D o and D m are the effective diffusivities of O 2 and CO 2 , respectively (m 2 ·s −1 ); and r o is the O 2 consumption rate (mol·m −3 ·s −1 ). r m is the CO 2 generation rate (mol·m −3 ·s −1 ) of coal oxidation, whose ratio to r o is 0.1 (Taraba et al., 2014). The oxidation reaction of residual coal in the gob is affected by working face advance. Based on Equation 2, r o can be expressed as where w is the influence coefficient of working face advance.
In addition, to ensure the normal operation of modeling, the following main assumptions are supplemented (Ma et al., 2020;Zhang et al., 2022): (i) The effects of moisture transport and evaporation are not considered. (ii) Thermal conductivity and heat capacity of coal and gas are independent of other parameters. (iii) The transfer of heat in the gob merely involves convection and conduction, and the radiative transfer is not taken into account. (iv) The gas in the gob is incompressible.
According to the actual conditions of the 5011 working face, the corresponding geometric model has been established. COMSOL Multiphysics software was applied to simulate the oxidation process of residual coal in the mine gob. The whole model consists of four parts: the mining area, the inlet airway, the return airway, and the working face. The cross-sectional area of both the inlet and return airways is 22.5 m 2 , and the cross-sectional area of the working face is 25 m 2 . The length of the inlet and return airways is 5 m. The length of the working face is set as 220 m, and the length and height of the gob area are 180 and 15 m, respectively. For the condition without CO 2 injection, grid independence was checked by carrying out three sets of grids. Eventually, after considering the CPU time and the sensitive parameter, 34,954 mesh blocks are divided among this model, including 32,389 mesh blocks in the gob area. The mesh blocks of the models for CO 2 injection were slightly larger than those without CO 2 injection. The time step was set as 1 h. As shown in Figure 7, the CO 2 injection port was a round tube with a radius of 0.2 m and a length of 0.2 m. Two parameters, injection flux (U, m 3 h −-1 ) and injection port burial depth (L, m), were considered to investigate the performance of CO 2 injection in preventing residual coal fires. The parameters used in the simulations are listed in Table 1. The boundary and initial conditions used in the mathematical model are illustrated in Table 2.

Analysis of the oxidation zone without CO 2 injection
As shown in Figures 8A, B, the porosities and permeabilities of the gob near the working face and pillars are larger than those of the  inner zone due to the mining process. Fresh air enters the working face from the inlet airway with a velocity of 0.7 m s −1 , and more than 90% of the ventilation flux passes through the working face and then flows into the return airway. The airflow streamlines in the gob area, normally from the windward zone near the inlet airway to the zone near the return airway. The airflow velocity near the working face is larger than that in the inner zone due to the smaller forcedconvection resistance.
With the deepening of the gob, the O 2 concentration gradually decreases. As shown in F8(e), on the side of the inlet airway, the fresh air from the intake airway enters the gob and flows toward the deep. However, gases with low O 2 concentrations flow from the deep to the return airway. Therefore, the O 2 concentration near the side of the return airway shows a more rapid decline than that of the intake airway.
This distribution of the oxidation zone has been observed in many other field measurements and simulations (Xia et al., 2015;Chu et al., 2019;Xu et al., 2020). Coal near the working face undergoes a strong oxidation reaction because of the abundant O 2 , accompanied by the release of heat. We can find that the high-temperature zone is located on the windward side of the oxidation zone, as shown in Figure 8C. It is due to the fact that O 2 was consumed in the high-temperature zone, causing fewer O 2 molecules to enter the interior of the gob. Therefore, it is difficult to transport enough O 2 to the inner zone, resulting in a low concentration of O 2 within the area. The distribution of CO 2 in the gob is the opposite of the distribution of O 2 because most of the CO 2 in the gob is derived from low-temperature oxidation of residual coal, as shown in Figure 8D. The overall CO 2 values are much lower than those of O 2 . The CO 2 concentration near the working face is lower than that in the inner zone due to the dilution of airflow seepage. Figure 8E is a three-dimensional distribution of O 2 concentration in the gob area. The O 2 concentration in the upper area in the same vertical direction is slightly higher than that in the lower area. The volume of the oxidation zone without CO 2 injection is calculated to be 121,000 m 3 . The average width of the oxidation band is 45 m. Figure 8F shows the range comparison of the oxidation zones (z=0.2 m) between the field measurement and simulation. The two results show good alignment. It indicates that the numerical simulation can reflect the performance of the self-ignition of the residual coal in the gob.

Analysis of the base case with CO 2 injection
The base case of U = 540 m 3 h −1 and L = 20 m was calculated, and the distributions of O 2 , CO 2 , and temperature in the gob area are shown in Figure 9. CO 2 injection has the properties of a displacement effect on the air leakage, a dilution effect on the O 2 distribution, and a cooling effect on the high-temperature zone. Therefore, the spontaneous combustion of residual coal is delayed or prevented. CO 2 injection can dilute the O 2 concentration. CO 2 concentration is much higher than O 2 concentration downstream of the injection port. The oxidation zone on the inlet airway side moved significantly toward the injection port, and the width of the oxidation zone decreased significantly. CO 2 was continually mixed with leaked O 2 , and the dilution effect on the O 2 distribution was gradually attenuated. The high-temperature zone under CO 2 injection was a smaller area than that without CO 2 injection. CO 2 also suppresses the maximum temperature in the gob. It indicates that CO 2 can restrict the self-ignition of coal. In addition, the area of the oxidation zone for CO 2 injection on the surface z = 0.2 m is 1300 m 2 , much smaller than that without CO 2

Types
Positions Values Frontiers in Earth Science frontiersin.org injection, which is 7400 m 2 . However, the volumes of the oxidation zone with and without CO 2 injection are 63,000 m 3 and 121,000 m 3 , respectively. The injection port is located near the bottom, causing a large amount of CO 2 cluster in the lower area and, therefore, not easily diluting O 2 in the upper area. As a result, the width of the oxidation zone in the upper area is much greater than that of the lower area.

Effect of the CO 2 injection flux
Based on the base case and the corresponding simulation without CO 2 injection, the effects of the CO 2 injection flux on the prevention of residual coal fires were analyzed. The assumed CO 2 injection fluxes are 180, 360, 540, 720, 900, and 1080 m 3 h −1 , respectively. Figure 10 shows the oxidation zone in the gob for different CO 2 injection fluxes at L=20 m. An increase in the injection flux weakens the O 2 concentration in the gob and increases the displacement effect. When U ≥ 180 m 3 h −1 , the higher pressure of the intake airway and gravity suppress the injected CO 2 , which flows along the lower side of the gob in the direction of airflow seepage, as shown in Figure 11. At U= 180 and 540 m 3 h −1 , CO 2 concentrations in the upper zone of the gob are much lower than those in the bottom zone. Because the CO 2 injection fluxes (≤1080 m 3 h −1 ) are smaller than the air flux from the intake airway (1050 m 3 min −1 ). CO 2 injection hardly affects the air leakage streamlines, and thus, the high CO 2 areas at different U values are extremely similar. With the increase in CO 2 injection flux, the overall distribution of the oxidation zone moves toward the working face. The oxidation zone in the lower gob area moves toward the upwind side of the injection port when U ≥ 0. As wind speeds increase, the distance between the iso-surfaces increases by 10% to 18%, reducing the volume of the oxidation zone. The volumes of the oxidation zone at U = 180, 360, 540, 720, 900, and 1080 m 3 h −1 are 89,000, 62,500, 43,000, 2,950, 24,000, and 23,000 m 3 at L=20 m, respectively. Therefore, when U < 900 m 3 h −1 , the volume of the oxidation zone decreases with the injection flux, and when U > 900 m 3 h −1 , the volume changes little as flux increases. This phenomenon also appears in other cases with different L values.  Figures 12, 13 show O 2 and CO 2 distributions in the gob with different burial depths at U = 540 m 3 h −1 . Once CO 2 is injected, the oxidation zone on the inlet side is significantly reduced. This is because the high CO 2 concentration near the injection port effectively dilutes O 2 . The closer the inlet port is to the working face, the closer the oxidation zone on the inlet side is to the working face, and the weaker the effect of the O 2 concentration in the upper area of the inlet side being displaced by CO 2 . It is due to the high leakage strength near the working face and a large amount of fresh air that suppresses CO 2 diffusion by the dilution effect. Therefore, as shown in Figure 13, the smaller L is, the lower the CO 2 concentration within the gob. At L = 10, 20, and 30m, the CO 2 concentration in the upper area of the gob is low. Moreover, as the burial depth of the injection port increases, the airflow leakage from the working face is small, and CO 2 can migrate to the upper area without high resistance. The volumes of the oxidation zone at L= 10, 20, 30, 40, 50, and 60 are 65,000, 43,000, 27,000, 23,000, 19,200, and

FIGURE 11
Simulated CO 2 distribution with different injection fluxes at L = 20 m.

FIGURE 12
Simulated oxidation zone with different burial depths at U = 540 m 3 h −1 .

FIGURE 13
Simulated CO 2 distribution with different burial depths at U = 540 m 3 h −1 . 18,500 m 3 at Q=540 m 3 h −1 , respectively. As L increases, the width of the oxidation zone in the upper area decreases rapidly. The volumes of the oxidation zone are almost identical at L = 50 and 60 m. In addition, the change in the burial depth of the injection port has little effect on the distribution of the oxidation zone on the return side.
5.5 Comprehensive analysis of coupling the CO 2 injection flux and location Figure 14 shows the effects of coupling factors, that is, the CO 2 injection flux and location, on the predicted volume of the oxidation zone. A decrease in the CO 2 injection volume or burial depth of the port will reduce the volume of the oxidation zone. When U < 720 m 3 h −1 , an increase in the burial depth significantly reduces the volume of the oxidation zone, especially at a lower injection flux. When U=1080 m 3 h −1 , the change in the position of the injection port does not have a significant effect on the distribution of the oxidation zone. At L = 10 m, the volume of the oxidation zone decreases linearly with the injection flux. However, when L > 30 m, the volume dramatically reduces with increasing U from 0 to 540 m 3 h −1 , but changes little as flux continues increasing. In practice, longer gas injection lines are avoided, which can make operations difficult. Smaller U values can save costs and reduce CO 2 concentrations in the upper corner and working face when achieving similar effects of suppressing coal spontaneous combustion for different U values. Therefore, based on the simulation results, the parameters of U = 540-720 m 3 h −1 and L = 30-40 m are recommended for the 5011 working face.

Conclusion
To accurately predict the effectiveness of CO 2 injection on spontaneous combustion in the coal mine gob. A mathematical simulation by COMSOL Multiphysics was applied in this study.
According to the conditions of the 5011 working face in the Tanyaoping coal mine, a 3-D mathematical model considering fluid-solid-thermal inter-relationship was established. The coal oxidation kinetic model was derived from the adiabatic oxidation test, where four stages were considered in the model. The apparent activation energies of the four stages are 143.25, 57.45, 32.89, and 19.48 kJ mol −1 , respectively. A good agreement between the field measurement and numerical simulation indicates the rationality of the mathematical model. CO 2 injection can dilute the O 2 concentration and prevent spontaneous combustion of coal. The coupling effects of injection flux, U, and the burial depth of the injection port, L, on the oxidation zone were investigated. The oxidation zone volume decreased with increasing U and L. When U or L is small, the upper width of the oxidation zone is wider, and as U or L increases, the upper width gradually decreases until it is no longer sensitive to parameter changes. An optimum injection plan is to keep the volume of the oxidation zone at a lower value in consideration of economic viability and operational parameters. The injection parameters of U = 540-720 m 3 h −1 and L = 30-40 m are recommended after considering the activity of the 5100 working face. The volume of the oxidation zone was estimated to be 18,500-23000 m 3 .

Data availability statement
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.