Coal and gas outburst, or outburst in short, is known as the uncontrolled ejection of coal and gas in underground coal mines, which occasionally happens along with mining activities. It is a destructive hazard that the probability of its occurrence has increased significantly by the increase of mining depth (Guo et al. 2017; Shadrin 2016; Yang et al. 2019). In Australia, 878 outburst accidents have been reported in history, and in China, 510 outburst events have been reported from 2001 to 2017, causing over 3576 tragic deaths (Black 2019a, b; Wang et al. 2019a, b). Figure 1 indicates the number of incidents reported in New South Wales, Australia, and China (Black 2019a, b; Sun et al. 2018). Outburst is a complex phenomenon because of the variety of its control parameters and finding an effective method to predict outburst probability using the most efficacious factors has been controversial.

Number of reported incidents in a Australia, New South Wales and b China (Black 2019a, b; Sun et al. 2018)

The effective factors contributing to outburst occurrence can be related to a large variety of coal and gas properties, including gas composition, gas content, gas pressure, temperature, adsorption/desorption rate, diffusion rate, absolute and relative permeability, porosity, coal seam thickness, coal type, coal strength, in situ stress, coal particle size, fracture toughness, and moisture content. Furthermore, the presence of geological structures and induced stress during mining can play a major role in triggering outbursts. Different methods have been suggested to determine a rational outburst threshold over the years to help the mining industry to manage outburst risks. For the Bulli seam in Australia, a linear relation between gas composition and gas content was proposed by Lama (1995) to define the outburst threshold, and for non-Bulli seams, Williams and Weissmann (1995) propounded the DRI900 index as the outburst threshold limit. Other countries use indices like ∆P_0–60, ∆P_Express, and K_T, which are based on the gas released from fresh coal sampled ahead of a working face (Black 2019a, b).

Apart from using the above methods to assess outburst risks, efforts were also made to explain the outburst mechanism. The theoretical model proposed by Litwiniszyn (1985) was one of the first attempts to explain outburst initiation. He assumed three-phase states (i.e., solid, liquid, gas) in a coal seam, in which the liquid phase is the result of capillary condensation of gas within small pores. By defining two kinds of shock waves within the micropores, i.e., compression and rarefaction, the author introduced rarefaction shock waves as the outburst trigger. Rarefaction shock waves have higher stress at the wave front compared to the behind, and Litwiniszyn (1985) also justified his outburst theory by assuming that stress changes induced by the propagation of rarefaction shock waves may cause the liquid phase suddenly transform into gas in pores, resulting in crushing the skeleton of the porous media and then outburst occurrence as a consequent.

Later, Paterson (1986) rejected this model because of its inability to explain outburst phenomena in non-coal geomaterials that have very weak or no adsorption properties. The lack of explanation about the source of shock waves and ignoring their attenuation due to the phase change were two other reasons that Paterson (1986) mentioned. Despite the disagreement, the author also assumed outburst as the structural failure of coal as a result of excessive stress and gas pressure and then modelled this phenomenon with three stages, i.e., gas flow, stress change and failure description. Paterson’s research can be counted as the first attempt to quantitatively model the initiation condition of outbursts (Paterson 1986). Later, in the work by Choi and Wold (2001), coal seams were considered gas reservoirs with dual porosity structure and the poroelasticity theory was used to simulate outburst initiation, which forms the basis of today’s studies on the numerical simulation of outbursts. Eventually, the energy aspect of outburst initiation and statistical analyses were added to the numerical models to yield a more comprehensive understanding of outburst initiation (Lu et al. 2019a, b; Soleimani et al. 2023; Wold et al. 2008; Xue et al. 2015; Zhao et al. 2022, 2020). With the development of numerical modelling methods and software packages, the capacity of using numerical simulation to investigate outbursts by considering a wide range of control parameters has been significantly improved over the years.

Given the recent development in new numerical modelling methods, it is crucial to have a comprehensive review of different methodologies, their strengths and weaknesses, and clarification of future research gaps on implementing them in outburst modelling. In this paper, the effective control parameters in outburst occurrence are first summarised into four categories (gas properties, coal properties, stress regimes and geological structures), and then various numerical methods to assess outburst occurrence are reviewed. Since the main focus of this paper is on numerical methods in outburst modelling, all governing equations to describe the involved physics in an outburst process and different possible scenarios (sorption effect, two-phase flow, stress effect on the permeability, etc.) are analysed. Previous studies on the numerical modelling of outbursts are divided into two groups (the iterative coupled group and the fully coupled group), and then their detailed modelling procedures, innovations and deficiencies are elucidated. Since most of the previous works have focused on the coupling of geomechanical and fluid flow equations and failure/damage evolution in porous media, there is a lack of analysis on quantitively assessing the moment of outburst initiation. In the last section of the paper, ideas for future works are presented from three perspective of analysis the phenomenon from particles viewpoint, conservation of energy and statistical methods. Figure 2 demonstrates the workflow of this paper.

Workflow of the paper

Generally, the key control parameters in an outburst event can be categorized into four groups: gas properties, coal properties, stress regimes, and geological structures (Fig. 3). What makes outburst studies complicated is that these parameters should not be considered separately since they are not independent of each other. For instance, gas content is related to porosity and sorption capacity while changes in the stress and gas pressure can affect both porosity and the amount of adsorbed gas. Furthermore, the presence of geological structures may cause different rock mechanical properties and thus different stress distributions.

Effective control parameters in coal and gas outburst

2.1 Gas properties

2.1.1 Gas content

One of the factors affecting the energy release potential of coal is gas content. Gas content is the summation of free and adsorbed gas in coal and is recognized as the key parameter used for determining outburst propensity (Lama 1995). In Australia, gas content is used alone to define a Threshold Limit Value (TLV) in coal mines to assess outburst risks. Although the proposed TLV defines a lower limit for CO₂-rich coal compared to CH₄-rich coal, numerical modelling and field data do not fully justify this (Choi and Wold 2002). Black (2018) compared gas content measurement data from outburst-prone locations in the Bulli seam and concluded that a 9 m³/t TLV can be used in coal mines regardless of gas composition (see Sect. 2.1.2 below).

Gas content is normally measured in three defined steps (Fig. 4). In the first step, the free gas volume Q₁ is estimated from the measurement data of the initial gas desorption rate at the field. Q₁ represents the amount of escaped gas or lost gas during sample retrieving and before sealing is available, which cannot be measured directly. In the second step, the desorbed gas volume (Q₂) is measured from the sealed intact samples. Then the coal samples are pulverized using a crusher bowl in the third step, and the fast-desorbed gas volume (Q₃) representing residual gas would be measured under atmospheric pressure. The summation of Q₁, Q₂ and Q₃ are called the measured gas content (Q_m) (Diamond and Schatzel 1998; Saghafi 2017).

Apparatus used to measure gas content (Saghafi 2017)

2.1.2 Gas pressure

Gas pressure is a key indicator of the potential of energy release and consequently the outburst proneness of a coal seam (Black 2018). In China, there is a TLV of 0.74 MPa for outburst-prone mines, while in Australia, Lama (1983) declared a 0.6 MPa threshold of gas pressure at five-meter deep from the coal face. However, this TLV has not been widely adopted in Australian mines. The observed sharp pressure changes versus depth (in Fig. 5) is believed to be one of the main reasons for outburst initiation.

Gas pressure and gas content gradients behind the coal face – Gemini seam, Leichhardt colliery (Black 2019a, b; Hanes 1995)

2.1.3 Gas composition

Investigations from different countries with reported coal and gas outbursts cases, like Australia, Canada, China, France, Turkey, etc., show that the dominant seam gas in an outburst event is methane (CH₄), carbon dioxide (CO₂) or a combination of both. Previous research shows that the sorption capacity of CO₂ is higher than CH₄ in coal, and thus CO₂-rich coals are in the risk of more severe outburst at the same reservoir pressure (Black 2018; Si 2015). On the other hand, at the similar level of gas content, CH₄-rich samples would have higher gas pressure compared to CO₂-rich ones, which means higher risks to trigger outburst initiation. In addition, the conclusions from Holda (1986) and Lama (1995) about the softening behaviour of coal in conjunction with the presence of CO₂ was questioned by Choi (2014), and later by Black (2018) after studying field data from different seams in Australia. They believed that there is no clear connection between gas composition and outburst occurrence (Fig. 6). Based on these investigations, although CO₂ is more preferably adsorbed by coal compared with CH₄, its role in affecting the initiation of outburst is unclear.

Recorded gas content and composition at outburst-prone sites in Australian underground coal mines (Black 2018)

2.1.4 Sorption capacity

The main difference between sandstone gas reservoirs and coal is the gas sorption behaviour. Gas sorption can be caused by adsorption on the surface of pores and cracks and even penetration into coal structures (absorption). This factor is related to other parameters, such as coal rank, gas type, pressure, temperature, moisture, and ash content, which is normally described by Langmuir isotherms.

Although sorption is the main gas storage mechanism in coal seams, its contribution to an outburst incident is debatable as its effect is associated with diffusion, which is a slower phenomenon compared to outburst initiation, which is an immediate process. But in coals with smaller particles, gas may migrate through a shorter path to the external environment, and thus the effect of desorption would be observed faster compared to larger coal particles (Li et al. 2019). The inverse relationship between the coal particle size and desorption rate could be the reason for the higher outburst propensity in approximate to geological structures with pulverised coal. Desorption rate is also the basic parameter for defining many outburst indices, for instance, the DRI, ΔP Index, IDR300, Ettinger’s Sorption/Desorption Index, Hargraves Emission Value, KT index, V30 index, etc. (Bodziony and Lama 1996).

2.1.5 Diffusion coefficient

Gas flow obeys different physics in the coal matrix and through cleats, as shown in Fig. 7. In the coal matrix, gas flow is driven by diffusion through micropores, which is controlled by the concentration difference due to gas desorption from the pore surface of coal.

Schematic of gas flow in coal (Ahsan 2006)

Gas sorption is a complicated process that causes difficulty in determining the diffusion rate of desorbed coal seam gas. Although various experiments have been suggested to determine the diffusion coefficient of coal (Cavé et al. 2009; Cormenzana et al. 2003; Xu et al. 2015), inconsistent measurement results with a wide range of variation made it difficult to study its effect on outburst initiation. Recently, micro-CT imaging was suggested as a reliable method to assess diffusion coefficient in coal. This method can be conducted under high pore pressure and confining stress that would yield a more accurate understanding of gas transport under pre-excavation stress distribution (Lu et al. 2020).

2.2 Coal properties

2.2.1 Coal strength

Coal and gas outburst can be defined as the structural failure of coal, which is followed by the ejection of failed coal particles by depressurised coal seam gas (Choi 2014; Paterson 1986). To prevent this incident, coal needs to withstand the drag force and pressure gradient force exerted by gas. The presence of tectonic structures, reduction of confining stress, and increase of abutment stress during mining can weaken the coal structure, which may explain outburst occurrence near geological structures or ahead of mining faces. The development of fractures/cleats and lack of ductility can also affect coal strength, and thus the post-peak behaviour and strain-softening of coal should be considered in outburst analysis. Furthermore, the weakening of pulverised coal and easier transportation of smaller coal particles by gas make the coal damage behaviour and fragmentation as an important factor to be considered in outburst studies.

Despite conventional laboratory methods to determine coal strength (e.g., uniaxial and triaxial testing), the extreme low strength of outburst-prone coal prevents it to be cored and tested in laboratory conditions. In China, as an alternative to the uniaxial compressive strength, the Protodyakonov coefficient (firmness index) is suggested for studying the strength of outburst-prone coals. This index measures coal vulnerability to be crushed by dropping a 2.4 kg hammer from the height of 600 mm. Particles smaller than 0.5 mm would be the representative value to be used for calculating the Protodyakonov coefficient (Black 2018).

2.2.2 Coal rank

Coal rank is the expression of coal maturation defined by its physical and chemical properties (Fig. 8). Although previous studies revealed that there is a certain correlation between the coal rank and outburst propensity (due to its effect on coal characteristics such as gas sorption), no specific thresholds regarding the required coal rank to trigger outbursts have been defined (Das et al. 1991; Hedberg 1980; Patching 1970). Beamish and Crosdale (1998) concluded that the higher vitrinite content leads to the higher outburst propensity, which also results in smaller particle sizes. Also, Australian mining records show that outbursts usually occur in bituminous coal mines (vitrinite reflection > 1.2%). But there are also reported outbursts in the Velenje coal mine, Slovenia, which mines lignite coal (Bodziony and Lama 1996; Si et al. 2015a, b).

Classification of coal rank proposed by American Society for Testing and Materials (ASTM) (Schweinfurth 2003)

2.2.3 Porosity and permeability

The pore size in coal can be categorised as Table 1, which also demonstrates the inverse relationship between the coal rank and pore size. Figure 9 illustrates the porosity response to pore pressure and stress change, as well as how the fluid and solid interacts with each other. A cubic relationship between porosity and permeability for coal cleats was suggested by Mckee and Hanson (1975).

The interaction of solid and fluid in porous media

Following gas diffusion, permeability is the key parameter controlling the second stage of gas flow in coal seams, which is driven by the pressure gradient through cleats. Lama and Bodziony (1998) concluded that coal with permeability higher than 5 mD is not outburst prone, whereas low permeability or tight coal seams, may lead to high pressure difference between the mining face and deep coal inside the wall, indicating much higher outburst risks. The gas pocket theory is based on this fact, which assumes outburst normally happens in geological structures with low permeability and high gas content (Lu et al. 2011). In-situ well-testing and core flooding are the two methods normally used for permeability measurement. Permeability at different confining stresses can be measured by the apparatus shown in Fig. 10.

Experiment setup for stress-dependent permeability measurement at UNSW

2.3 Stress

2.3.1 In-situ stress

In-situ stress is one of the main parameters controlling the dynamic failure of coal, particularly if a large amount of elastic strain energy is stored in coal. This parameter can be dominative that the ejection of coal particles might happen even without the involvement of any gas, which is known as coalbursts or rockbursts in general. Rock failure in rockbursts normally results in the ejection of the elastic strain energy with high stress on brittle rock materials. Not all of this strain energy would be converted to kinetic energy for the ejection of coal particles. Plastic deformation, fracture opening, and seismic wave emissions also consume a certain portion of the strain energy. The presence of gas and its contribution to either failure or ejection of coal make outburst a much more complicated process.

In-situ stress can be analysed in two categories: vertical and horizontal stresses. Vertical stress is determined by the mining depth. Deeper coal mines are considered more prone to outburst because of more vertical stress being exerted on coal, with higher gas pressure, lower permeability, more maturity of coal, and higher virgin rock temperature. Mining depth plays a significant role in steeply dipping seams and flat seams which has the vertical stress as the major principal stress is the vertical stress (Bodziony and Lama 1996). On the other hand, the strength of coal will increase by increasing the confining pressure. In addition, the relationship between coal's mechanical properties and confining stress can also be strongly influenced by the presence of gas and its pressure, making the study of this phenomenon more complex (Ding et al. 2021).

It is acknowledged that areas with high tectonic stress, where coal may undergo mechanical degradation, are more susceptible to outbursts by mainly affecting horizontal stress. In the case of soft coal, when the high-stress-induced mechanical degradation process causes block coal to fracture into small particles, an outburst can be easily triggered by gas energy. On the contrary, in the case of hard coal, although high stress cannot break solid coal into particles, a significant amount of elastic strain energy can be stored, which can also lead to dynamic hazards (Li et al. 2007).

2.3.2 Mining-induced stress

Coal extraction changes the stress distribution around mine openings. The reduction of confining stress during mining reduces coal strength and expands the mining damage zone, resulting in fractures at the face and corners of mine openings. The superimpose of in-situ and mining-induced stresses can largely increase the outburst propensity. This phenomenon has been commonly observed in rockburst or coalburst mines, which can be categorised into three types based on the main causes of high stress concentration. The first category is “high in-situ stress” rockbursts occurring near geological structures. The second category is “high mining-induced stress” rockbursts influenced by factors such as coal pillars and isolated working faces. The last category is “high in situ stress plus high mining-induced stress” rockbursts involve both geological structures and coal pillars contributing to critical stress (Fig. 11).

Schematic of stress contribution in rockbursts/coalbursts (Guo et al. 2017)

Cai et al. (2021) studied on coal bursts triggered by fault structures under both static and dynamic stresses caused by mining activities and discussed the interaction between mining activities and geological faults, which plays a key role in fault reactivation and subsequent coal bursts. The study proposed two types of fault reactivations: mining-induced quasi-static stress (FRMSS)-dominated and seismic-based dynamic stress (FRSDS)-dominated. These reactivations were validated through experiments, numerical modelling, and in situ microseismic monitoring. In this paper, four conceptual models of fault reactivation induced by underground coal mining activities. First when mining activities are far from the fault, and there is no direct interaction between fault stress and abutment stress. Fault reactivation is unlikely to occur, but some local deformation or transient slip might be triggered by far-field mining-induced seismicity. The second model was mining activities approach the fault, causing a superposition of fault stress and abutment stress, resulting in a high static stress concentration in the fault pillar. This change in stress state triggers fault reactivation dominated by mining-induced quasi-static stress (FRMSS). The third model covered mining activities shift away from the fault, leading to unloading of the abutment stress in the horizontal direction and gradual loading in the vertical direction. This can trigger an FRMSS-dominated fault reactivation if the fault pillar is not wide enough. The last model was about Mining activities shift parallel to the fault. If the fault pillar is wide enough, there may be no significant interaction between fault stress and abutment stress. Otherwise, fault reactivation can be triggered by an extra driving force, such as mining-induced seismicity, resulting in an FRSDS-dominated fault reactivation.

The results indicate that fault-induced coal bursts are triggered by the combination of high static stress in the fault pillar and dynamic stress from fault reactivation. Experimental validation using cylindrical sandstone samples confirmed the generation of stick–slip phenomenon in the fault plane. Numerical modeling and microseismic monitoring demonstrated high stress concentration in the fault pillar, and seismic events were observed along the fault plane and in the surrounding medium.

2.4 Geological structures

Shepherd et al. (1981) introduced geological structures (Fig. 12) as the main parameter in outburst occurrence. These structures contributed to about the 98% of outburst cases reported in Australia (Lama 1995). The most common geological structures in coal mines are folds, dislocation systems (like faults, joints and discontinuities), and magmatic intrusions. These structures can cause the concentration of tectonic stress and increase the accumulated elastic strain energy, and gas pressure gradient may also increase because of gas accumulation in low permeability barries (as explained by the gas pocket theory shown in Fig. 13).

Illustration of geological structures in outburst-prone coal seams (Shepherd et al. 1981)

Schematic of an outburst induced by the gas pocket structure (Li 2001)

The traditional methods to assess outburst risks are mostly experiment-based, which measure the total gas content or gas released from coal samples during a specific time. In Australia, Lama (1995) proposed the TLVs for coal mines based on previously occurred outbursts in the Bulli seam. He stated that gas content is the key factor in outburst occurrence and proposed a linear TLV based on gas content and gas composition (Fig. 14). Any coal seam conditions above the TLV line will be considered outburst prone.

Proposed outburst threshold limits based on outburst events in the Bulli seam (Lama 1995)

The volume of gas released in the first 30 s of the third stage of gas content measurement (Q₃ measurement) from a 200 g sample is called the Desorption Rate Index (DRI). DRI has a linear relationship with the gas content Q_m, and the slope of which depends on gas composition. Williams and Weissmann (1995) reported that DRI900 is aligned with the suggested thresholds for CH₄ and CO₂ rich coals in the Bulli seam, and with further study, they reported that DRI900 could be also applied as the outburst TLV for non-Bulli seams. A recent study on the collected data during 2017 from Bulli seam revealed that the criteria defined based on this seam has slightly shifted compared to the results reported in 1995 (Black 2019a, b).

Ettinger (1952) proposed Ettinger’s Sorption/Desorption Index, which can be referred as the foundation for all emission rate methods. Using the original method proposed in 1952, samples with higher adsorption rate would be considered outburst prone. The procedure was to crush the sample (particle size: 0.25–0.50 mm), dry it at 60 °C, and then vacuum it to remove all remaining gas. The amount of gas release in the percentage of sorption capacity at the absolute pressure of 1 atm and 30 °C was plotted against time. This method was later modified in 1953, and the rate of gas emission while the pressure increases for a sample in an enclosed chamber was taken as the new index, i.e., gas emission index (Bodziony and Lama 1996).

To determine outburst risk, ΔP_0–60 index is a method to measure the pressure build-up during the first 60 s of desorption while the sample is sealed in a desorption chamber. In Germany, 15 mmHg (2 kPa) is considered as the outburst TLV. The presence of ash can affect this index, so the amount of ash content should be corrected by Eq. (1) (Black 2018; Bodziony and Lama 1996).

K_T index refers to the slope of the log–log plot of the gas desorption rate versus desorption time to represent changes in desorption rate. Paul (1977) proposed the ΔP_express index to address the issue of lacking equipment for determining K_T in underground mines. In that method, 70 g of coal samples with 0.25–0.50 mm particle size for bituminous coal (or 2–3 mm particle size for anthracite coal) are enclosed in a chamber to be evacuated for 2 min. Then methane is injected into the chamber and let the pressure raise to 2 bar. By instantly connecting a manometer to the chamber, the gas pressure after 1 min is considered the ΔP_express index (Bodziony and Lama 1996; Esen et al. 2018). Threshold valued for Ettinger, ∆P_0–60 and K_T indices are shown in Table 2.

Other methods that are based on gas emissions from coal samples in the laboratory or during drilling in the field include the V₃₀ index, V_index, Hargrave's emission rates, gas emission rate, and gas flow index (Bodziony and Lama 1996). They are all different varieties of the abovementioned indices and will not be elaborated here.

According to a recent report from CSIRO, in 2019 China National Coal Mine Safety Administrator represented a procedure through which 3 steps as a combination of different methods were defined. In the first step, four indices are used to determine the risk of outburst. When all four indices fall within the specified ranges in Table 3, the coal mine is considered outburst prone. But, for gas pressure and coal firmness if the threshold values meet the values in Table 4 the risk of outburst occurrence is considered.

After determination of outburst prone mines using these four indices, a map illustrating the outburst risk within these mines is generated using gas content or gas pressure. During this stage, the threshold values need to be determined through experiments. However, in the absence of conducted tests, the threshold value for relative pressure is set at 0.74 MPa, while for gas content, it is considered to be 8 m³/t for unstructured zones and 6 m³/t for structured zones.

Finally, the outburst risk ahead of the working face is evaluated using localized prediction. This involves measuring the volume of borehole cuttings, determining the rate of desorbed gas from borehole cuttings, and the initial rate of gas flow from the borehole (Qu et al. 2020).

With the development of high-capacity computational technologies, numerical modelling has become an essential approach for solving engineering problems. Since underground excavations and fluid flow in porous media were not exceptions, numerical modelling is also a powerful tool to analyse coal and gas outbursts.

Choi and Wold (2001) suggested studying the outburst process in 3 steps, namely pre-initiation, initiation, and post-initiation. As the first step, pre-initiation is the study of coal behaviour while quasi-static deformation happens. The initiation step is when the dynamic failure happens and can be counted as the last part of the pre-initiation modelling. The post-initiation step is the study of outburst intensity as the dynamic failure and fragmentation of coal.

4.1 Fluid flow and solid mechanics coupling in porous media

The first attempt of theoretical research on the dual-porosity media with single-phase fluid flow can be attributed to Barenblatt et al. (1960). In that study, fractures are high permeable channels that are responsible for most fluid flow. On the other hand, pores provide storage space to supply fluid for fractures. They assigned two pressures and two velocities to each pore structure (one for pores and one for fractures). Also, they considered fluid transfer between pores and fractures and the mass conservation law for both pore structures. After Barenblatt et al. (1960), other studies were conducted to understand the behaviour of dual-porosity media. However, these studies lacked investigations into the coupling behaviour of solid mechanics and fluid flow, and the interaction of fluid and solid was not considered (Kazemi 1969; Odeh 1965; Warren and Root 1963). Since effective stress is related to the vertical stress and fluid pressure in the pore space, any fluid movement can cause changes in effective stress, and consequently pore deformation and pore volume change. Thus, the fluid pore pressure changes are also responsible for the fluid movement in porous media (Sibernagel 2007). To analyse such material, models were proposed with the mixture of the coupled poroelasticity theory and dual porosity concept, as suggested by Biot (1941) and Barenblatt et al. (1960), respectively. Therefore, to study the geomechanical behaviour of coal seams, the force equilibrium equation was used to analyse the stress exerted on a unit bulk volume of solid, which would have the internal force equilibrium instantly after any stress changes. (Sibernagel 2007). Also, to study the fluid flow behaviour, coal is assumed as a porous medium containing water and gas, which is normally described as a dual-porosity system. Choi and Wold (2001) mentioned that a coal seam can be treated like a gas reservoir. But the key difference between coal seams and conventional gas reservoirs is the gas sorption characteristics of coal, which empowers it with a different gas storage mechanism to act as both the reservoir and source rock. The fluid transport mechanism in the coal matrix is diffusion, controlled by both sorption and permeability due to the critical mass transfer rate between fractures and the matrix. The gas flow in cleats obeys Darcy’s law, and the gas sorption behaviour can be indicated by Langmuir isotherms (Harpalani and Schraufnagel 1990). Coal permeability changes are related to effective stress and adsorbed gas, which should be also considered in calculations (Aminian et al. 2004; Clarkson and Bustin 1999a, b; Wu et al. 2010).

4.2 Governing equations

Wilson and Aifantis (1982) proposed the theory of consolidation with dual-porosity and single-phase flow constitutive equations as Eqs. (2)–(4).

where \({\sigma }_{ij}\) indicates the component of the stress tensor, \({\varepsilon }_{ij}\) indicates the component of the strain tensor, \({\varepsilon }_{\text{kk}}\) is the volumetric strain, G and \(\lambda\) are Lamé Parameters, \({\delta }_{ij}\) is Kronecker Delta, \(\alpha\) is Biot Coefficient, \(\beta\) is a coefficient based on compressibility and porosity, µ is the dynamic viscosity, ω is the transfer coefficient between fracture and matrix (Khalili et al. 1999; Wu et al. 2010), P represents pressure and k represents the permeability. Subscripts m and f indicate the coal matrix and fracture, respectively. In some studies, the term containing the total strain was not mentioned and porosity change has been used in the equation as an alternative (Wu et al. 2010). These equations can be used for every dual-porosity medium, but for coal, the stress equation has an additional volumetric strain induced by sorption (Eq. (5)), which is defined as Eq. (6) (Wu et al. 2010). In these equations, K represents the bulk modulus, and ε_L and P_L are the Langmuir volumetric strain and Langmuir pressure constant, respectively.

In Eqs. (2)–(4), the components of the stress tensor and strain tensor are indicated by \({\sigma }_{ij}\) and \({\varepsilon }_{ij}\), respectively. \({\varepsilon }_{kk}\) represents the volumetric strain, while G and \(\lambda\) are known as the Lamé Parameters. The Kronecker Delta is represented by \({\delta }_{ij}\), while the Biot Coefficient is represented by \(\alpha\), and the coefficient based on compressibility and porosity is represented by \(\beta\). Dynamic viscosity is represented by µ and the transfer coefficient between the fracture and matrix is represented by ω (Khalili et al. 1999; Wu et al. 2010). Pressure is represented by P and permeability is represented by k, with subscripts m and f denoting the coal matrix and fracture, respectively. In some studies, the term containing the total strain was not mentioned and porosity change has been used in the equation as an alternative. For porous media like coal with sorption characteristics, an additional volumetric strain induced by sorption (\({\varepsilon }_{\text{s}}\)), as defined by Eqs. (5) and (6), should be included in the stress equation (Wu et al. 2010). In these equations, the bulk modulus is represented by K, the Langmuir volumetric strain by \({\varepsilon }_{\text{L}}\), and the Langmuir pressure constant by \({P}_{\text{L}}\).

In a porous medium like coal with sorptive characteristics, an additional term named sorption-induced stress (\({\sigma }_{\text{s}}\)) should be included in the stress equation (Eq. (5)). The most common way to define this term in the outburst study is based on the phenomenological nature and does not conform to the laws of thermodynamics. Two other methods are equilibrium thermodynamic and non-equilibrium thermodynamic. The former approach can be used for gas pressures in a state of equilibrium, but the model does not take into account the change in pressure over time. In the coal system after the transfer of mass and energy with the surrounding environment, it cannot be brought back to the primary state, so it would be a matter of non-equilibrium thermodynamics (Espinoza et al. 2013, 2016; Lv et al. 2022). Lv et al. (2022) compared the three methods and expressed that the non-equilibrium thermodynamic method is the best one regarding the realistic estimation and convenience of the calculation. Furthermore, Soleimani et al. (2023) applied the phenomenological (Eq. (6)) and non-equilibrium thermodynamic methods (Eq. (7)) into assessing the outburst occurrence. Although the results of both methods were compatible, the higher amount of sorption-induced stress calculated by the phenomenological method led to observing the outburst in the earlier stages of the excavation. Also, in Eq. (7), \(z\) and \(R\) are the compressibility factor and gas constant, \(T\) represents the temperature, \(C\) denotes the coupling coefficient, and \({n}_{\text{L}}\) is the Langmuir coefficient.

Equations (2)–(4) represent the fluid flow in a dual-porosity and dual-permeability medium. In coal, the permeability difference between the matrix and fractures creates a dual-porosity and single permeability medium. Furthermore, Fick's diffusion equations can be applied to calculate the matrix pressure. In this case, Eq. (8) expresses the mass transfer per unit volume between the matrix and fracture (q_m), indicating the mass change of gas presenting in the coal matrix over time (Lim and Aziz 1995). M_g,m is the amount of gas in the matrix, which is the summation of free gas and adsorbed gas in the pores (Eq. (9)) (Langmuir 1918). In Eqs. (8) and (9), D_e is the effective diffusion coefficient, c is the concentration of free gas, φ is the porosity, M_W and V_M are gas molecular weight and volume, R is the gas constant, ρ_c is the coal density and V_L represents the Langmuir volume constant.

The effect of moisture can be added to Eq. (5), as shown in Eq. (10) (Chen et al. 2012):

In the above equation, λ_m is the adsorption capacity decay coefficient with moisture content. M is the moisture content in the percentage of unit weight, \({\varepsilon }_{\text{L}m}\) represents the maximum moisture sorption-induced swelling strain and m_L is the Langmuir strain constant for moisture sorption strain at the half of \({\varepsilon }_{\text{L}m}\). In addition, m_c represents the equilibrium moisture content, which is also in the percentage of unit weight.

In the presence of moisture, two Darcy equations (Eqs. (11) and (12)) are governed for either water or gas flow in cleats. The coupling between these two equations is based on the capillary pressure, and relative permeabilities that are associated with the saturation of each phase.

In Eqs. (10) and (11), ρ, S, q, h and k_r represents density, saturation, mass production rate, the height of the coal seam and relative permeability, respectively, with the subscript g for gas and w for water.

4.3 Coupling fluid flow with stress using permeability correlations

There are multiple factors that can control the permeability of coal, such as the fractures and pore structure, stress state and mechanical properties, swelling/shrinkage due to gas adsorption/desorption, damage, failure, etc. The initial permeability of coal is largely dependent on the fractured nature of coal (Fig. 15), i.e. the connectivity of cleats. Changing the gas pressure regime and stress state in coal (by excavation or gas drainage) will alter the permeability of coal. Thus, to fully couple the multiple physics involved in outburst modelling, permeability should be expressed as a function of effective stress.

Anisotropic cleat system (Liu et al. 2018; Pan and Connell 2012)

As shown in Fig. 15, the porosity of fractures in an anisotropic medium can be calculated by Eq. (13). Also, considering that gas flow in each direction is analogue to the flow between parallel plates, Eq. (14) for permeability calculation can be used (Bai and Elsworth 2000; Pan and Connell 2012).

To simplify the calculation, the medium can be assumed as isotropic and the matrix deformation can be neglected in comparison to the aperture change of cleats (i.e., the spacing between fractures \({a}_{i}\) is a constant while the fracture aperture \({s}_{i}\) is a variable). Therefore, the popular cubic relation between porosity and permeability can be concluded (Eq. (15)). \({k}_{f0}\) and \({\varphi }_{f0}\) represent the permeability and porosity of fractures at the initial stage, respectively.

In outburst modelling, the effect of stress on porosity and permeability is of great importance. Since the basic definition of porosity is the pore volume divided by the bulk volume, its changes can be expressed as a function of bulk and pore volumetric strains. Another basic definition is compressibility, which represents the volumetric change as a response to stress change (Durucan and Edwards 1986; Harpalani and McPherson 1986; Shi and Durucan 2004).

When it comes to coal, the gas adsorption/desorption effect should also be considered in permeability calculation, and thus, both mechanical and sorption-induced strain on permeability change should be included simultaneously in outburst modelling. To this end, a Langmuir-based equation is usually used to calculate the sorption-induced strain for coal, and then this effect will be reflected in permeability change (Harpalani and Schraufnagel 1990; Levine 1996). Another suggested method is using surface potential energy calculated by thermodynamic approaches (Pan and Connell 2007; Vandamme et al. 2010). The permeability models are usually based on the elastic state of coal, while during the complex outburst process, due to excavation-induced stress and gas pressure changes, inevitability, there will be plastic deformation. In addition, the changes in existing fractures and the evolution of new cracks would lead to damage and failure in the coal structure. Thus, damage-based permeability models should be used in outburst simulation, which will be introduced below.

Si et al. (2015a, b), introduced a coefficient in the stress-based permeability equation to independently consider fracture permeability. Later, Chen et al. (2016) suggested applying a logistic growth function as the coefficient in the equation that relates permeability to the exponential of compressibility (Eq. (16)). They introduced a variable L as an “S” shape function of deviatoric stress (Eq. (17)) and then fitted that function using experiment data.

In Eqs. (15) and (16), \({c}_{\text{f}}\) is the cleat compressibility, \(\tilde{\sigma }\) and \({\sigma }_{\text{d}}\) are mean effective and deviatoric stresses, respectively, and ζ, η, γ are fitting parameters.

Zheng et al. (2018) applied the scalar damage variable (\(D\)) into the cubic relationship of permeability and porosity (Eq. (18)). \({\alpha }_{\text{k}}\) is the coefficient of damage-permeability effect and can be assumed to be 5 (Zhu and Wei 2011).

To describe the stress–strain relationship in a deep coal seam under excavation, Zheng et al. (2018) defined damage initiation criteria for tensile and compressive failure using the maximum tensile stress and Mohr–Coulomb criteria (Fig. 16), respectively, as defined in Eqs. (19) and (20), where, \({f}_{\text{t}0}\) and \({f}_{\text{c}0}\) represent uniaxial tensile and compressive strengths, respectively, and \(\theta\) is the internal friction angle.

The stress–strain relationship in coal (Xin et al. 2021; Zheng et al. 2018)

Then the damage variable was calculated using Eq. (21), where \({\varepsilon }_{\text{t}0}\) and \({\varepsilon }_{\text{c}0}\) are the maximum principal strains in tension and compression after damage occurrence, respectively.

An et al. (2013) and Lu et al. (2019a, b) modified the Palmer and Mansoori model by using the equivalent plastic strain (γ^P) for coal in the post-peak stage, which contains both the strain-softening and residual stages (Eq. (22)). In their model, there is a permeability reduction over the elastic stage, due to the decrease of cleats aperture caused by stress increase. Then in the post-peak stage, due to the evolution of new cleats, permeability will increase at first, but it would remain mostly constant in the residual stage.

In Eq. (23), γ^P* is the equivalent plastic strain at the beginning of the residual stage, \(\xi\) is the sudden increase coefficient of permeability, \(\Delta {\varepsilon }_{\text{m}}^{\text{s}}\) is the change in the adsorption strain of matrix (\({\varepsilon }_{\text{m}}^{\text{s}}=\frac{{\varepsilon }^{\text{s}}}{1-{f}_{\text{m}}}\)), and f_m represents the internal swelling ratio of matrix adsorption deformation.

Generally, for the numerical modelling of multiphysics interactions in a dual-porosity medium, four approaches are used (Settari and Walters 2001).

1. (1)

   Decoupled models In this method, pore compressibility as a scalar value is considered in the fluid flow equation, which forms the basis of many reservoir simulators. In this approach, fluid flow and solid deformation are not actually coupled, and the response of porous media to stress change is not captured (Minkoff et al. 1999a, b).

2. (2)

   Explicitly coupled models In this method, fluid flow is calculated in each timestep, and geomechanical computation is performed on a different timestep, which is based on the changes in pore volume obtained in the previous timestep. This approach can be counted as a specific case of iterative coupling that contains one iteration per timestep (Dean et al. 2006; Minkoff et al. 1999a, b; Minkoff et al. 2003; Settari and Walters 2001)

3. (3)

   Iteratively coupled models In this method, the fluid flow and geomechanics calculations are performed sequentially. In each timestep, the calculated pore-volume at the end of nonlinear iterations couples with the fluid flow and deformation (Dean et al. 2006; Fung et al. 1994; Settari and Mourits 1994).

4. (4)

   Fully coupled models In this method, equations of fluid flow and geomechanics are solved together at once. This technique (which is also called implicit coupling) is the most stable, but slowest in comparison to the explicit and iterative methods (Aifantis 1980; Dean et al. 2006; Lewis and Sukirman 1993; Li and Zienkiewicz 1992; Tortike and Ali 1987)

Dean et al. (2006) compared the results of explicit, iterative and fully coupled models through four different scenarios. This included two simple reservoir depletion scenarios for single-phase flow with different displacement boundary conditions, a three-phase flow scenario in which a production well and an injection well are in two opposite corners, and a scenario originally proposed by about a soft reservoir surrounded by stiff non-pay regions. In the first three scenarios, the model results from all three methods are consistent. However, in the last one, the reservoir simulator was unable to consider the geomechanical effect on the boundaries, and the same results were only obtained by using tight iterations and small timesteps in the iteratively and explicitly coupled methods. In terms of outburst simulation, to capture all interactive physics, the most commonly used methods are the iteratively and fully coupled models, which will be reviewed in detail in this paper.

5.1 Fully coupled outburst models

In fully coupled outburst models, geomechanical and gas flow equations are solved at the same time using one simulator. One of the early fully coupled models was developed by Paterson (1986), which simulated gas flow, stress and failure as rock physical responses associated with outburst initiation. In the gas simulation part, a single-phase, isotherm gas flow that obeys Darcy’s law was considered, and permeability was assumed to be an average value, independent of the effective stress. For stress description, the equilibrium equations for a perfectly elastic structure were used. For the failure part, tensile failure was considered the main reason for outburst occurrence. By defining the maximum pressure gradient that coal can resist at each point, rock failure and stress equations were fully coupled. The author also suggested that outbursts can be prevented by reducing the gas pressure gradient ahead of the mining face below a certain threshold value. Despite the deficiencies in the oversimplified assumptions, such as the lack of considering coal sorption behaviour, moisture effect, and stress-dependent permeability, this study opens an active research area of using numerical models to analyse the complex behaviour of outburst initiation.

Xu et al. (2006) considered the three outburst stages suggested by Paterson (1986) and used the RFPA2D-GasFlow simulator to model the induced outburst during crosscutting or intersecting a fresh coal seam during the drilling, which is known as a high risk pattern for the outburst occurrence. In their model, Darcy’s law was used for gas flow modelling, and Terzaghi’s law was used to calculate effective stress. For rock failure description, they considered two aspects, namely the damage induced stiffness degradation and damage induced permeability increase. Based on the numerical modelling results, they divided the outburst occurrence process into four stages: (1) Stress concentration ahead of the coal face, (2) Coal/rock fracturing and splitting induced by stress, (3) Crack propagation driven by gas pressure, and (4) Coal ejection carried by high-pressure gas (Fig. 17). They concluded that the key parameters causing outburst initiation are in-situ stress, gas pressure, and the properties of coal and rock. This model successfully simulated the spatial variability of outburst initiation, and more importantly, the dynamic failure process of coal triggered by stress with the superimposing of gas pressure has been well-captured. The main drawback was that the gas desorption effect in coal was not considered and lack of calculations regarding the mass transfer between the coal matrix and cleats.

Results of modelling coal and gas outburst due to mining in a coal seam (Xu et al. 2006)

Fan et al. (2017) used COMSOL Multiphysics to develop a fully coupled stress-seepage-damage outburst model (Fig. 18). They simulated the dynamic evolution of outbursts and analysed the interaction of coal–gas properties, geological structures and mining disturbance. They used Darcy’s law for gas flow simulation and updated the adsorbed gas volume based on Langmuir isotherms in the calculation. However, the matrix pressure was not updated in the simulation and the same pressure value was used for the fracture pressure. But the fluid flow mechanism is not the same in the matrix and fractures, so their pressure values are not the same while the mass transfer happens between them. They also applied the damage elastic modulus in the solid mechanics analysis and used the maximum tensile stress and Mohr–Coulomb criteria for tensile and shear damage, respectively. Finally, they studied the size of the damage zone in front of the mine opening with and without the presence of geological structures (Fig. 19). They concluded that the mining disturbance and presence of geological structures lead to a larger damaged area and thus increase the outburst risk and intensity. However, there was no time factor of mining reported in this work, and the progressive coal extraction leading to outburst initiation was not simulated, which is the key to mimicking the gradual build-up of an outburst condition.

The interaction of different physics proposed by Fan et al. (2017)

The area of outburst damage zones under different geological conditions (Fan et al. 2017)

More recently, Zhao et al. (2020) applied COMSOL Multiphysics to conduct outburst modelling by considering gas diffusion and the effect of deformation on coal matrix porosity (Fig. 20). They verified their model by comparing the model results with the gas production rate of a wellbore located in a coal mine in China. Then they studied the damage area as the representative of outburst-prone zones. The effects of geological structure, Young’s modulus, and overburden stress were also analysed using the proposed model (Fig. 21). In that paper, the matrix pressure gradient at the working face was more than the fracture pressure. The reason was mentioned to be the high permeability at the released stress region ahead of the excavation. But this seems to be unlikely given diffusion is known as a much slower process compared to Darcy’s flow. Also, other researchers state that one cannot expect the matrix pressure gradient to be more than the fracture pressure (Lu et al. 2019a, b). Furthermore, higher permeability results in more pressure decline, but in that study, the mass transfer between the fracture and matrix was assumed to be constant, so less pressure gradient was calculated for higher permeabilities.

The interaction of different physics used by Zhao et al. (2020)

The evolution of outburst damage zones under different mining conditions (Zhao et al. 2020)

Table 5 indicates the summary of fully coupled studies.

5.2 Iteratively coupled outburst models

In this method, the geomechanical and gas flow models were solved separately by either using different simulators for each part or using one simulator that supports both physics. The former is more common, in which the model results from each part need to be exchanged between two simulators over defined time steps.

Since 2000s, CSIRO in Australia has conducted a series of studies on outburst modelling. Choi and Wold (2001) applied the concept of dual-porosity gas reservoirs with multi-component coal seam gas into outburst modelling in coal mines. They used Darcy’s law to model the two-phase fluid flow in cleats. Also, gas desorption was modelled using Langmuir adsorption isotherms. They also used Fick’s law for gas diffusion modelling in the coal matrix, and so as to calculate the mass transfer between matrix and cleats. However, they did not study the change of matrix pressure and its effect on the mechanical behaviour of the model.

They used the in-house software packages SIMED and FLOMEC to couple geomechanics and fluid flow in outburst modelling, in which the stress-dependent permeability and Biot coefficient were considered. After that work in 2004, they further developed their model to simulate the outburst post-initiation stage by assuming coal dynamic failure with fragmentation, and considered a quasistatic yielding failure criterion. Parametric studies on outburst control parameters were also attempted, whereby gas content, geological structure, and coal strength were considered as the most effective parameters (Fig. 22) (Choi and Wold 2004).

Schematics of the coupling process in modelling gas outburst evolution (Choi and Wold 2004; Wold et al. 2008)

Xue et al. (2011) proposed a coupled geomechanical and fluid flow approach using two commercial software to model outburst initiation. Their study was limited to outburst initiation and did not cover outburst evolution or the mechanism of rock/coal fragmentation. The commercial geomechanical simulator FLAC3D was adopted to analyse mining-induced stress based on the classic finite difference method. They defined coal as a porous and fractured medium filled with gas and water, and later analysed its deformation using the concept of continuum mechanics, but no criterion for outburst damage evolution was proposed.

For the fluid flow part, the COMET3 code was used to simulate a dual-porosity single permeability system. Fluid flow was analysed in two categories, gas sorption and mass transport in the cleats and matrix. Adsorption/desorption of gas as a result of pressure change was studied using the Langmuir’s adsorption model with a monolayer assumption. For the mass transport, Darcy’s law was used for flow simulation in fractures, and Fick’s law was used for gas diffusion through the coal matrix. However, the matrix pressure change was not considered. In that study, the temperature was considered uniform and constant at the coal seam, and the flow between the matrix and fractures was assumed as pseudo steady state during the whole process.

Two coupling routines developed in C +  + and FISH language provided by FLAC3D were used to facilitate data exchange between these two software packages. The main function of the coupling routines is to transfer mutual parameters (permeability and pore pressure) between the two simulation packages (Fig. 23). In each interval, the effective stress was calculated in the geomechanical model and then the corresponding permeability was calculated based on a pre-defined stress-permeability relationship. This new permeability would be transferred to the fluid flow model. Using the basic equations of mass conservation for gas and water along with the Langmuir adsorption model, pore pressure would be calculated. New pore pressure was then returned to the geomechanical model to calculate the effective stress, and thus the iteration cycle would be repeated, and large deformation was considered as the outburst indicator (Fig. 24).

The coupling process between FLAC3D and COMET3 for outburst modelling (Xue et al. 2011)

Mesh deformation and pressure contour maps during 5 steps of roadway development (Xue et al. 2011)

Si (2015) attempted to develop outburst models to numerically replicate two prevalent outburst types, i.e., the dynamic failure outburst from floor seams and gas pocket outburst. The author applied outburst modelling to actual geological settings and explained the outburst cases that took place in field conditions. The geomechanical simulator FLAC3D together with the reservoir simulator ECLIPSE were used to model the geomechanical and fluid flow behaviour in gas outbursts, respectively. In that study, the criterion of outburst occurrence was assumed at the moment that the pore pressure exceeded coal cohesion in each element (Fig. 25).

Numerical modelling of outburst developed by Si (2015)

Lu et al. (2019a, b) developed a 3D outburst model using COMSOL Multiphysics. Although this software is capable of simultaneously solving both geomechanical and fluid flow equations (i.e., a fully coupled model), they used separated models for these two physics. They solved the solid mechanics' equations first under stationary assumptions, and then solved the time-dependent fluid flow equations. Therefore, their work is categorized as an iterative coupled model here. Three physics (solid mechanics, gas diffusion and gas flow through fractures) were considered, and the results of stress and pore pressure were studied (Fig. 26). A linear approximation with the equivalent plastic strain was adopted to estimate the post-peak permeability. The modelling results were then used to analyse the involved energy in an outburst event.

Stress distribution along the direction of coal extraction (Lu et al. 2019a, b)

Zhao et al. (2022) aimed to understand the response parameters during outburst events using thermos-fluid–solid coupling. They employed the Drucker-Prager criterion and the Mohr–Coulomb criterion, along with the maximum tensile stress criterion, to evaluate yield and tensile failure. For the fluid flow, they combined Darcy’s law, gas content, and the gas equation of state. The authors discovered that during outbursts, there was an increase in stress field change, a shift in abutment stress, and tensile failure (Fig. 27).

Stress distribution since before the outburst occurrence until the stress is stable conducted by Zhao et al. (2022)

Soleimani et al. (2023) used COMSOL Multiphysics and used the concept of conservation of energy to assess the possibility of the outburst occurrence during excavation through a coal seam and tried to find the key triggers. They assumed different mechanical properties for the different sections of the coal and conducted the calculations for different gas pressure values. The authors observed a relationship between the damaged zone and the intensity of the outburst occurrence which was indicated using energy balance concept (more explanation on Sect. 6.2). They monitored that although outburst occurrence is more likely to be observed during excavation close to the coal affected by the geological structures, higher gas pressure can cause the outburst in the earlier stages and be more intense (Fig. 28).

comparison between the outburst-prone region and damaged area conducted by Soleimani et al. (2023)

Table 6 indicates the summary of iteratively coupled outburst modelling studies.

As a dynamic multi-parameter phenomenon, the assessment of the coal and gas outburst is complex and requires techniques that cover all the key parameters and represent the reality of this incident. Numerical modelling is introduced as a method that can assist with that end. In this approach, solid mechanics and fluid flow simulation in coal seams are combined. The fluid flow should consist of both flow through the coal cleats and diffusion through the coal matrix. In the geomechanical model, the effect of excavation and geological structures should be studied.

Despite the previous efforts, no study considers all these parameters, and simplifications were used due to restrictions of software packages and running time. The most significant deficiency in previous models is the lack of comparison between field data and simulation results. Although the mathematical equations (coupling of solid mechanics and fluid flow, permeability, and damage equations) have been validated by limited experimental data, rare attempts were made regarding the validation of the outburst initiation model using field observations (Chen et al. 2016; Xu et al. 2006; Zhao et al. 2020). This is largely due to the difficulties in collecting first-hand field data on past outburst incidents. In addition, there is no attempt to relate the experimental method of outburst assessment (Sect. 3) with numerical modelling results. The experimental methods give precious information regarding gas ejection, but none of the previous studies used that information in numerical modelling or correlate it with modelling results. Therefore, there is no connection between the suggested indices and existing numerical models. By reviewing all the efforts that have been made, it seems that outbursts cannot be studied using just one method. A combination of different disciplines could lead to a better understanding of the phenomenon. In the following sections, we will explore some ideas that may open up new research perspectives.

6.1 Coal and gas outburst from particles point of view

In the models discussed above in 5, coal was studied as a continuum medium. Cundall and Strack (1979) first suggested using the discrete element method (DEM) in rock mechanics. Since coal and gas outburst is mostly observed in coal affected by geological structures in which coal has granulation properties, the properties such as porosity can be studied from the particles’ point of view (Guo et al. 2021). Furthermore, the coupling of the DEM and computational fluid dynamic (CFD) can be used as the method of coupling geomechanical rock failure and fluid flow, leading to more detailed outburst modelling, especially simulating the dispersion of released gas and coal particles in the ventilation airstream after outburst initiation. But unfortunately, this has not been widely applied so far. Using this DEM-CFD coupled approach, solid can be studied on the microscale as discrete particles, while fluid flow is studied on the mesoscale as a continuum medium. The drag forces causing the exchange of momentum and the effect of particle motion on fluid flow are used to couple the involved physics (Norouzi et al. 2016). One pioneering work was conducted by Xue et al. (2015), which used a DEM-LBM (Lattice Boltzmann method) coupled model to investigate all outburst stages. Their prediction for the outburst occurrence was based on the cavity shape by considering the formation of a large void with a small mouth as the outburst occurrence signature. The fact that the authors used a small-scale model to simulate a large-scale phenomenon (e.g., an outburst) needs to be carefully assessed. Thus, further efforts are required to use this method in field-scale studies and validate the coupling method using available experimental and field data.

Du et al. (2021) used the DEM-CFD coupling method to simulate the post-initiation stage of the outburst. The authors used EDEM and FLUENT software and coupled them explicitly. They simulated the transport of pulverized coal using the Hertz-Mindlin contact model in EDEM software, while FLUENT was used to simulate the pressure distribution and gas flow. The authors then compared the results with their test results from a previous work and found a 10% difference between the experiment and numerical simulation results (Wang et al. 2019a, b). However, that study only focused on the dynamic behaviour of the outburst and did not attempt to investigate the transfer from the quasi-static to the dynamic stage (outburst initiation), which is the most critical aspect of outburst prediction.

Another attempt for studying the outburst from the particles’ viewpoint was conducted by Ma et al. (2020) using the FEM with Lagrangian integration points (FEMLIP) and generating a 2D model with one-step excavation, as shown in Fig. 29. FEMLIP combines the benefits of both conventional Lagrangian FEM and Eulerian FEM by solving the distortion of meshes during deformation and tracking the geomaterials’ internal variables. In this study, the outburst initiation was defined as the transition of solid to fluid, recognized by the second-order work criteria introduced by Prunier et al. (2009) and Li et al. (2016) to describe the transition during the failure of granular materials. This definition is aligned with the 0.74 MPa pressure criteria, although it is worth noting that the pressure criteria may differ for coal with different properties, as found in Soleimani et al. (2023). The comparison between FEM, FEMLIP and DEM methods, their benefits and limitations are mentioned in Table 7.

Outburst modelling using FEMLIP by Ma et al. (2020)

In the previously mentioned studies, particles were assumed to be spherical, even though the particles' shape could affect the occurrence of outbursts (as explained in Sect. 6.2). This assumption adds several extra points that need to be considered, such as the reconsideration of contact forces and detection, particle rotation, and friction forces.

6.2 Conservation of energy during coal and gas outburst

The previous numerical models are generally used to describe the pre-initiation stage of the outburst occurrence, but outburst initiation is known as the moment that the quasi-static deformation becomes dynamic (Choi and Wold 2004). Due to the large variety of involved parameters and complexity of the outburst event, it has been suggested to analyse outburst occurrence from the energy point of view and use the energy conservation law to explain the process (Cai and Xiong 2005; Gray 1980; Hodot 1966; Jiang and Yu 1996; Li et al. 2012; Valliappan and Zhang 1999; Wen et al. 2002; Zhao et al. 2016).

During outburst initiation, two groups of energy are involved, supplier energy and consumer energy. The supplier energy is the combination of elastic energy and the expansion energy of desorbed and free gas. This energy group should overcome the required energy to break coal, in addition to supplying the minimum kinetic energy to move coal particles. Also, inevitably, there would be some energy loss during this process (Lu et al. 2019a, b; Nie et al. 2019; Wu et al. 2020).

Lu et al. (2019a, b) combined the energy analysis with the numerical modelling, but they only went through the supplier energy calculation. For the initiation moment, they used the time of recorded sound during an outburst event in a Chinese coal mine as the outburst initiation time in their energy analysis. But they did not use the coal and gas properties from the same coal mine, which may lead to inaccurate results. Furthermore, In the thermos-fluid–solid coupling model, Zhao et al. (2022) also computed the elastic energy as a function of stress and observed that tensile failure occurred during outbursts aligned with a 3–5 m shift in the stress concentration ahead of the mining face. Soleimani et al. (2023) conducted a numerical model and used the results to calculate both supplier and consumer energies. The authors defined the exceedance of supplier energies from the consumer energies as the moment of outburst occurrence and found reliable results compared to other research, such as Xue et al. (2021) work in which supplier energies for six coal seams with the reported outburst incidents were calculated. The deficiency of this work was using semi-empirical and laboratory data (Nie et al. 2019) to estimate consumer energies. Guo et al. (2021) used the concept of energy conservation to calculate the elastic, rotation kinetic and translational kinetic energies and the consumption energies (i.e., the external work) from the particles’ viewpoint.

In all the studies on the outburst considering the coal as an aggregation of particles, they were assumed granular. Luo et al. (2018) studied the effect of particle shape and the total surface area of particles after crushing by defining the equivalent diameter on the crushing energy involved in the outburst occurrence. They broke down coal samples into a set of uniformly shaped and sized particles and considered the sieve diameter through which these particles could pass as the equivalent diameter of that group. They analysed the crushing energy for fine particles (with a mean equivalent diameter of less than 0.075 mm) and concluded that coal samples from tectonically affected zones (known as outburst-prone coal) were easier to be crushed than coal in geologically normal zones. In their experiments, by applying the same crushing energy, outburst-prone samples were converted to much smaller particles.

Also, they analysed the equivalent diameter for different particles’ shapes and calculated the newly formed surface areas based on these diameters (Table 8). These results revealed that different particle shapes can change the total surface area to up to 5 times the spherical-shaped particles. Figure 30 shows that, at the same surface area, the particle shape can effectively change the crushing energy.

The effect of coal particle shapes on the crushing energy (Luo et al. 2018)

6.3 Statistics in coal and gas outburst

Stochastic modelling can improve the accuracy of coal and gas outburst assessment from two perspectives: the estimation of properties within the coal seam and the assessment of the likelihood of outburst occurrence based on different values of involved parameters.

The first perspective is useful in applying the effect of heterogeneity and conducting more realistic modelling. Wold et al. (2008) made the first attempt in this area. They applied statistical distributions on rock/coal properties to monitor the effect of heterogeneity on outburst occurrence. The researchers found that the gamma probability distribution was suitable for fitting the core permeability data measured in the lab. This led them to use the gamma distribution to define permeability input for their reservoir models. Later, using permeability measurement data from field well-testing, they defined another gamma distribution for permeability input by rescaling the core data. They also defined a gamma distribution for strength data measured by Rock Strength Device (RSD) tests. They then used the Monte Carlo sampling method to define the spatial distribution of permeability and strength across the coal seam (as shown in Fig. 31). The researchers concluded that the variation of these two parameters can significantly affect outburst occurrence and suggested that this approach is an effective method for assessing the variability in outburst analysis.

Comparison of the pore pressure distribution when using well-testing permeability input and core permeability measurement input (Wold et al. 2008)

Later on, Zhao et al. (2020) used the same concept to study the heterogeneous effect of Young’s modulus (Fig. 32). They defined a probability density function using the Weibull distribution (Eq. (24)) (Fan et al. 2018; Zhu and Tang 2004; Zhu and Wei 2011). \(\overline{E }\) represents the average value of Young’s modulus, and \({m}_{\text{h}}\) is the homogeneity coefficient, which means by increasing its value the homogeneity of coal will increase. Their study also showed that heterogeneity can affect the distribution of stress, damage evolution, gas pressure and coal permeability (Fig. 33).

Distribution of coal elastic modulus when \(\overline{E }\)=2.7 GPa and \({m}_{\text{h}}\)=6 (Zhao et al. 2020)

Effect of heterogeneity on the damage evolution induced by the outburst (Zhao et al. 2020)

Besides the abovementioned studies using probability distribution to assign properties to the coal, geostatistics is another method that could be used for that purpose. Geostatistical methods, particularly Kriging, are commonly used to estimate properties in unsampled locations. Kriging can be performed through three approaches: ordinary, simple, and universal Kriging. These methods can accurately estimate the value of the required parameter based on the number of sampled locations and the distance from those locations (Srivastava 2013).

As mentioned, statistics can be applied for the prediction of outburst occurrence other than properties estimation. The assessment of the possibility of the outburst occurrence using statistical methods has been conducted through machine learning algorithms. The advancement of Artificial Intelligence (AI) and Machine Learning techniques has led researchers to conduct studies using these methods for assessing the risks of mining hazards, such as coal and gas outburst events. Both classification and regression methods have been utilised. The classification methods examine whether a mine is prone to outbursts or not, while the regression methods estimate values such as emitted gas volume (Anani et al. 2022). The following paragraphs summarize some of these algorithms used in the coal and gas outburst.

An attempt towards using regression methods was conducted by Bodlak et al. (2018) who used the random forest method, which is an ensemble technique that combines multiple decision trees to create a more robust prediction model. They employed the XGBoost algorithm, a variant of gradient boosting that is designed to handle large datasets and deliver high accuracy in a computationally efficient manner. The researchers used this approach to predict the mass of rock, the volume of gas, and the length of collapsed and/or damaged workings. They employed cross-validation with RMSE as the scoring method, and they verified the predictions using out-of-sample data from the same region.

Comparison of different algorithms to classify the possibility of the outburst occurrence was done by Xie et al. (2019). They conducted a study on the IFOA-GRNN algorithm, a combination of the Improved Fruit Fly Optimization Algorithm (IFOA) and the General Regression Neural Network (GRNN) algorithm to optimise the model’s parameters and make predictions at the same time. FOA is a swarm intelligence-based optimisation algorithm that uses the behaviour of fruit flies to search for the optimal solution to a problem. The movement of each fly is influenced by the attraction to food (the optimal solution) and the repulsion from other flies. IFOA is an improved version of FOA, where additional techniques are used to improve the performance of the algorithm. The main difference between FOA and IFOA is that IFOA incorporates additional techniques, such as dynamic adjustment of the control parameters to improve the efficiency and accuracy of the optimisation process. They compared the results with GRNN, FOA-GRNN, and Back Propagation (BP), which is a supervised learning algorithm used to train artificial neural networks by adjusting the weights in the network to minimise error. They defined five categories of no-outburst, small, medium, large, and extra-large outbursts. Their model consisted of four levels (geological factors, coal structure factors, gas factors, and operation factors) with a total number of fourteen factors. They compared the results of four models using Fitting degree (R²), Average Error (AE), Mean Squared Error (MSE), and Risk Level Prediction Accuracy (RLPA) methods, and concluded that IFOA-GRNN is the most accurate one.

Later on, Xie et al. (2022) conducted a sensitivity analysis on the results of the Apriori-Bayesian algorithm to prevent coal and gas outburst accidents by analysing data related to past accidents and deducing accident causes (Fig. 34). The Bayesian network was modelled by Netica software to find the causal relationships between the influential factors. In this type of network, nodes represent variables, and the edges represent the causal relationships between them. Nodes can be categorised as parent, child, or intermediate nodes, depending on their position in the network and their relationships with other nodes. Parent nodes directly influence the values of child nodes, while intermediate nodes are influenced by parent nodes and, in turn, influence other nodes in the network. The Apriori algorithm was employed to identify the reasons behind coal and gas outburst accidents and to establish a map of the influence paths between these accident causes. In that paper, the outburst phenomenon was studied within five aspects: unsafe conditions, unsafe acts, individual safety ability, safety management system, and safety culture, by considering 84 outburst incidents in China.

Apriori-Bayesian algorithm used by Xie et al. (2022) to prevent coal and gas outburst

The main issue with the usage of machine learning methods in the coal and gas outburst is lack of the data and proper validation, which needs to use new data different from the samples used to train the algorithm. However, as the safety standard is improving over the years and less outburst events can be expected in the future, it is difficult to obtain more field data for outburst research given the limited database (Anani et al. 2022; Xie et al. 2022). A suggestion to solve this problem is to use machine learning in conjunction with numerical modelling to generate more possible scenarios causing outburst incidents.

Coal and gas outburst is a serious safety issue, and the propensity of the coal mines toward its occurrence should be studied, and it is important to find out about the triggers and the outburst initiation process. Due to the complexity of its nature, and considering the power of computers and the wide range of software packages computational methods could be the best available method to cover the most effective parameters. To that end a coupled geomechanical and fluid flow model is necessary and most of the efforts so far have used the continuum mechanics and Darcy’s law, few attempts have been done to study the phenomenon based on an aggregation of particles’ behaviour. Even in those methods, the assumption was spherical particles, whereas their shape could make difference in the required energy to crush the coal. Furthermore, the outburst could be studied from the conservation of energy and stochastic modelling viewpoints. Due to the complexity of the coal and gas outburst’s nature, it seems that the combination of methods could lead into a better assessment of the hazard.