Numerical simulation of coalbed methane under the mining conditions based on dynamic virtual well and variable permeability field

After coal seam mining,stress change and rock movement in adjacent rock strata, make permeability and other fields change in adjacent coal seam. At the same time, due to the exploitation of coal seams, material exchange may occur in the adjacent coal seams. It also caused desorption of near coal seam near the mining working face. Recent years, with the development of coalbed methane mining technology under pressure relief in mining areas, it is necessary to carry out numerical simulation under coalbed methane mining conditions to realize the optimization of coalbed methane extraction simulation under mining conditions. The coupling effect of dynamic virtual well and variable permeability field advancing with the mining face is used to simulate the output of coalbed methane under mining conditions: the dynamic virtual well is used to advance with the mining face to characterize the influence of the working face on the material exchange of adjacent coal seams in the mining process; the use of variable permeability field to realize the influence of the permeability of the adjacent coal seam during the advancement of the working face. Finally, through the dynamic coupling of variable permeability field and dynamic virtual well with the mining face, the rapid pressure relief effect of high speed and short time on the adjacent layer under the influence of mining is characterized. The simulation results of the computer simulation software written by the algorithm show that when the mining effect is not considered, that is, when the algorithm does not consider the dynamic virtual well production and the variable permeability field, it is basically consistent with the commercial software simulation results, but for the consideration of mining condition is a situation that commercial software cannot simulate, and the simulation results of this algorithm are more consistent with the actual monitoring data of the mine. This method achieves numerical simulation of CBM under mining conditions by improving the mature algorithm, which can provide effective guidance for future CBM capacity prediction under mining conditions.


Introduction
Coal-bed methane reservoir is a kind of unconventional gas reservoir, which is very different from conventional gas reservoir in terms of storage method and extraction mechanism. Coalbed methane has the characteristics of autogenous self-storage formation, and its main component is methane, which is stored in the coal matrix system with low permeability in the form of adsorption. In addition to matrix systems, cuttings systems are widely developed in coal seams, which are the main seepage channels of CBM reservoirs. Coal seams are generally characterized by low porosity, low permeability and strong non-homogeneity. Coalbed methane is an effective alternative to conventional energy and is widely distributed worldwide. The three countries with the largest coalbed methane reserves in the world are Russia, Canada and China. (Qun 2003;Qin 2006;Song et al. 2013). The current coalbed methane extraction technology is changing from the old way of underground extraction to a combination of surface development and underground extraction. The combination of surface development and underground extraction is usually accompanied by the mining of coal seam workings, which can be divided into two categories: protected seam mining and unloading coalbed methane surface well extraction. Surface direct well extraction of CBM from the extraction area has been successfully used in the United States, Australia, Ukraine and other countries, and has become a common method for developing CBM (Li et al. 2009;Song et al. 2005;Zhang and Cheng 2015;. In the 1990s, as the mechanism of coal-bed methane transport and storage continued to be studied in depth, the non-equilibrium adsorption model was generally accepted. This model considers that coal-bed methane is mainly adsorbed on the inner surface of the matrix and obeys Darcy seepage in the cut fracture, which is currently represented by the Comet model and the Coalgas model. In 2001, Scott and Larry (2001) proposed a three-pore dual-permeability model, which is more realistic than the dual-porosity single-permeability model. Meanwhile ARI  has launched its latest product, COMET3, the most advanced software in all COMET series simulation software. The software is a model of triple pore, double permeability and multiple gas co-adsorption, which is widely used in coalbed methane exploration and development. There are many studies on the numerical simulation of CBM development without considering the mining motion, and there are mature commercial software for the numerical simulation of conventional CBM (Peng et al. 2015;Guo et al. 2012;Chen et al. 2004;Philpot et al. 2013). Since COMET3 and COALGAS are developed based on the characteristics of coal reservoirs under simple geological conditions (Ding et al. 2014;Parra et al. 2006;Han and Sang, 2009), they are not applicable to coal reservoirs under complex geological conditions with strong non-homogeneity and large anisotropy. While the traditional commercial reservoir numerical simulation software such as Eclipse, CMG, tNavigator, etc. can also simulate the development of CBM, which can realize the simulation from single medium, double medium (double pore single percolation, double pore double percolation), triple medium (triple pore single percolation, triple pore double percolation) and single component, double component to multi-component, however, there is still no detailed consideration of the mining conditions and no modules or keywords are set for the mining conditions.
As a reservoir full of natural fractures, the permeability of coal reservoirs is also determined by many fracture characteristics, such as segregation density, length, and directional distribution. Gray (1987) proposed the first permeability model for coal reservoirs, taking into account compaction and matrix shrinkage effects, and suggested that stress variation is exponentially related to permeability. Harpalani and Chen (1995) proposed a model based on matchstick theory that relates permeability changes to matrix dimensions. The model was enriched by Ma et al. (2011). Levine (1996) proposed a model to calculate the variation of fracture opening and correlated permeability with fracture opening. Robertson and Christiansen (2006) gave a model for permeability under variable stress conditions. Although CBM scholars have presented insights into the complex permeability change patterns of CBM, no relevant reports on strain changes and permeability change patterns in CBM reservoirs under mining conditions have been published. CBM simulation model is built on the basis of reservoir permeability variation with stress and gas desorption equation. Luo (2003) proposed a three-dimensional model of coalbed methane reservoirs and discussed the method of solving the model in full implicit association. In the same year, Qun (2003) gave a mathematical model as well as a numerical model of gas-water two-phase flow in three-dimensional, dual-pore, non-equilibrium adsorption, and proposed steady-state conditions in coal-bed methane reservoirs, and developed the corresponding numerical simulation software for coal-bed methane reservoirs. Zhang and Tong (2008) established a numerical simulation model of three-dimensional, non-equilibrium adsorption, coupled gas and water two-phase flow in dense coalbed methane reservoirs under the proposed steady-state conditions considering the initiation pressure gradient, and gave a fully implicit finite-difference format and numerical solution method for the model. Meng et al. (2018) established a multi-factor coupled mathematical model for the dynamic change of coal reservoir permeability by considering the effects of stress sensitivity, matrix shrinkage effect, and slip effect. Sun et al. (2018) developed a multi-factor coupled gas-water phase effective permeability model for CBM production process based on the P-M model and the Kerry relative permeability model, which quantified the effective permeability of the gas/water phase as a function of pressure and saturation. Peng et al. (2019) studied in detail the effect of coal dust on the early drainage and development permeability of coalbed methane reservoirs, and established a mathematical model of the four phases of gas-water-fluidsolid-static describing the fluid transport pattern in coal seams, which was solved using finite-difference equations. Danesh et al. (2022) studied the dynamic stress distribution in CBM reservoir development and established a mathematical model of CBM with coupled fault stresses and fault interactions, and the effect of matrix shrinkage on reservoir stress generation was considered in their paper. In summary, although many scholars have made outstanding contributions to the numerical simulation of CBM development, no one has conducted a detailed study on the numerical simulation of unloading and extraction development of CBM.
The current commercial software does not have the capability to do specific calculations for this important CBM development method, and no one has yet developed a numerical simulation method that takes into account the mining conditions. Considering the importance of surface well extraction technology for CBM development, there is an urgent need to develop and solve numerical simulation methods for CBM extraction considering the influence of extraction to adapt to different development geological conditions. For the numerical simulation of coalbed methane in surface extraction, especially under the conditions of mining influence, the influence on the permeability field, seepage field and adsorption and desorption of the adjacent coal seam during the advance of the working face during the mining process must be considered. In this paper, a mathematical model for recovering CBM with dynamic virtual wells and dynamic variable permeability coupled with multi-physical fields of working face advancement is developed, which seeks to simplify the actual model, and the model was also solved and the productivity, pressure distribution and gas content saturation field of the mining development area were calculated and compared with the actual mine production data. At the same time, we borrow the current mature model to accurately reflect the main mechanism, and realize the numerical simulation analysis and optimization of coalbed methane development under mining conditions through numerical simulation calculation.

Physical model
Mining operation includes dynamic changes of this layer and adjacent layer caused by tunnelling roadway and coal face. This study mainly focuses on the surface extraction of the depressurized coalbed methane in the mining area. At this time, the mining affected area can be divided into the mining affected area of the coal seam and the adjacent mining affected area. As for the mining layer, with the advance of the mining face, the coal seam in the heading roadway is extracted, and the regional stress field near the working face changes during mining, which changes the permeability field and pressure field. Based on the analysis of the field data, the whole grid system is set unchanged when the model is established, and the flow boundary of the dynamic virtual well is used to replace the actual mining roadway. The energy and material exchange between the mining roadway and adjacent layers in the mining process is characterized by the production of the dynamic virtual well, so as to simulate the pressure relief effect of the mining work on the adjacent layers. Determining the coefficient of permeability change based on the stress field change during mining, coal seam quality jointly, and thus characterizing the permeability as a function of the distance from the mining workings. Eventually, the dynamic change of mining workings advance position with time is expressed as the effect of the mining process on the permeability field of the neighbouring layers, and the mining conditions are finally simplified as the coupling effect of dynamic virtual wells and dynamic permeability. As shown in Figure 1. The basic hypotheses are: 1. There are two-phase fluid (gas and water) in the coal reservoir and two phase fluid in porous media, namely Newton fluid, and the flow is laminar flow movement, in agreement with Darcy's law; 2. The diffusion process of gas in coal base is non-equilibrium quasi steady process in agreement with the first law of Fick; 3. PVT characteristics of the two-phase fluid are in accordance with the law of micro compressibility, without considering the gas dissolved in water; 4. In the process of mining working face, the permeability field of adjacent layers is affected by the mining working face, and the mining may cause pressure relief of adjacent layers and desorption of coalbed methane.

Model establishment
The whole mathematical model consists of three elements: control equations, diffusion equation of gas from matrix to fissure and auxiliary equations. Control equations. Based on the simplification of the previous model, the dynamic virtual well production is used to characterize the effect of pressure relief (water production) on the protective layer in the process of mining face advance; The dynamic permeability field is used to characterize the influence of mining conditions on the permeability field change in the process of working face advance. Dynamic permeability is determined as a function of the mining distance from the mining face based on the seam stress conditions and the seam quality, while this function moves with the advance of the mining face. Finally, the basic mathematical model of the following fractures can be established: In this formula, c s is the permeability adjustment coefficient in fractures. For the mining affected area according to the distance from the mining working face set as a function of changes greater than 1, for the mining unaffected area set to 1. where q p is the mass gas source and represents the gas diffusing from the matrix into the fissure, and only methane is present in the CBM reservoir. q s denotes the dynamic virtual well volume flow rate of the mining operation, and the location of the virtual well changes position over time together with the advancement of the mining workings.
The analysis of equation 1 shows that if the permeability adjustment coefficient in the whole model is set to 1 and the flow rate of the dynamic virtual well is set to 0, the model can be changed into the traditional CBM seepage model.
Diffusion equation of gas from matrix to fissure. It is generally accepted that only single-phase gas is diffused into the fracture of coal matrix block at pseudo-steady state because water cannot enter into the micro pores of matrix block (assuming that only methane exists in the coal matrix), taking matrix gas desorption and diffusion processes as diffusion of the pseudo-steady state. Pseudo-steady diffusion is in agreement with the first diffusion law of Fick.
According to the Langmuir equation, the coal bed gas mass that concentrate in per cubic metre of coal is V L P m P L +P m , this is the concentration of adsorbed gas, which is: The total concentration c px is based on the total volume of the matrix: The pseudo-steady diffusion process is based on the first law of Fick, and the time change rate of the total concentration c px in the matrix block is believed to be directly proportional to the difference c pxc p, which is in accordance with the first diffusion law of Fick: Auxiliary equations. In addition to the representation of the fluids (coal seam gas and water) migration process with partial differential equations according to conservation of mass, further auxiliary equations to support this model is necessary, including saturation equation and capillary force equation: In numerical simulation, well is treated as the special boundary conditions. If there is a well in the grid (i, j, k), volume flow rate is Q v , it can be directly substituted into the seepage equation as the source and sink terms. Production wells are negative and injection wells are positive. If wells produce at a given flow pressure, it is necessary to use grid pressure P i,j,k , and flow pressure P wf to represent Q v , usually regarded as pseudo-steady state flow, and plane radial flow formula of producing well or injection well of l phase fluid meet: Where: ln r e r w + S K = K fx K fy , r e = 0.28 [(K fx / K fy ) 1/2 Δy 2 + (K fy / K fx ) 1/2 Δx 2 ] 1/2 (K fx / K fy ) 1/4 + (K fy / K fx ) 1/4 In the process of numerical simulation, the outer boundary condition is another boundary condition, which is usually divided into two kinds of situations, the constant pressure boundary and the closed boundary (this paper uses the closed boundary). Expressed as the following mathematical expression: Boundary condition can also be written as: The formula expresses the normal derivative of pressure on Γ 2 boundary, with the numerical identity as zero, indicating that boundary and adjacent grid pressure is equal, or there is no the exchange of material, namely a closed boundary.
Mining roadway can be treated as variable permeability strips, which mainly changes permeability in the fracture. (In a given period of time, it is considered as a constant, and the roadway advancement is simulated by the change in permeability over time. When the composition of the roadway is connected to the well, the outer boundary of the roadway is treated as an equivalent injection well.).
K f (x, y, z, t)| Γ 1 = c(x, y, z, t) × K, t = t 1 , t 2 . . ., Known function. it is assumed that the roadway is advancing along a certain direction, and the distance and velocity are known .

Simulation results and field applications
In order to verify the accuracy of the proposed method, this paper takes Dingji coal mine as an example. Firstly, the simulation of CBM development without considering the mining conditions was carried out, and the results were compared with the conventional Comet software without considering the mining conditions to verify the accuracy of the method proposed in this paper. Then, a virtual well was added and the permeability field was changed according to the stress variation calculation to perform the numerical simulation considering the mining conditions, while other parameters of Dingji coal mine remained unchanged.

No consideration for mining
In accordance with the method employed in this paper, the numerical simulation is carried out by using Fortran programming language, and dingji coal mine in Huainan, Anhui province is chosen as the theoretical example. The mine has low permeable coal seam, and the key parameters are as follows: Dingji Coal Mine working face is in −826 horizontal west mining area, with the length of oblique of 208m, and its trend length is 1780m. The 11-2 coal seam is mainly mined with an average thickness of 2.45m. Direct top is sandy mudstone with an average thickness of 0.8 m, old top is fine sandstone with an average thickness of 10.5 m. The protected coal seam is the upper coal seam13-1; average thickness is 3.1 m; gas content is 5.25 m 3 /t; distance is about 80m. The permeability of coal bed is 1 × 10 −3 μm 2 , the original formation pressure is 3.77MPa; porosity is 2%.
The above data are used to calculate the basic parameters through computer simulation calculation. In order to verify the reliability of the algorithm, a basic simulation without considering mining conditions was first carried out. The simulation range is a 2,700m × 2,700m rectangular area with a zone grid of 9 × 9 × 1 and nine Wells evenly deployed. In the commercial numerical simulation software such as Eclipse, CMG, tNavigator and Comet, the mining conditions in the CBM development process were often not considered, so that the predicted production data could not be well matched with the actual production data. When the mining conditions are not considered, the time comparison curves of CBM production and pressure calculated based on this algorithm and commercial software comet are shown in Figure 2 and Figure 3.
From Figure 2, we can see that the peak daily gas production of the block in the method applied in this paper is high at the beginning and decreases after a period of high production, gradually lower than the daily production simulated by Comet, and lower than Comet in the subsequent period; and from Figure 3, it can be seen that the average reservoir pressure of the blocks in the calculation method in this paper decreases at a lower rate than the average pressure calculated by Comet, and after about 1000 days the average pressure of the calculation method described in the text is slightly lower than that simulated by the commercial software Comet. As a whole, the simulation results of the equation solving method used in this paper have the same trend as the adaptive fully implicit calculation method used by the commercial simulation software Comet and others, and the computational results do not differ much from each other. Therefore, it can be considered that the algorithm in this paper is basically the same as the comet simulation results, which indicates that the present algorithm is reliable.

Consider mining
The simulation when mining is considered is still performed in Dingji coal mine, and the initial parameters such as formation temperature and pressure, porosity permeability are the same as Figure 2. Comparison of daily gas production simulation results by the commercial software with this paper method without mining being considered.  when mining is not considered. When the introduction of the mining conditions was carried out, the mining face were advanced at a speed of 5 m/d (the same speed as the dynamic virtual well advance). In addition, the dynamic virtual wells produce with constant water production and the water yield is 200 m 3 /d. Variable permeability treatment is adopted near the working face. The permeability variation coefficient is jointly determined according to the stress field of the mine and the coal seam quality. The permeability within 10m of the treated working face is adjusted to the initial permeability 400 times (same treatment for vertical grids); the permeability is adjusted to 100 times of the initial permeability for the range of working face before and after 10-20m (the vertical grid is in the same way); permeability is adjusted to10 times of initial permeability for the range of working face before and after 20-50m range (the vertical grid is in the same way). As the work surface is moving, permeability of the mining active area (before and after 50m) is 1.5 times to the original permeability.
Based on the above given data and algorithms, using computer calculations, the calculation results considering mining conditions are drawn into the following CBM production curve over time (Figure 4), gas saturation distribution of reservoir at different times during the advancement of working face ( Figure 5), reservoir pressure distribution at different times during the advancement of working face ( Figure 6) and reservoir conductivity coefficient distribution at different times during the advancement (Figure 7).
It can be seen from Figures 4 to 7 that the simulation results of the improved coalbed methane seepage model are more consistent with the actual production data after processing by using the coupling method of dynamic virtual wells and variable permeability field. Affected by the advance of the mining face, the effective stress of the reservoir changes, so that the reservoir begins to desorb in advance within about 50m of the working face before the well is opened, and the gas saturation increases. At 50m after the advance of the working face, due to the re compaction of the reservoir and the reduction of pressure after coalbed methane desorption, the pressure of the gas reservoir remains relatively stable. Under the dual action of the advance of desorption of the working face and the production of production wells, the gas saturation presents two effects: higher gas saturation within the sweep range of the working face and lower production saturation near the gas well, the reservoir outside the advancing and spreading range of the working face is still dominated by normal pressure drop, saturation distribution and reservoir conductivity coefficient distribution.

Discussion
In the process of CBM numerical simulation technology, the inherent seepage mechanism, seepage law and adsorption and diffusion mechanism are the core of the underlying logic and mathematical equations for model establishment and subsequent software compilation. Therefore, the exploration and analysis of the laws and mechanisms of oil and gas seepage and diffusion adsorption in coalbed methane reservoirs is a long-term and important research topic.
Langmuir equation is selected as the adsorption capacity of coal bed gas in this model, which is widely used in the development of coalbed methane and shale gas reservoir, and has a good characterization effect on the desorption capacity of coalbed methane and shale gas. The diffusion of coalbed methane in the matrix of coal seam is based on the classical the first diffusion law of Fick. It mainly characterizes how the adsorbed CBM enters the larger fractured pore system from the smaller pore system with the decrease of pressure in the CBM matrix of the pore-fracture dual medium decreases. And it is the most commonly used characterization method in the development of coalbed methane and similar shale gas. There is no reliable mechanism model for the study of dynamic permeability under mining condition. The dynamic permeability field in this paper changes with the advancing distance of the mining face. At different time points, the permeability of each grid point in the mining area is different. An empirical formula is adopted in this model, which may have different application effects in different coalbed methane reservoirs.
In the CBM development industry, there is no clear understanding of the permeability change rule caused by the change rule of the structure and form of the coal matrix reservoir under the influence of mining, and the model building of the dynamic permeability change under the conditions of stress, temperature and coal seam desorption needs to be continued.
Under the mining condition of coalbed methane reservoir, with the change of pressure, because the coal seam itself belongs to organic matter, the coalbed methane desorption will lead to the change of particles in the coal seam, which has an impact on the numerical simulation. Different from inorganic minerals, coal is a complex mixture of organic and inorganic substances and will expand under low stress (Chen et al., 2021), leading to a great degree of influence of stress changes on numerical simulation. The development of coalbed methane is sometimes affected by formation water. Formation water will cause blockage in low-permeability coalbed methane reservoirs, which will affect coalbed methane reservoir property and make it difficult to develop the reservoirs. In addition, during the development of coalbed methane, due to the desorption of coalbed methane in the coal seam and the large elastic coefficient of the coal itself, the interlayer interference between multiple coal seams also has a great impact on the development, which may affect the development due to the great differences in the physical properties of reservoirs at different layers and depths in the development process. From the simulation results in this paper, the presence or absence of mining conditions has a significant impact on the effective development of CBM. From the simulation results in this paper, the presence or absence of mining conditions has a significant impact on the effective development of CBM. In fact, the main factor that can affect the development of CBM is the change of seepage performance in the coal seam as the reservoir temperature and pressure change, so the above-mentioned sensitive influencing factors should ultimately be attributed to the permeability performance of the coal seam. What kind of development approach can modify these influencing factors so as to effectively increase and efficiently utilize the coal seam permeability is the key to development research.

Conclusion
1. Combining the method of dynamic virtual wells with variable permeability field, CBM numerical simulation under mining condition is realized. 2. Mining conditions of CBM numerical simulation not only need to put into considerations mining working face impact on the permeability field, but also need to consider the mining of caused by the mining conditions area desorption to better characterize the mining conditions of CBM numerical simulation. 3. The mining face will cause the coal seam to desorb in advance, and the reservoir permeability increases first and then decreases before and after the work, resulting in a significant increase in the initial production of the gas well compared with the non mining well. After the advance of the mining face, although the reservoir permeability decreases, due to the existence of a large amount of desorbed free gas in the reservoir within the affected range, the gas well can still be produced with high production, but it will decline rapidly.