Rockburst and Gas Outburst Forecasting using a Probabilistic Risk Assessment Framework in Longwall Top Coal Caving Faces

A probabilistic risk assessment framework was developed to mathematically represent the complex engineering phenomena of rock bursts and gas outbursts for a heterogeneous coal seam. An innovative object-based non-conditional simulation approach was used to distribute lithological heterogeneity present in the coal seam to respect their geological origin. The changing mining conditions during longwall top coal caving mining (LTCC) were extracted from a coupled numerical model to provide statistically sufficient data for probabilistic analysis. The complex interdependencies among abutment stress, pore pressure, the volume of total gas emission and incremental energy release rate, their stochastic variations and uncertainty were realistically implemented in the GoldSim software, and 100,000 equally likely scenarios were simulated using the Monte Carlo method to determine the probability of rock bursts and gas outbursts. The results obtained from the analysis incorporate the variability in mechanical, elastic and reservoir properties of coal due to lithological heterogeneity and result in the probability of the occurrence of rock bursts, coal and gas outbursts, and safe mining conditions. The framework realistically represents the complex mining environment, is resilient and results are reliable. The framework is generic and can be suitably modified to be used in different underground mining scenarios, overcoming the limitations of earlier empirical indices used. Parameters causing rockbursts and gas outbursts were linked along with their influences and interdependencies into a probabilistic risk assessment framework. Dynamically updated system feedback from the numerical model was fed into the framework to represent the current stress state in retreating mining and estimate the probability of the occurrence of rockbursts and gas outbursts. Statistically significant data were used to quantify the probability of rockbursts and gas outbursts using Monte Carlo simulation. Parameters causing rockbursts and gas outbursts were linked along with their influences and interdependencies into a probabilistic risk assessment framework. Dynamically updated system feedback from the numerical model was fed into the framework to represent the current stress state in retreating mining and estimate the probability of the occurrence of rockbursts and gas outbursts. Statistically significant data were used to quantify the probability of rockbursts and gas outbursts using Monte Carlo simulation.


ERR n
Energy released at nth excavation step)J/m 3 ) G Shear modulus (GPa) h block

Height of the coal block (m) k h
Horizontal permeability (m 2 ) k n Permeability at nth excavation step (m 2 ) k v Vertical permeability (m 2 ) K Bulk modulus (GPa) l block Length of the coal block (m) m f Mass of ejected coal (kg) p g Gas pressure in the coal matrix (MPa) p n−1 Gas pressure at n−1th excavation step (MPa) p n Gas pressure at nth excavation step (MPa)

P atm
Atmospheric pressure (Pa) P L Langmuir pressure (MPa) s New surface area (m 2 /g) t n−1 Time at n-1th excavation step t n Time at nth excavation step V a The volume of adsorbed gas (m 3 ) V c The volume of coal ejected (m 3 ) V ej Ejection velocity (m/s) V free The volume of free gas (m 3 ) V gas Volume of gas in coal (m 3 ) V L Langmuir volume (m 3 /tonne) V m Volume of coal (m 3 ) w Specific energy of coal (J/m 2 ) w block Width of the coal block (m) W adsorbed Adsorbed gas energy (J/m 3 ) W c Crushing energy (J/m 3 ) W free Free gas energy (J/m 3 ) W gas Total gas expansion energy (J/m 3 ) W k Transportation energy (J/m

Introduction
Rockbursts and gas outbursts are complex engineering phenomena having non-linear dependencies on several parameters.According to several multi-factor theories proposed and widely acknowledged by researchers, rock bursts and gas outbursts result from the combined effects of stress, gas pressure, and mechanical properties of coal (Wang and Xue 2018;Cheng et al. 2021).They are special hazards in terms of their suddenness, frequent occurrences, and high consequences (Liu et al. 2016;Zhou et al. 2016c;Agrawal et al. 2021).As per the most widely accepted energy hypothesis for rockburst and gas outburst occurrences, the strain energy combined with the gas expansion energy causes violent failure (Cao et al. 2019a;Dai et al. 2019;Tu et al. 2019Tu et al. , 2021;;Lei et al. 2020Lei et al. , 2021;;Cheng et al. 2021;Xue et al. 2021).Despite several years of research, the complex rockburst and gas outburst mechanisms remain not well understood, making them difficult to forecast accurately (Zhou et al. 2020).
Hazard forecasting can be classified into long-term and short term (Pu et al. 2019).Long-term forecasts are useful during the planning stage and early stages of development when limited data are available, to assess the feasibility of the project, make suitable design changes as well as guide future mining operations (Li et al. 2019).These methods include the use of empirical indices, traditional risk assessment approaches, and predictive models to ascertain the likelihood of hazards at excavation sites.Short-term forecasting predicts the time, location, and intensity of damage as the mining progresses using continuous and reliable field monitoring data (Pu et al. 2019;Dou et al. 2014;Liu et al. 2016;Tang et al. 2016;Yin et al. 2016).Short-term forecasting is data-intensive and needs continuously updated values to reliably monitor the dynamically changing mining conditions.Research presented here focuses on long-term forecasting.
Empirical indices were mostly developed based on the analysis of coal/rock behaviour to different confinement, gas pressure, and stress conditions during loading/unloading experiments in a laboratory environment.These indices have several limitations in terms of practical application in the field and cannot be used in complex and dynamically changing environments to realistically forecast hazard potential (Li et al. 2019).Different indices may have different critical values that depend on the mining conditions and vary from one mine to another.These indices may identify a working face to be safe but accidents have occurred at such safe faces, suggesting that they may provide a false sense of safety (Tang et al. 2016).
Traditional risk assessment approaches include analytical hierarchical process, bowtie diagrams, cause-consequence analysis, decision analysis, event-tree analysis, fault-tree analysis, interval analysis, multi-risk analysis, etc. (Cosgrove and Hudson 2016).These approaches are based on the analysis of event initiation and sequences leading to system failure and calculate the probability of different failure outcomes using logic gates.They are not suitable for dynamic systems, nor for those that have complex non-linear interdependencies among failure processes as they may need significant pre-processing and still may not represent a realistic scenario (Goldsim 2017).They also suffer from a lack of response or feedback from the system, which can be detrimental in the case of rapidly evolving systems (Mattenberger et al. 2015;Goldsim 2017).In addition, traditional risk assessment approaches do not provide an immediate actionable result (Mattenberger et al. 2015) and suffer from having statistically insufficient data (Maleki 1995).
The development of predictive coupled numerical models, representing realistic conditions leading to rockbursts and gas outbursts in coal mining, provides statistically sufficient data overcoming the limitations of traditional risk assessment approaches.Researchers have proposed several forecasting indices based on numerical simulation analysis.Using the energy analysis, Salamon (1984) proposed an energy release rate ( ERR ).Mitri et al. (1999) used mining- induced strain energy density and energy storage index to calculate a burst potential index.Vieira and Durrheim (2002) advanced the concept of ERR to put forward incremental energy release rate ( ERR i ) during mining.Beck and Brady (2002) and Wiles (2006) proposed the concept of local energy release density ( LERD ) to reflect total energy change during the pillar failure process.Wiles (2002) proposed the local energy release rate ( LERR ) and Heal et al. (2006) pro- posed rockburst vulnerability index ( RVI ) to forecast rock- burst occurrences.
The numerical modelling-based indices should be used with caution as these were mostly developed for a site-specific scenario and do not incorporate the spatial variability and complexity involved in rockbursts and gas outbursts.Furthermore, most of the indices were developed based on static simulation results and cannot reasonably reflect the dynamic process of actual events.The geological conditions, rock properties, and mining design vary in different parts of the world, so the threshold values for these indices should be suitably modified to suit the geological conditions of the mining area.Mitri et al. (1999) found that the ERR values in Canadian mines were lower than that for South African mines experiencing rockburst problems.Most of these indices have ignored the heterogeneity occurring in geological formations or have over-simplified its incorporation using normal or Weibull distribution (Wang et al. 2017(Wang et al. , 2019)).
Several researchers have used these monitoring data to predict rockburst and gas outburst occurrences (Zhang et al. 2014(Zhang et al. , 2013;;Tang et al. 2016;Zhou et al. 2012;Mutke et al. 2015).Jia et al. (2015) combined acoustic emission, electromagnetic, and microseismic data to propose a multiagent approach for rockbursts prediction.Li et al. (2016) integrated electromagnetic radiation and microseismic data to predict rockbursts.Si et al. (2015a) combined microseismic and seismic tomography monitoring to identify excessive gas emissions occurring due to geological anomalies.These methods are characterised by their simple, low labour cost and rapid operation.However, they suffer from large prediction errors due to lithological heterogeneity in coal seams, the anti-interference performance of the signals and the precision of data identification (Tang et al. 2016;Pu et al. 2019).
These intelligent methods establish models using existing data (training datasets) to analyse the weight ratio of several factors for their quantitative calculations.With the change in excavation areas with different mining conditions, the weight ratio of various factors used to forecast rockbursts and gas outbursts may change (testing datasets).The lack of original data to represent the changed mining conditions makes determining the new weight ratios difficult.Furthermore, when such updated information is lacking, predictions using such algorithms and mathematical formulations may fail.Intelligent methods present difficulties such that they are time-consuming, difficult to interpret, computationally extensive, and suffer from overfitting and lack of transparency.These models are highly dependent on the training dataset quality and their outputs may be mutative (Zhou et al. 2020).
To overcome the limitations of approaches developed earlier, and indices, this research developed coupled numerical models of a longwall top coal caving (LTCC) face by distributing different coal lithotypes using an object-based non-conditional simulation approach to respect the geological origin of heterogeneity present in the coal seam at Coal Mine Velenje.The data representing the stochastic behaviour of abutment stress, pore pressure, the volume of total gas emission and incremental energy release rate were extracted from the coupled numerical model and fed into a probabilistic risk assessment (PRA) framework to accommodate the uncertainty and variability occurring at each excavation step in the dynamically changing mining environment.This paper describes a PRA framework developed in GoldSim focused on risk occurrence rather than consequences (e.g.cost or other implications).The PRA framework incorporates Monte Carlo simulation to forecast the occurrence of rockbursts and gas outbursts and uses a retreating LTCC mining panel representing conditions prevalent at Coal Mine Velenje to illustrate its implementation.

Coal Mine Velenje
The Velenje coal basin lies in a synclinal valley bounded by two major faults, the Sostanj and Smrekovec faults (Fig. 1).The lithological sequence was controlled by the movement of the Periadriatic fault system causing a high subsidence rate.The simultaneous deposition of paleo-forests, bush swamps, dwarf plants and vegetation by open water formed the heterogeneous clastic sediments more than 1000 m thick in the basin.Forests compacted to a lesser degree as compared to dwarf plants at the same time under different depositional environments leading to the formation of different coal lithotypes with a varying abundance of xylites, detrites, and mineral matters (mostly composed of alumino-silicates and carbonate minerals) throughout the lignite deposit (Markic and Sachsenhofer 2010).Differential compaction led to the bowl shape of the coal seam as illustrated in Fig. 2a.The seam extends along WNW-ESE direction and is approximately 8.3 km long and 1.5-2.5 km wide, having a thickness of up to 165 m at the centre and pinches out towards the margins (Kanduc and Pezdic 2005;Likar et al. 2012;Si et al. 2015b;Durucan et al. 2019).
The coal seam is mined using a combination of multilevel mining and LTCC mining known as the Velenje Mining Method (Jeromel et al. 2010;Likar et al. 2012).The coal deposit is divided into 10-20 m mining levels depending on the coal seam thickness (Fig. 2b).The bottom 3-4 m is cut by a shearer under hydraulic power supports and the 7-17 m Fig. 1 Geology of the Velenje coal basin and surrounding area (after Brezigar 1986) top coal is caved and recovered in front of the supports to allow steady face advance (Fig. 2c) (Si et al. 2015b;Cao et al. 2018;Durucan et al. 2019).Seam gas at Coal Mine Velenje is a mixture of CO 2 and CH 4 varying from 98% CO 2 -2% CH 4 to 0% CO 2 -100% CH 4 in different parts of the mine.

Coal Seam Lithological Heterogeneity
To represent lithological heterogeneity in the coal seam, a three-dimensional fork-shaped solid was created with varying extent along the X-, Y-and Z-axis and lateral offset along the Z-axis to represent the deposition by river streams over a long time to form an area with different xylite distributions for use in the numerical models (Fig. 3).The geological setting dictates the use of an object variable (ellipsoids) to represent the lithological heterogeneity present in the coal seam at Coal Mine Velenje.An object-based non-conditional simulation was used to randomly drop ellipsoids in the solid structure to simulate different lithotypes.The shape of the ellipsoids was varied by changing the ratio of semi-major axes along the X-, Y-and Z-axes.The length of the ellipsoids was varied using a Gaussian distribution.The Poisson point process was used to independently and randomly distribute ellipsoids in the model (Chiles and Delfiner 2012).A random dose was given to each ellipsoid and the summation rule was applied to monitor the locations where ellipsoids overlap each other.The overlapping locations were considered as xylite, while locations where no ellipsoids were present were considered as detrite.The random distribution facilitates mimicking different orientations, shapes and sizes of xylite present due to different depositional environments that were found by petrographic analysis by Markic and Sachsenhofer (2010).

Energy Analysis
The energy hypothesis based on the strain energy, free gas energy, gas expansion energy, crushing energy and transportation energy was used to determine triggering conditions to forecast rockbursts and gas outbursts in the PRA framework.

Strain Energy
The energy release rate ( ERR ) proposed by Salamon (1984) has been used by several researchers to evaluate rockburst potential in deep underground mining (Zhang et al. 2017).Originally developed for hard rock mines, the index has been widely used in coal mining in the US (Sears and Heasley 2009).Salamon (1984) noted a good correlation between ERR and rockburst occurrences.ERR can be calculated as (Vieira and Durrheim 2002), where ERR is the energy release rate (J/m 3 ), zz is the verti- cal stress acting on the solid face before it is mined (MPa), and d z is the face convergence due to mining (mm).
(1) Vieira and Durrheim (2002) proposed a probabilistic methodology to calculate rockburst risk by applying a simple technique to remove multiple counts of ERR at the same position by introducing incremental ERR ( ERR i ).ERR i is defined as the energy released between excavation step n − 1 ( ERR n−1 ) and excavation step n ( ERR n ).It is given as, The ERR i index does not consider initial convergence and in situ energy that can be stored in the rock before extraction starts.Vieira and Durrheim (2002) considered rock to be elastic except at the mining face and simulated a hypothetical case.In the research presented in this paper, a strain-softening model was used to represent coal failure, the model was equilibrated before excavation to account for in situ energy, and a retreating LTCC mining was modelled as practised at Coal Mine Velenje (Si et al. 2015b).

Gas Energy
The gas expansion energy present in coal is related to the free gas and adiabatic expansion of gas desorbed from the coal structure.The volume of free gas present in the coal seam can be calculated as (Tu et al. 2019(Tu et al. , 2021)), where V f ree is the volume of free gas present in the cracks and fractures (m 3 ), V m ( = l block × w block × h block ) is the vol- ume of coal ( l block is the length of the coal block (m), w block is the width of the coal block (m), h block is the height of the coal block (m)) (m 3 ), is the coal porosity (%), p g is the gas pressure in the coal matrix (MPa), P atm is the atmospheric pressure (MPa) and is the ratio of specific heats at constant pressure over constant volume for methane and carbon dioxide ( = 1.31).
The volume of adsorbed gas can then be calculated as, where V a is the volume of gas adsorbed in the coal (m 3 ) and V gas is the volume of gas released (m 3 ).
Free and adsorbed gas energy can be calculated as (Tu et al. 2019(Tu et al. , 2022;;Lei et al. 2020;Xue et al. 2021), (2) where W f ree is the free gas energy (J/m 3 ), W adsorbed is the adsorbed gas energy (J/m 3 ), W gas is the total gas expansion energy (J/m 3 ), V a is the total volume of gas adsorbed at n th excavation step (m 3 ), V m is the volume of coal ejected (m 3 ), and p n is the gas pressure in the coal seam at n th excavation step (MPa).

Crushing Energy
Several physical processes like abrasion, crack propagation, crushing, fracture, friction, and vibration require energy.The energy dissipation during coal crushing includes heat energy and sound energy losses (Lei et al. 2020).The major energy dissipation occurs in crushing the coal which is directly proportional to the new surface area of the pulverised coal (Dai et al. 2019;Cheng et al. 2021).Due to the presence of several cracks, fissures, and fractures in the coal before an outburst, it is difficult to identify an exact increase in the surface area.The crushing energy required can be calculated as (Cheng et al. 2021), where W c is the crushing energy required per unit volume of coal (J/m 3 ), s is the new surface area (m 2 /g), w is the specific energy of coal (J/m 2 ), and is the density of coal (kg/m 3 ).

Transportation Energy
As a rockburst or gas outburst propagates, broken coal particles move from their original positions into the mine openings (Xue et al. 2021).Researchers proposed a horizontal projectile (Lei et al. 2020) or parabolic motion (Cao et al. 2019a) of the broken coal particles.Previous researchers suggested an ejection velocity in the range of 8-50 m/s for a violent failure to occur (Hosseini et al. 2015;Zhou et al. 2016c;Gale 2018).Using the kinetic energy equivalent, the (5) minimum energy required for violent ejection of coal/rock into the mine workings can be calculated as, where W k is the transportation energy required for violent ejection (J/m 3 ), m f (= V m × ) is the mass of the coal/rock block being ejected (kg), and V ej is the ejection velocity of the coal/rock mass (m/s).

Conditions for Rockbursts and Gas Outbursts
Unconfined compressive strength, tensile strength and total energy stored per unit volume of coal/rock are good indicators to forecast rockbursts and gas outbursts (Nussbaumer 2000).For rockbursts and gas outbursts to occur, several conditions need to be simultaneously satisfied, these are: strength conditions, energy conditions and rockburst/gas outburst tendency (Dou et al. 2018;Cai et al. 2019).

Strength Conditions
The essential strength condition implies that rockbursts may occur when the stress acting on the coal at nth excavation step exceeds the uniaxial compressive strength (  n >  c ) such that the coal has failed (Wang and Park 2001;Cai et al. 2016;Dou et al. 2018).Similarly, for gas outbursts to occur, the difference between the gas pressure at the exposed face to the atmospheric pressure should be more than the tensile strength of coal (Lei et al. 2021;Tu et al. 2021).

Energy Conditions
The minimum energy accumulated in the coal seam should overcome the required crushing energy and additional energy losses to cause rock failure (Wen et al. 2016;Dou et al. 2018;Canbulat et al. 2019).The accumulated energy should be sufficiently higher than the required transportation energy to cause the violent ejection of broken coal particles (Cai et al. 2016;Fedotova et al. 2017;Canbulat et al. 2019).

Bursting Tendency
The ability of a coal/rock to store strain energy and release them instantaneously when a failure occurs is determined by the brittleness index of coal ( c / t ) (Peng et al. 1996;Xu et al. 2020).The energy accumulated in the coal due to retreat longwall mining at n th excavation step ( ERR i ) and free gas energy ( W f ree ) should be more than the minimum crushing energy ( W c ) and the ( 9) transportation energy ( W k ) (Wen et al. 2016).Based on the energy hypothesis, different triggering criteria were proposed, Rockbursts

Gas outbursts
Coal and gas outbursts If the above triggering criteria are not met, the energy stored is insufficient to eject the broken coal violently, then rockbursts and coal and gas outbursts are unlikely.For gas outbursts, the free gas energy should be more than the crushing energy required as adsorbed gas takes time to desorb from the coal surface to flow into the cracks.If the free gas energy does not overcome the crushing energy required, it inhibits rapid desorption, and thus gas outbursts may not occur.In such a scenario, if the strength conditions are met, it may lead to quasi-static coal failure and excessive gas emissions. (10

Coupled Geomechanical and Gas Flow Modelling
A retreating LTCC mining model of 500 m × 500 m × 150 m dimension with a grid size of 5 m × 5 m × 3 m along the X-, Y-, and Z-axes, respectively, was constructed in FLAC 3D .The X-axis was set along the length of the longwall panel, Y-axis was set across the width of the longwall panel, and Z-axis was set along the vertical direction passing through the origin at the bottom of the model.Different rock layers were simulated to represent the lithology at Coal Mine Velenje.The lithology considered for modelling had a 15 m mining level, overlain by 60 m goaf of previous mining levels and 15 m clay, and underlain by 60 m coal on the floor (Fig. 4).
The fork-shaped solid representing the heterogeneous zone was inserted in the FLAC 3D model of the coal seam, assuming full penetration to generate 10,703 heterogeneous zone elements spanning the mining level and underlying coal seam (Fig. 5a).Detrite and xylite zones were assigned depending on the heterogeneity distribution (Fig. 5b).The longwall face modelled cuts through the heterogeneous zone having different xylite blocks at each excavation step.A representative cross-section along XX (Fig. 5c) and YY (Fig. 5d) shows the interaction of heterogeneous zones with the longwall panel.
The baseline geomechanical properties of different formations at Coal Mine Velenje, taken from experimental work carried out by the former Velenje engineers and published literature, were assigned to the model (Table 1).Mohr-Coulomb strain-softening model was implemented in FLAC 3D to represent the mechanical behaviour of coal and coal-measure rocks (Itasca 2017).The model was constrained with stress boundaries on both sides (X-and Y-axes), and the bottom (Z-axis) was fixed.The model was placed below − 305 m, with an overburden density of 2360 kg/m 3 initialised in the model.The load corresponding to 305 m (7.19 MPa) thick overburden was applied on the top surface.The model was gravity loaded and reached initial equilibrium before excavation started.
In the coupled model, the coalbed methane module of ECLIPSE 300 with two different coal regions was used to represent different gas adsorption behaviour for detrite and xylite (Schlumberger 2017).Several sub-routines were written in FLAC 3D to facilitate the exact assignment of reservoir properties, gas pressure gradient, and Langmuir properties to each lithotype present in the model (Table 2 and Fig. 6a).Si et al. (2015b) monitored the gas pressure dynamics at different mining levels and found that the gas pressure was maximum in the first mining level but due to   K is the bulk modulus (GPa), G is the shear modulus (GPa), c is the cohesion (MPa), is the angle of internal friction ( o ), t is the tensile strength (MPa),c r is the residual cohesion (MPa), and tr is the residual tensile strength (MPa) subsequent over-mining, the seam gas pressure dropped to around 0.8 MPa in the modelled mining level.They also determined the underlying seam gas pressure referring to borehole data and found it to be stable at around 1.6 MPa before extraction commenced.
In the coupled numerical simulation workflow, the pore pressure ( p n−1 ) calculated at excavation step ( n − 1 ) by ECLIPSE 300 was passed to FLAC 3D to calculate a provisional stress state ( n ′ ) at the excavation step ( n ).The provi- sional stress state ( n ′ ) and pore pressure ( p n−1 ) was used to calculate the permeability ( k n ) at the excavation step ( n ) in FLAC 3D .The permeability ( k n ) was passed to ECLIPSE 300 to calculate the updated pore pressure ( p n ).The simulated fluid time, t n−1 to t n , represent coal extraction time at the excavation step ( n ).The updated pore pressure ( p n ) was fed back to FLAC 3D to re-equilibrate and calculate the actual stress state ( n ) at the excavation step ( n ).The actual stress state ( n ) and pore pressure ( p n ) was used as an input for the next excavation step ( n + 1 ).Further details of the coupled model can be found in Si et al. (2015c).
To represent LTCC mining with a face advance rate of 5 m in each excavation step, first, the lower 3 m thickness of the coal was extracted (nulled) to represent cutting by the longwall shearer and equilibrated to create the initial face; next, the top coal (12 m) left during the previous step was instantaneously removed; and finally, the complete 15 m mining level was reinstated with an elastic goaf material property and solved to equilibrium to represent one excavation step in LTCC (Fig. 5b).The process was repeated for several excavation steps, advancing the longwall face through the heterogeneous zone with a varying abundance of xylite, to represent soft conditions and hard conditions of retreating longwall mining.

Probabilistic Risk Assessment Framework
In probabilistic risk assessment (PRA), any complex engineering system that can be described quantitatively using mathematical equations can be modelled (Mattenberger et al. 2015;Song and Yang 2018).Fixed as well as stochastic m is the matrix porosity (%), c is the cleat porosity (%), k h is the horizontal permeability (m 2 ), k v is the vertical permeability (m 2 ), P L(CH 4 ) is the Langmuir pressure for CH 4 (MPa), P L(CO 2 ) is the Langmuir pressure for CO 2 (MPa), V L(CH 4 ) is the Langmuir volume for CH 4 (m 3 /tonne), and V L(CO 2 ) is the Langmuir volume for CO 2 (m 3 /tonne) parameters can be implemented in the framework to represent uncertainty.The Monte Carlo simulation approach is used to propagate the uncertainty in input throughout the system and the output (Mattenberger et al. 2015;Goldsim 2017;Song and Yang 2018).Multiple triggering criteria can be defined and monitored to analyse the system's response to the stochastic variations.The PRA calculates the probability of the occurrence of unlikely but high-consequence outcomes as the proportion of realisations where the triggering criteria were met (Goldsim 2017).The results obtained from PRA analysis are quantitative, reliable and provide a meaningful alternative to traditional risk assessment approaches (Song and Yang 2018).
In this research, a PRA framework was developed to estimate the probability of rockbursts and gas outbursts in retreating LTCC mining (Fig. 7).The fixed parameters considered include dimensions of the longwall panel (the length, width, and height of the coal block mined), physical properties of coal (density, porosity, cohesion, the angle of internal friction, and specific energy), and initial conditions and parameters of gas migration (the initial gas pressure in the coal matrix, atmospheric pressure, and adiabatic constant).These values do not change in a longwall panel and are mostly fixed in the evaluation.To incorporate the dynamic changes and system feedback during mining, the vertical stress acting on coal, incremental energy release rate, and total volume of gas emission at the face for each excavation step were fed into the model as stochastic variables.Pore pressure did not show much variation and thus a fixed value was fed into the model.Fixed and stochastic parameter values were used to estimate several secondary parameters that were needed to calculate energy and stress states using Eqs.(1)-( 9).Several triggering criteria that control the occurrence of rockbursts and gas outbursts (Eqs.( 10)-( 12)) were assigned in the framework.The probability of a rockburst occurrence was calculated without taking gas expansion energy contribution into account (Eq.10).The probability of a gas outburst was calculated without considering the strain energy contribution (Eq.11) and the probability of coal and gas outbursts was calculated using Fig. 7 A schematic PRA framework showing different parametric calculations and triggering criteria to forecast rockbursts, gas outbursts and coal and gas outbursts in retreating LTCC mining, where l block is the length of the coal block, w block is the width of the coal block, h block is the height of the coal block, is the coal porosity, c is the cohesion, is the angle of internal friction, ρ is the density, s is the new surface area, w is the specific energy of coal, p g is the initial gas pressure in the coal matrix, P atm is the atmospheric pressure, γ is the adiabatic constant, V ej is the ejection velocity of the coal/rock mass, n is the vertical stress and p n is the pore pressure in the coal seam at n th excavation step, V gas is the volume of gas in coal, ERR i is the incremental energy release rate, V m is the volume of coal, c is the uniaxial compressive strength, W c is the crushing energy required per unit volume of coal, V free is the volume of free gas present in the cracks and fractures, m f is the mass of the coal/rock block being ejected, V a is the volume of gas adsorbed from the coal, W free is the free gas energy, W adsorbed is the adsorbed gas energy, t is the tensile strength, and W k is the transportation energy required for violent ejection Fig. 8 Flowchart representing the PRA framework developed Fig. 9 The PRA framework implemented in GoldSim both strain energy and gas expansion energy contributions (Eq.12).It was considered risk-free when none of the triggering criteria was true.A flowchart representing the PRA framework developed is presented in Fig. 8.
The PRA framework developed (Figs. 7, 8) was implemented in GoldSim, a dynamic and probabilistic risk assessment software, by linking influences and non-linear dependencies to represent the complex engineering system (Fig. 9).Depending on the nature of input values (fixed or stochastic), suitable data elements were selected, and distributions were assigned.Expression elements were used to graphically link relevant data elements, represent their interdependencies and calculate values.As per Eqs.(1-2), the stochastic energy release rate ( ERR ) obtained from the numerical model was used to calculate the energy release rate ( ERR i ).The length ( l block ), width ( w block ) and height ( h block ) of coal block were multiplied to represent coal volume ( V m ).The coal volume ( V m ), porosity ( ), gas pressure in the coal matrix ( p g ), atmospheric pressure ( P atm ) and adiabatic constant (γ) were used to calculate the free gas present in the coal ( V f ree ).Fol- lowing Eq. ( 4), the stochastic volume of gas released ( V gas ) calculated from the numerical model and the volume of free gas ( V f ree ) were used to calculate the volume of adsorbed gas present in the coal seam ( V a ).
As per Eq. ( 5), the coal volume ( V m ), the volume of free gas ( V f ree ), adiabatic constant (γ), atmospheric pressure ( P atm ), and final gas pressure ( p n ) were used to calculate the free gas energy ( W f ree ).The volume of adsorbed gas ( V a ), adiabatic constant (γ), atmospheric pressure ( P atm ), final gas pressure ( p n ) and coal volume ( V m ) were used to calculate the adsorbed gas energy ( W adsorbed ) from Eq. ( 6).The density of coal (ρ), specific energy of coal ( w ), and new surface area ( s ) were used to calculate the crushing energy of coal ( W c ).The coal volume ( V m ) and density (ρ) were used to calculate the coal mass ( m f ).The coal mass ( m f ), stochastic ejection velocity ( V ej ) and coal volume ( V m ) were used to calculate the required transportation energy ( W k ).The stochastic value of vertical stress ( n ) acting on coal seam at each excavation step was taken from the numerical model.The cohesion ( c ) and the angle of internal fric- tion ( ) were used to calculate the unconfined compressive strength ( c = 2ccos ∕(1 − sin ) ).The tensile strength ( t ) value was taken from the input to the numerical model.The correlations between different elements are represented using black arrows.The correlations and dependencies were verified in the model, as incorrect dependencies may result in inappropriate system responses, which can be difficult to debug in complex systems, increase computational times and give erroneous results (Mattenberger et al. 2015).
Status elements were used to define several triggering conditions.The conditional output from the status element (true/false) was represented using green arrows.The unconfined compressive strength ( c ), stress ( n ) and tensile strength ( t ) were used to calculate the strength conditions (  n >  c , c / t <14).As per Eq. ( 10), the trans- portation energy ( W k ), crushing energy ( W c ), incremental energy release rate ( ERR i ), and strength conditions were used to calculate the rockburst potential.The atmospheric pressure ( P atm ), final gas pressure ( p n ), free gas energy ( W f ree ), crushing energy ( W c ) and tensile strength of coal ( t ) were used to calculate the gas outburst conditions ( W f ree >W c , p n − P atm >  t ).The adsorbed gas energy ( W adsorbed ), crushing energy ( W c ) and gas outburst conditions were used to calculate the gas outburst potential from Eq. ( 11).Following Eq. ( 12), the crushing energy ( W c ), transportation energy ( W k ), incremental energy release rate ( ERR i ), and free gas energy ( W f ree ) were used to calculate the coal and gas outburst potential when rockburst potential is false and strength conditions are true.The crushing energy ( W c ) and adsorbed gas energy ( W adsorbed ) were used to calculate the quasi-static coal failure and excessive gas emissions potential ( W adsorbed > W c ) when gas outburst potential, rockburst potential and coal and gas outburst potential were false, while strength conditions were true.The safe mining conditions were concluded when the rockburst potential, gas outburst potential, coal and gas outburst potential, and strength conditions were false.
Monte Carlo simulations using Latin Hypercube sampling were used to generate independent and equally likely realisations of the model, uniformly spanning the values of stochastic parameters to represent uncertainty.The model was then simulated through time at appropriate timesteps to monitor the future state of the model.Result elements were used to monitor the probability distribution.The values of fixed parameters fed into the numerical model and the stochastic parameters obtained from the results of numerical modelling were fed into the framework using the data elements and all calculations were performed in GoldSim to run different scenarios and ultimately give the probability of rockbursts, gas outbursts and coal and gas outbursts.
Based on field monitoring data, several researchers have confirmed that most microseismic activities were observed within 100 m ahead of the longwall face (e.g., Cao et al. 2020).Therefore, the vertical stress acting on the coal seam, the displacement of the roof due to mining, pore pressure, and the incremental energy release rate were calculated within 100 m ahead of the active longwall face.Histograms of the vertical stress and strain energy were plotted, and a suitable distribution was fitted to calculate the stochastic variation in these parameters as the LTCC face retreats.The volume of total gas emission was calculated for each excavation step and fed into the model as stochastic parameters with the Poisson distribution.The values of these parameters vary at each excavation step and are also sensitive to the change in values of geomechanical and reservoir properties.These uncertainties were fed into the GoldSim model as carbon dioxide emission rate, and f total gas emission rate, at different excavation steps stochastic parameters.To realistically incorporate the uncertainty in the stochastic parameters throughout the model and the output, the model was simulated through time at a timestep of 1 day for 100 days using Monte Carlo simulation and the entire model was run for 100,000 independent and equally likely realisations in GoldSim.

Parametric Analysis of the Effects of Lithological Heterogeneity
The coupled geomechanical and gas flow modelling results were used to provide system feedback by generating statistically sufficient data to represent dynamically changing mining conditions at the LTCC.Fixed parameters taken from the coupled models are l block =5 m, w block =5 m, h block =3 m, p g =0.7 MPa, P atm =0.085 MPa, and =2360 kg/m 3 .The crushing energy of coal was calculated using values s = 0.015 m 2 /g and w = 10 J/m 2 taken from the literature (Cheng  et al. 2021).The velocity of ejection ( V ej ) was implemented as a stochastic having a Poisson distribution with a mean value of 10 m/s.

Effect of Xylite Distribution
The xylite distribution may vary throughout the deposit.To accommodate this variability, three different distributions were considered to represent 30%, 60% and 90% xylite content in the heterogeneous zone (Fig. 3b-d).These distributions exhibit different face heterogeneity levels at the longwall face as it retreats through the heterogeneous zone.The layer properties were kept the same as listed in Tables 1  and 2. The peak vertical stress acting ahead of the longwall face was compared at the 25th excavation step where the variability was maximum.The abutment stress increased from ~ 15.62 MPa for 30% and 60% xylite to ~ 18.21 MPa for the 90% xylite scenario (Fig. 10a).The pore pressure accumulated ahead of the longwall face at the 25th excavation step increased from ~ 0.88 MPa for 30% and 90% xylite to ~ 0.92 MPa for 60% xylite (Fig. 10b).However, as the pore pressure variation was only ~ 0.04 MPa, the change was assumed to be insignificant.A large volume of coal failed initially due to the vertical stress redistribution at the start of LTCC mining.After the fifth excavation step, the longwall face entered the heterogeneous zone with the volume of the failed zone depending on the geomechanical properties of different lithotypes.The variation in the volume of the failed zone is less for 30% xylite distribution because a relatively small volume of strong xylite is present in the model.The variation increases with the increase in xylite percentage which can be attributed to the increased volume of strong xylite blocks as the longwall retreat through the heterogeneous zone (Fig. 10c).The emission rate of methane (Fig. 10d) and carbon dioxide (Fig. 10e) across the longwall face varies inversely with the face heterogeneity and with the xylite percentage.The minimum gas emission rate occurs at the 18th excavation step for 30% xylite distribution and at 15th excavation step for 60% and 90% xylite distributions which corresponds to the high face heterogeneity at the active longwall face.A large volume of gas is released as the longwall retreats past the hard conditions to soft conditions as the face heterogeneity decreases.It can also be seen that the volume of gas released decreases sharply with the increase in xylite distributions in the heterogeneous zone (Fig. 10d-f) which is natural as xylite has lower porosity and significantly lower gas content (Fig. 6a) and permeability.Vertical stress acting ahead of the longwall face (within 100 m) roughly followed a normal distribution with a mean value of 9.45 ± 0.5 MPa (Fig. 11a).The pore pressure acting ahead of the longwall face does not show any appreciable change due to mining and mostly has a value of ~ 0.86 MPa, thus a fixed value was selected, and it was not stochastically varied in the PRA model (Fig. 11b).The variation in ERR i is large and roughly follows a log-normal distribution.Figure 11c shows a representative plot of the ERR i distribution for the change in xylite percentage.The face heterogeneity at the active longwall face keeps changing at every excavation step influencing the volume of gas emissions (Fig. 11d).The mean and standard deviation of ERR i and the volume of total gas emission at each excavation step for xylite percentage is listed in Table 3.
Figure 12a shows the probability of rockburst occurrence at different excavation steps for the change in xylite percentages.The variation in the probability is pronounced in the heterogeneous zone.For 30% xylite distribution, the heterogeneous zone has fewer xylite blocks and thus the influence of heterogeneity is less pronounced with rockburst probability peaking at the 25th excavation step only which corresponds to the situation when mining has retreated past the hard conditions and the face heterogeneity decreased to a very low value.For 60% xylite distribution, the number of xylite blocks in the heterogeneous zone is more and affects the stress and gas accumulation.As the mining retreats through the heterogeneous zone, the probability of rockbursts peak can be seen at the 10th excavation step which corresponds to the sudden increase in face heterogeneity (Fig. 11d) and at the 25th excavation step, when the mining retreats past the heterogeneous zone.For 90% xylite distribution, a single peak is seen around the 15th excavation step that corresponds to maximum face heterogeneity (Fig. 11d).Similarly, the probability of coal and gas outburst occurrence due to the change in the xylite distribution follows a similar trend as that for rockbursts (Fig. 12b).However, the probability of occurrence is relatively less than that for rockbursts.This can be attributed to the fact that only realisations where ERR i and free gas energy are more than the crushing energy was considered.For 30% xylite distribution, the probability is very low with a single peak at the 25th excavation step, as mining retreats past the hard conditions and face heterogeneity drops.For 60% and 90% xylite distribution, a similar trend was observed as for rockbursts, which is dependent on the face heterogeneity.The low probability can be related to the fact that the pore pressure (~ 0.86 MPa) developed in the coal seam was not sufficient to overcome the tensile strength of coal (0.92 MPa) and thus tensile failure was not observed in the model.
The probability of gas outbursts is nil in all excavation steps, which can be attributed to the fact that the gas content in the coal seam is not high, thus the crushing energy required for gas outbursts is not met.The probability of occurrence of rockbursts and coal and gas outbursts hazard is very low (in the order of 10 -3 ).It can be judged that it is mostly safe to continue retreat mining, but the operators should be more vigilant when the face heterogeneity varies significantly.
To verify the accuracy of the GoldSim model, 100,000 random scenarios considered were also evaluated in terms of the effect of % xylite distribution in the coal seam on the probability of rockburst and gas outburst occurrence using MS-Excel.The PRA framework implemented in the Gold-Sim was replicated in Excel by specifying the input parameters assigned to and stochastic variations obtained from the numerical models as input values for calculations.Verification analysis has shown general consistency between the two models, with only slight deviations attributed to rounding off errors of significant digits.Therefore, the GoldSim model can be used to determine the probability of rockbursts, gas outbursts and coal and gas outbursts for other scenarios.The analysis of different xylite distributions shows that 60% and 90% xylite had a very similar response to face advance and change in reservoir conditions due to coal production, which is quite different from the 30% xylite scenario.As 90% xylite distribution represents an extreme scenario that is unlikely to occur in most coal seams, 30% and 60% xylite scenarios were considered for further parametric analysis.Since no appreciable change in pore pressure was observed, the effect of changes in lithology on pore pressure was not analysed further.

The role of Geomechanical Properties of Xylite
Further analyses were conducted by increasing the stiffness of xylite.The cohesion, angle of internal friction, tensile strength, and residual cohesion of xylite blocks were increased while the reservoir properties of xylite were kept unchanged (Table 4).The properties of other layers in the model were also kept fixed as in Tables 1 and 2 (referred to as the base case).This resulted in an increase in the UCS of xylite from 11.26 MPa for the base case to 16.32 MPa for Case 1 and 24.13 MPa for Case 2.
The peak vertical stress acting ahead of the longwall face varied in the range ~ 15.64-16.89MPa for 30% xylite (Fig. 13a) and in the range ~ 15.86-17.47MPa for 60% xylite (Fig. 13b) for the increase in stiffness of xylite.For 30% xylite, the peak vertical stress occurs further inside the longwall on solid coal, showing fractured coal ahead of the face.For 60% xylite, the peak vertical stress occurs ahead of the face on solid coal for the base case (UCS = 11.26MPa) and Case 1 (UCS = 16.32 MPa), and some fractures are induced at the face.As the stiffness of xylite is increased to a very high value in Case 2 (UCS = 24.13MPa), the vertical stress peaks further towards the faceline.This can be attributed to the fact that the UCS of xylite is higher than the maximum stress acting on the face (~ 16.17 MPa), thus no appreciable fractures are induced in coal, and it can withstand the stress abutment acting on the face.
As Fig. 13c, d illustrates, the volume of failed zones is inversely proportional to the stiffness of xylite.As the stiffness of xylite increases from the base case, the volume of failed zone decreases.The volume of the failed zone is also sensitive to face heterogeneity as longwall retreats through the heterogeneous zone.It decreases as the face heterogeneity increases.For 30% xylite, as the longwall retreats past the hard condition (high face heterogeneity), the volume of failed zone increases to a higher value than that at the start of extraction, suggesting more fractures are induced due to mining-induced stress abutment (Fig. 13c).For 60% xylite, the effect of stiffer xylite is more prominent as the volume of failed zone decreases significantly in the case of very stiff xylite (Case 2).Even after crossing the hard conditions, the volume of failed zone remains low as compared to the value at the start of extraction for Case 2 (Fig. 13d).These observations suggest that a relatively high distribution of stronger xylite can control the volume of the failed zone, thereby preventing high permeability pathways from being formed, blocking gas movement towards the face, and reducing emissions at the active face.The results exemplify that the volume of the failed zone depends strongly on the spatial distribution, geomechanical properties and the relative abundance of xylite in the heterogeneous zone.The rate of total gas emission obeys the face heterogeneity while retreating through the heterogeneous zone (Fig. 13e, f).In general, with the increase in stiffness of xylite, the rate of gas emission decreases sharply, which can be attributed to the lesser volume of failed zones.As the longwall passes through the hard conditions, a sudden spike in the rate of total gas emission occurs.This spike is steeper and higher for Case 2 as compared to that for Case 1 or the base case for 30% xylite (Fig. 13e).For 60% xylite, the rate of total gas emission decreases sharply from around 120 m 3 /min to around 10 m 3 /min as the stiffness of xylite is increased, demonstrating the barrier effect of a large spatial distribution of stiff xylite to gas flow.This is followed by a sudden increase in gas emission rate once the face heterogeneity reduces (Fig. 13f), which matches the observation of previous researchers who have suggested that strong, low permeability xylite acts as a barrier to gas flow.
Vertical stress acting within 100 m ahead of the longwall face roughly followed a normal distribution with a mean value of 9.45 ± 0.5 MPa for both 30% and 60% xylite distri- butions (Fig. 14a, b).A small peak is found around 15 MPa that corresponds to the maximum stress abutment acting at the active face, which increases in magnitude when the stiffness of xylite is increased.
The effect is more pronounced for 60% xylite distribution (Case 2).The pore pressure acting ahead of the longwall face does not show any appreciable change due to mining and mostly has a value of ~ 0.86 MPa (Fig. 14c, d).The variation in ERR i follows a log-normal distribution.A representative plot of the ERR i distribution for 30% xylite distribution is shown in Fig. 14e and for 60% xylite distribution is shown in Fig. 14f.The mean and standard deviation of ERR i and the volume of total gas emission at each excavation step for the change in geomechanical properties of xylite are listed in Table 5.
Figure 15 shows the probability distribution for rockbursts and coal and gas outbursts with the change in the geomechanical properties of xylite.For 30% xylite distribution (Fig. 15a, b), the probability of the occurrence of rockbursts and coal and gas outbursts are very low at all excavation steps except for the 25th step which corresponds to the soft mining conditions when the longwall face has retreated past the zone of high face heterogeneity.It can also be seen that for very stiff xylite (Case 2), the probability is higher for both rockbursts and coal and gas outbursts as compared to that for Case 1 and the base case.Similarly, for 60% xylite distributions (Fig. 15c, d), the probability of rockbursts and coal and gas outbursts peaks at the 25th excavation step with stiff xylite (Case 2) having the maximum probability.The probability of rockbursts and coal and gas outbursts increases as the stiffness of xylite is increased.The small fluctuations observed in the heterogeneous zone confirm the dependence of rockbursts and coal and gas outbursts on face heterogeneity.The probability is relatively higher for the 60% xylite scenario, however, the value is in the order of 10 -3 , suggesting it is mostly safe to continue mining in the representative scenarios.

The Role of Gas Storage Capacity of Detrite
The role of the gas storage capacity of detrite on the probability of rockbursts and coal and gas outbursts were also analysed by assigning different Langmuir parameters for detrite.The geomechanical properties of xylite were maintained as they were for Case 2 (Table 4), referred to as the base case this time.Two new Langmuir isotherms were constructed with higher CO 2 sorption/desorption rates for detrite (Table 6 and Fig. 6b).It can be seen from Fig. 16a that, for the scenarios with increased gas storage capacity (therefore desorption rate) for detrite, the rate of total gas emission increases further as compared to the cases in Fig. 13e and f for both 30% and 60% xylite distributions.The impact is more prominent for 60% xylite where the rate of total gas emission decreases sharply to around 10 m 3 /min following the increase in face heterogeneity and increases rapidly to around 150 m 3 /min once the hard conditions are crossed by the retreating longwall face (Fig. 16b).Increased gas holding capacity for detrite in a high lithological heterogeneity scenario further suggests that, if potential coal and gas outburst conditions are created by much stronger xylite layers dominating the seam, the probability for coal and gas outbursts will be increased further.
The geomechanical properties were not changed, thus, the stress distribution obtained for this case remains the same at 9.45 ± 0.5 MPa for both 30% and 60% xylite distributions (Fig. 14a, b).The pore pressure acting ahead of the longwall face showed some change due to the change in Langmuir parameters.A slight increase in pore pressure can be seen but it mostly has a value around ~ 0.86 MPa (Fig. 17a, b).A representative plot of the ERR i showing log-normal distribu- tion for 30% xylite is shown in Fig. 17c and for 60% xylite distribution is shown in Fig. 17d.The mean and standard deviation of ERR i and the volume of total gas emission at each excavation step for the change in Langmuir properties are listed in Table 7.
Figure 18 shows the probability distribution for rockbursts and coal and gas outbursts with the change in Langmuir properties at different xylite percentages.The probability of rockbursts and coal and gas outbursts is very low for  18a, b).The change in Langmuir properties does not affect the probability of rockbursts and gas outbursts for a very stiff xylite.For the 60% xylite distribution (Fig. 18c, d), the probability of rockbursts and coal and gas outbursts shows some variation.It peaks just after entering the heterogeneous zone as the face heterogeneity starts to increase.For the hard conditions (when face heterogeneity is high), the probability remains low mostly because the change in Langmuir properties does not have a significant effect and the gas emission rate is low.

Conclusions
This paper presented a new probabilistic risk assessment methodology to realistically assess the complex rockburst and gas outburst occurrences using stress conditions, energy conditions and the bursting tendency of coal.The limitations of earlier developed approaches were overcome by accommodating feedback from a coupled geomechanics and gas flow model and implementing uncertainty in the mining parameters.The PRA framework developed is easy to adapt to other underground mining scenarios and offers the flexibility to modify the triggering criteria to suit the mining conditions.The results obtained from the PRA are representative of the actual field observations which give the   2.5 15 25 Goaf 0.63 30 0.52 0.63 0.52 0.01 0.30 3 × 10 -14 3 × 10 -14 ----confidence for it to be used as a versatile tool to forecast the probability of hazards in underground mining by suitably incorporating the mining conditions.The research described in this paper suggests that, given the heterogenous coal seam characteristics present at Coal Mine Velenje, the probability of rockbursts is generally higher as compared to the probability of coal and gas outbursts.This can mostly be attributed to the low pore pressure developed in the coal seam, which is less than the tensile strength of coal, thus, tensile failure is restricted in most cases.However, high-stress abutment acting ahead of the longwall face does not rule out the chance of excessive gas emission at the longwall face.The triggering criteria considered in the PRA framework suggest that the longwall mining is safe to conduct, however, this needs to be tested in the field by correlating it with other field monitoring data.The PRA framework can be further improved to provide shortterm/real-time analysis of the probability of rockbursts, gas outbursts and coal and gas outbursts by incorporating the field monitoring data into the model.
It was found that the total volume of gas emission is highly sensitive to face heterogeneity.The emitted gas includes free gas as well as adsorbed gas.The volume of free gas calculated was significantly low as compared to the volume of adsorbed gas per unit volume of coal extraction.The free gas energy was calculated to be between 2 and 4% of the adsorbed gas energy.This matches the observations of Cheng et al. (2021), who also found adsorbed gas energy to be the main reason behind gas outbursts and concluded that gas desorption energy is the major contributor to gas expansion energy.As the free gas energy was found to be less than the crushing energy (0.354 MJ/m 3 ), the probability of gas outburst occurrence is unlikely.In the case of coal and gas outbursts, ERR i may combine with the free gas energy to cumulatively exceed the crushing energy and thus a low probability for coal and gas outbursts was observed, which suggests the framework is reliable.
It was also evident from the analysis that since a large volume of adsorbed gas is present and in most situations the strength criteria were met, excessive gas emissions are highly likely.This observation made from the PRA framework supports the field conditions reported by Si et al. (2015a), who found that excessive gas emissions are a regular phenomenon at Coal Mine Velenje.

Fig. 2 aFig. 3 a
Fig. 2 a Velenje coal deposit, b schematic of multi-level and LTCC mining (after Si et al. 2015b) and c LTCC mining method at Coal Mine Velenje (after Jeromel et al. 2010)

Fig. 4
Fig. 4 Lithology and gas pressure gradient implemented in the Coal Mine Velenje model

Fig. 5
Fig. 5 A mining model developed with a full penetration of the heterogeneous zone in FLAC 3D b 90% distribution of xylite in the heterogeneous zone, c a cross-section across XX along which the vertical

Fig. 6 a
Fig. 6 a Pure CO 2 and CH 4 isotherms for Velenje lignite used in the heterogeneous coupled model, b additional Langmuir isotherms considered for detrite

Fig. 10
Fig. 10 The influence of xylite percentages on a vertical stress distribution at 25th excavation step, b pore pressure distribution at 25th excavation step, c volume of failed zones, d methane emission rate, e

Fig. 11
Fig. 11 Plots for the change in xylite percentage a normal distribution of vertical stress acting within 100 m ahead of the longwall face, b pore pressure variation acting within 100 m ahead of the longwall

Fig. 12
Fig. 12 The probability of a rockbursts and b coal and gas outbursts at different xylite distributions

Fig. 13
Fig.13The influence of geomechanical properties of xylite on a vertical stress distribution at 25th excavation step for 30% xylite, a vertical stress distribution at 25th excavation step for 60% xylite, c the

Fig. 14
Fig. 14 Variation in the a vertical stress for 30% xylite distribution, b vertical stress for 60% xylite distribution, c pore pressure for 30% xylite distribution, d pore pressure for 60% xylite distribution, e energy release rate at 23rd excavation step for 30% xylite distribution, and f energy release rate at 22nd excavation step for 60% xylite distribution ◂

Fig. 15
Fig. 15 Probability distribution for a rockbursts, b coal and gas outbursts for 30% xylite distribution, c rockbursts and d coal and gas outbursts for 60% xylite distribution with the change in geomechanical properties of xylite

Fig. 16
Fig. 16 The influence of reservoir properties of detrite on total gas emission rates at different excavation steps for the change in Langmuir properties a 30% xylite, and b 60% xylite

Fig. 17
Fig. 17 Variation in the pore pressure for a 30% xylite distribution, b 60% xylite distribution and energy release rate c at 29th excavation step for 30% xylite distribution, and d at 19th excavation step for 60% xylite distribution for the change in Langmuir parameters

Fig. 18
Fig. 18 Probability distribution for a rockbursts, b coal and gas outbursts for 30% xylite distribution, c rockbursts and d coal and gas outbursts for 60% xylite distribution with the change in Langmuir parameters

Table 2
Si et al. 2015c)ir properties of different layers/ lithotypes considered for numerical simulation (afterSi et al. 2015c)

Table 3
Stochastic input for the change in the energy release rate (log scale) and gas emission due to variation in xylite distributions in the coal seam

Table 5
Stochastic input for the change in the energy release rate (log scale) and gas emission due to variation in geomechanical properties

Table 6
Coal heterogeneity scenarios with varied Langmuir properties for detrite while most other layer properties are kept fixed v P L(CH 4 ) P L(CO 2 ) V L(CH 4 ) V L(CO 2 )

Table 7
Stochastic input for the change in the energy release rate (log scale) and gas emission due to variation in Langmuir properties