Coupled Heat and Mass Transfer Modelling of Coal Self-Heating in Longwall Goaf Areas with Spatially Variable Permeability

1. Introduction

Coal self-heating in goaf zones is one of the most significant fire hazards in underground coal mining. The process results from the oxidation of residual coal in contact with oxygen transported by airflow migrating through the porous rock mass. This phenomenon may lead to spontaneous combustion, posing a serious threat to the safety of mining operations.

In the years 2020–2024, a total of 13 endogenous fires caused by spontaneous coal combustion occurred in underground coal mines in Poland, including five in active mining areas [1]. These incidents directly affected the safety of mining personnel; in 2024 alone, 1166 miners were evacuated from operational areas, including 18 using self-rescue equipment. These data clearly demonstrate the importance of developing effective methods for preventing fire hazards associated with coal self-heating.

The mechanism of coal self-heating has been widely investigated in the literature [2,3,4,5,6,7]. The studies indicate that the process is governed by coupled phenomena, including gas flow through porous media, mass transport, chemical reactions, and heat accumulation [8,9,10]. The oxidation of coal is a surface process whose intensity depends on temperature, oxygen availability, and the structure of the porous medium. The two gaseous components, CO and CO₂, produced as a result of the self-heating of coal, change the intensity of the gas flow in the open pores of the goaf layer due to the carbon they contain and thus their higher molar masses. The oxidation of hard coal is in fact a surface reaction. In many cases, computational fluid dynamics (CFD) methods have been applied to simulate the distribution of oxygen concentration, temperature, and gas composition in the goaf [11,12]. Recent studies have focused on the identification of spontaneous combustion zones, fire early-warning systems, and advanced numerical simulations of coupled gas flow and heat transfer processes in goaf environments [13,14]. Comprehensive reviews have also highlighted the growing importance of integrated modelling approaches for fire prevention and risk assessment in underground coal mines [15]. These models provide valuable insights into the development of self-heating zones and fire risk, nonetheless, have often been based on simplified assumptions, such as isothermal conditions (approximately 25°C) or steady-state airflow and constant methane supply. In many cases, the permeability of the goaf has been assumed to be constant, which limits the ability of such models to accurately represent real mining conditions, where the permeability of the caving zone depends strongly on the geological composition and fragmentation of the rock mass [16].

The aim of this study is to develop a mathematical model describing the initiation of fire in goaf zones resulting from coal self-heating under non-isothermal and transient conditions. The model accounts for airflow filtration, continuous methane emission, and the possibility of inertisation using nitrogen. In addition, it incorporates coupled chemical reactions and heat transfer processes occurring in both gas and solid phases. The present work extends the analysis to include non-isothermal conditions and chemical reactions. These processes involve heat generation and accumulation, as well as changes in gas composition over time in both the gas and solid phases. As temperature increases, the rate of oxidation of both carbon and methane intensifies, leading to an increased risk of fire development. The proposed mathematical model simulates the non-stationary process. The formulation requires the definition of initial conditions describing the spatial distribution of key physical fields, including pressure, temperature, molar concentrations, and gas velocity. The solution obtained in previous studies is used as the initial state for the simulations presented in this work [2,5,17].

It should be emphasised that the permeability distribution used in the present model was determined based on in-situ mining data obtained in previous studies [2,5]. Considerable progress has been made in the numerical modelling of gas flow and spontaneous coal combustion in longwall goaf areas, particularly using CFD-based approaches [9,12,15]. These studies have improved the understanding of coupled transport and reaction processes in porous mining environments. However, many existing models employ simplified permeability distributions and do not explicitly represent spatial variability derived from field data.

From a thermal engineering perspective, coal self-heating in a goaf represents a transient reactive heat and mass transfer problem in a heterogeneous porous medium. Heat generation due to chemical reactions, conductive and convective heat transport, species transport, gas filtration, and heat exchange between the gas and solid phases are strongly coupled. Consequently, predicting the development of local temperature increases requires the simultaneous solution of the governing thermal, flow, species transport, and reaction equations.

Although substantial progress has been achieved in modelling spontaneous combustion in goaf areas, there remains a need for numerical approaches that simultaneously account for transient gas flow, non-isothermal heat transfer, coupled heterogeneous and homogeneous chemical reactions, and spatially variable permeability derived from in-situ data. The present study addresses this need by developing a coupled mathematical model of heat and mass transfer in a reactive porous medium representing a longwall goaf.

The principal contribution of this work is the simultaneous prediction of gas- and solid-phase temperature fields, gas composition, pressure, and gas velocity under transient operating conditions. In contrast to approaches based on uniform permeability distributions, the present model incorporates spatially variable permeability reflecting the geological composition and fragmentation of the caving zone [2,5,16,19]. The model is applied to U-type and Y-type ventilation configurations to analyse how different flow conditions affect oxygen transport, heat generation, temperature development, and the spatial evolution of zones susceptible to coal self-heating.

Although the model was developed specifically for longwall goaf environments, the mathematical formulation represents a more general class of coupled heat and mass transfer problems in reactive porous media. Application to other coal-bearing porous systems would require appropriate modification of the geometry, permeability distribution, initial conditions, and boundary conditions.

2. Description of the Mathematical Model

The system under consideration is represented as a horizontally oriented rectangular prism corresponding to the goaf formed after coal seam extraction. The characteristic dimensions of the prism range from 200 to 600 m. Part of its boundaries is adjacent to ventilation galleries, along which airflow generated by the main ventilation system is transported, while the remaining boundaries are in contact with undisturbed rock mass. For modelling purposes, transport of mass and heat in directions perpendicular to the main airflow is neglected. The airflow within the galleries induces filtration of gases through the collapsed rock mass containing coal residues originating from mining operations. The present numerical analysis focuses on the possibility of fire initiation resulting from the oxidation of these coal admixtures. It should be noted that the actual average height of the goaf ($H\cong 7,0$) his relatively small compared to its horizontal dimensions; therefore, the problem is treated as two-dimensional. In contrast to three-dimensional analyses commonly reported in the literature [18,20,21], two-dimensional approaches are also well justified and have been successfully applied in similar studies [22]. As a consequence of this assumption, variations of the analysed intensive quantities along the vertical direction, i.e. along the goaf height H, are neglected.

2.1. Model Assumptions and Reaction Kinetics

The goaf is treated as a heterogeneous system consisting of rock fragments of various sizes containing carbon and a gas phase filling the pore space. At the initial stage of the process, the gas phase is composed of a mixture of ventilation air and methane released from the surrounding rock mass. As the gas mixture filters through the goaf, oxygen contained in the air reacts with carbon present in the solid phase. This is a heterogeneous, exothermic surface reaction that leads to local increases in temperature and changes in the composition of the gas mixture over time. The main product of this process is carbon mono and dioxide. At elevated temperatures (above approximately 400°C), additional reactions may occur in the gas phase, including the combustion of carbon monoxide and methane. These reactions are associated with a significant increase in heat generation and may lead to rapid temperature growth, locally exceeding 500°C, thereby increasing the risk of fire and explosion. The exchange of mass and heat takes place between two interacting phases: the solid phase (rock fragments containing carbon) and the gas phase. The oxidation of carbon occurs at the surface of solid particles, whereas the combustion of carbon monoxide and methane takes place in the gas phase and is treated as a homogeneous volumetric process.

For modelling purposes, the goaf is represented as a layer of loosely packed spherical particles with an average diameter d (input data for the calculation algorithm). Each particle is assumed to have a homogeneous composition consisting of carbon and inert mineral matter. Although the actual particle size distribution is variable, a representative diameter is used to simplify the calculations. Due to point contacts between particles, direct heat and mass exchange between solids is neglected. All interactions occur between the gas phase and the particle surface. The gas in the immediate vicinity of each particle is assumed to be locally homogeneous in terms of temperature and composition, which allows the formulation of heat transfer within a particle in spherical coordinates.

The open porosity of the goaf is estimated using the Blake–Kozeny approach [23], which relates permeability to porosity. Based on previously determined permeability values obtained from in-situ mining data [2,3,5,19], local porosity values are calculated by numerical inversion of this relationship. Using the local porosity, the surface development coefficient is determined, leading to the expression:

$$\xi ={\displaystyle \frac{6(1-\epsilon )}{d}},m^2·m^{-3}$$

(1)

The gas phase consists of the following components: nitrogen (N₂), oxygen (O₂), carbon monoxide (CO), carbon dioxide (CO₂), water vapour (H₂O), and methane (CH₄). The composition of the mixture is described using mass fractions w_i, where $i=1,2,\mathrm{...},N$ and the index assigned to a given component, and appropriate molar masses $M_i$, accordingly:

• nitrogen from the air N₂ - i=1, M₁=0.028 $kg·{mol}^{-1}$
• oxygen from the air O₂ - i=2, M₂=0.032 $kg·{mol}^{-1}$
• carbon monoxide (a product of carbon oxidation in the air) CO - i=3, M₃=0.028 $kg·{mol}^{-1}$
• carbon dioxide (a product of the combustion of carbon monoxide and methane) CO₂ - i=4, M₄=0.044 $kg·{mol}^{-1}$
• water vapour (a product of methane combustion) H₂O - i=5, M₅=0.018 $kg·{mol}^{-1}$
• methane CH₄ - i=N= 6, M₆=0.016 $kg·{mol}^{-1}$.

In addition, the calculations use the value of the atomic mass of carbon $M_s = 0.012$ kg/mol. The mass fractions satisfy the obvious relationship

$$\sum _{i=1}^N{\omega }_i=1$$

(2)

The molar concentrations of individual components are calculated based on the local gas density:

$${\gamma }_i={\displaystyle \frac{\rho {\omega }_i}{M_i}},i=1,2,\dots N, mol·m^{-3}$$

(3)

The density of the gas mixture is determined using the ideal gas law presented in the form of [2]:

$$\rho ={\displaystyle \frac{p}{RT{\sum }_{i=1}^N\frac{{\omega }_i}{M_i}}}\hspace{1em}$$

(4)

The total local pressure is composed of atmospheric pressure p₀ (a value that is approximately constant in the mining area) and pressure differential p* generated by the ventilation system ($p=p^{*}+p_0)$. The pressure gradient constitutes the driving force for gas flow through ventilated galleries and the goaf. The value of this parameter reaches $p^{*}=0$ at the end of the galleries, i.e. at the point where the gas mixture stream leaves the analysed area. Its highest value occurs at the inlet cross-section of the gallery. The molar fractions of individual components are defined as

$${\omega }_i={\displaystyle \frac{{\gamma }_i}{{\sum }_{i=1}^N{\gamma }_i}}$$

(5)

For the gas mixture considered, the dynamic viscosity (m) is determined as one of the key parameters required for simulation calculations. Its value is calculated using Wilke’s method [23], which accounts for the contribution of individual components of the gas mixture.

The model considers the course of three main irreversible chemical reactions related to the oxidation of hard coal occurring in the goaf area, as well as those occurring exclusively in the gas phase, i.e. related to the combustion of CO and CH₄:

(A)

carbon oxidation (heterogeneous surface reaction): $2\mathrm{C}+{\mathrm{O}}_2\to \mathrm{C}\mathrm{O}$

(B)

carbon monoxide combustion (gas-phase reaction): $\mathrm{C}\mathrm{O}+{\displaystyle \frac{1}{2}}{\mathrm{O}}_2\to \mathrm{C}{\mathrm{O}}_2$

(C)

methane combustion (gas-phase reaction); $\mathrm{C}{\mathrm{H}}_4+2{\mathrm{O}}_2\to \mathrm{C}{\mathrm{O}}_2+2{\mathrm{H}}_2\mathrm{O}$

Due to oxygen deficiency in the goaf, it is assumed that carbon oxidation leads primarily to the formation of carbon monoxide.

According to Cygankiewicz (2018), the oxidation process of carbon proceeds through three distinct stages, depending on temperature [24]:

preliminary self-heating $T_0<T<T_{kr}$,

secondary (accelerated) self-heating, $T_{kr}<T<T_{sp}$,

proper combustion of carbon. $T>T_{sp}$.

The oxidation process begins at the initial temperature T₀=298.0 K, with characteristic transition temperatures T_kr=360.0 K and T_sp=600 K. According to Cygankiewicz (2018) [24], the kinetics of this process can be described using a generalized Arrhenius-type expression:

$${\dot{\mathit{R}}}_{\mathit{A}}=-{\mathit{k}}_{\mathit{A}}{{\mathit{\gamma }}_2}^{\mathit{\alpha }},\hspace{1em}\mathit{m}\mathit{o}\mathit{l}{\mathit{O}}_2·{\mathit{m}}^{-3}{·\mathit{s}}^{-1}$$

(6)

where the exponent represents the reaction order. The reaction rate constant $k_A$ in the study [24] was defined replacing the quotient $E_A/R$ (where$E_A$ is the so-called reaction activation energy (A)), empirically denoted by the quantity T_s, K;

$$k_A=k_A^0\mathit{exp}({\displaystyle \frac{T_s}{T}})$$

(7)

Cygankiewicz (2018) provides numerical values of the kinetic parameters for all three temperature ranges corresponding to the different stages of the oxidation process. These parameters exhibit discontinuous changes when the local temperature exceeds the threshold values separating these ranges.

The above experimental data were adopted as input parameters in the present simulation model. The reaction rate expression (${\dot{R}}_A$) was further modified to represent a pseudo-homogeneous reaction (A) occurring at the outer surfaces of the spherical particles forming the goaf. Additionally, the carbon content in the solid phase was incorporated into the reaction rate formulation, as the presence of carbon is a necessary condition for reaction (A) to occur.

$${\dot{R}}_A^{*}={\displaystyle \frac{\xi }{{\xi }_s}}{f}^{\alpha }{\dot{R}}_{A}mol{O}_2·m^{-3}{·s}^{-1}$$

(8)

where ${\xi }_s$ denotes the surface development coefficient in the layer of crushed coal, determined based on experimental data reported by Cygankiewicz (2018) [24]. Its value is several orders of magnitude higher than the values obtained from calculations based on Eq. (1). The quantity (f) represents the mass fraction of carbon in the rock material within the goaf and (a) is an empirical exponent corresponding to the order of reaction (A) with respect to carbon content. In the present study, its value is assumed to be equal to 1. For simplicity, (f) is defined as the ratio of the mass of carbon contained in a single particle to its total mass at the initial stage of the oxidation process. This quantity is assumed to vary only with time, while remaining spatially uniform within an individual particle. Its decrease over time does not affect the particle size but leads to a reduction in its mass and an increase in closed porosity. This effect is taken into account in the subsequent part of the model. Taking into account the stoichiometry of reaction (A), the rate of carbon consumption within a single particle is described by the following differential equation:

$${\displaystyle \frac{d{\rho }_s}{dt}}=-{\displaystyle \frac{M_s}{d·{\xi }_s}}f{\dot{R}}_{A}kg·m^{-3}·s^{-1}$$

(9)

where (t) denotes time, (r_s) is the initial density of the solid phase, and the reaction rate is defined by Eq. (6).

Equation (9) forms part of a system of differential equations solved numerically in this study. Its solution provides the spatial and temporal distribution of the carbon mass fraction within the goaf (f=f(x,y)). The evolution of this quantity over time reflects the progressive consumption of carbon at different locations within the goaf during the fire process.

The subsequent reactions considered in the model, i.e. reactions (B) and (C) defined above, occur exclusively in the gas phase, provided that both the temperature and the molar fractions of the reactants are sufficiently high.

Their rates are expressed as:

−

for reaction (B) (carbon monoxide combustion rate);

$${\dot{R}}_B={\displaystyle \frac{\partial {\gamma }_3}{\partial t}}=-k_B{\gamma }_3({\gamma }_2{)}^{\mathrm{0,5}}molCO·m^{-3}·s^{-1}$$

(10)

−

for reaction (C) (the rate of methane combustion)

$${\dot{R}}_C={\displaystyle \frac{\partial {\gamma }_6}{\partial t}}=-k_C{\gamma }_6({\gamma }_2{)}^2, molC{\mathrm{H}}_4·m^{-3}·s^{-1}$$

(11)

In the detailed algorithm, the reaction rate expressions given by Eqs. (6), (10), and (11) are reformulated using Eq. (3), allowing them to be expressed as functions of the mass fractions of individual components and the local density of the gas mixture. The resulting expressions are not presented here due to their complexity.

Due to the volumetric nature of reactions (B) and (C), their rates are expressed in mol·m⁻³·s⁻¹. These rates are also described using Arrhenius-type relationships (Eq. (7)), with activation energy values taken from physicochemical data.

2.2. Mass Exchange During the Oxidation Process

In addition to filtration and chemical reactions (A), (B), and (C), mass transport in the goaf includes internal sources associated with methane inflow from the surrounding rock mass and carbon monoxide generation as a product of reaction (A). These components are introduced into the gas mixture in their pure form and act as mass sources, modifying both the composition of the gas mixture and the filtration flow. Methane emission from the adjacent rock mass is treated as a continuous process. The volumetric intensity of this source $\left({\dot{v}}_6\right)$ is determined based on technological data and is assumed to remain constant throughout the simulated period. The methane point source is defined as:

$${\dot{s}}_6^0={\displaystyle \frac{{\dot{v}}_6{\rho }_6^0}{H}}\hspace{0.33em},kgC{H}_4·m^{-3}·s^{-1}$$

(12)

where the density of pure methane under standard conditions (${\rho }_6^0$), i.e. at a temperature of $T_0$ and total pressure $p_0$is determined using Eq. (4) adapted for a single-component gas.

Similarly, the inflow of carbon monoxide is considered as a source term directly related to reaction (A). Based on its stoichiometry, the mass source of carbon monoxide is defined as:

$${\dot{s}}_3^0=-2{\dot{R}}_A^{*}{M}_3,\hspace{0.33em}kgCO·m^{-3}·s\hspace{0.33em}$$

(13)

The source terms corresponding to methane and carbon monoxide (pure components) are treated as internal mass sources supplying the gas mixture flowing through the goaf.

The influence of these sources on the mass fraction (i) of a given component in the gas mixture is described by:

$$\dot{s_i^0}=\dot{s_g^0}\left({\omega }_i^0-{\omega }_i\right), kg·m^{-3}{s}^{-1}$$

(14)

where the source term represents the mass supplied per unit volume and time, and the composition of the supplied gas is taken into account. In limiting cases, the supplied medium may consist of a pure component, for which $({\omega }_i^0=1$), or may not contain the considered component, for which (${\omega }_i^0=0.$).

The expression defined by Eq. (9) is used to formulate the mass balance for each component separately. In the case of carbon monoxide and methane, their concentrations increase due to internal source terms (an increase in the values of ${\omega }_3$ and ${\omega }_6$), whereas the mass fractions of the remaining components decrease accordingly to satisfy the overall mass balance. For volumetric reactions (B) and (C), no additional mass source term is introduced ($\dot{s_g^0}=0$), as these reactions satisfy the principle of mass conservation and occur entirely within the gas phase. As a result, they only alter the chemical composition of the gas mixture, leading to an increase in the concentration of products accompanied by a corresponding decrease in the concentration of reactants. The transport of gas mixture components during the fire process is described by a system of coupled partial differential equations supplemented with appropriate boundary conditions. These conditions are defined at the interfaces between the goaf and the surrounding rock mass, as well as along the vertical surfaces of the ventilation galleries through which the air flows. The system of governing equations is supplemented by energy balance equations in order to determine the temperature field within the goaf. All equations are strongly coupled, resulting in a nonlinear problem. The part of the mathematical model related to heat transfer, which enables the determination of temperature fields within the domain of dimensions (X, m) and (Y, m) is presented in the following section. The first set of coordinates corresponds to the direction perpendicular to the longwall face, whereas second is aligned with the direction of coal seam extraction.

To describe mass transport numerically, the following system of differential equations is solved:

a)

Filtration Equation

The gas flow through the goaf is described by a filtration equation obtained by combining Darcy’s law with the continuity equation, including total internal mass sources ($S_m)$. Assuming quasi-stationarity of the density field, the governing equation for pressure is given by:

$${\displaystyle \frac{\partial }{\partial x}}\left(K_x\rho {\displaystyle \frac{\partial p^{*}}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(K_y\rho {\displaystyle \frac{\partial p^{*}}{\partial y}}\right)+S_m=0$$

(15)

where $K_x={\kappa }_x/\mu$, $K_y={\kappa }_y/\mu$ and the coordinates of the Cartesian system are within the following ranges: $0\le x\le X$and $0\le y\le Y$. The formulation of Eq. (15) assumes quasi-stationarity of the density field, which is justified by the relatively short relaxation time of gas density compared to the time scale of the self-heating process. As a result, the accumulation term is neglected, and the equation takes an elliptic form [5].

The total source term ($S_m)$ is defined as the sum of methane inflow and carbon monoxide generation:

$$S_m={\dot{m}}_6+2{\dot{R}}_A^{*}{M}_{3}kg·m^{-3}{s}^{-1}$$

(16)

Due to oxygen consumption in reaction (A), the carbon monoxide source term represents a net mass contribution to the gas mixture.

Based on the pressure field obtained from Eq. (15), the linear velocity vector components of the gas mixture are calculated using Darcy’s law:

$$v_x=-K_x{\displaystyle \frac{\partial p^{*}}{\partial x}}{ v}_y=-K_y{\displaystyle \frac{\partial p^{*}}{\partial y}}m·s^{-1}$$

(17)

The resulting velocity field is required for the analysis of convective transport of gas components.

• b) Convection–diffusion Transport Equations

In order to determine the time-dependent mass fraction fields of individual components in the gas mixture, the components of the mass flux density vector are defined in the Cartesian coordinate system.

$${\dot{m}}_{xi}=v_x\rho {\omega }_i-\epsilon D_i{\displaystyle \frac{\partial \left(\rho {\omega }_i\right)}{\partial x}}\mathrm{and} {\dot{m}}_{yi}=v_y\rho {\omega }_i-\epsilon D_i{\displaystyle \frac{\partial \left(\rho {\omega }_i\right)}{\partial y}}$$

(18)

where the first term represents convective transport and the second term describes diffusive transport according to Fick’s law. The inclusion of porosity accounts for the heterogeneous structure of the goaf.

In order to formulate the local mass balance for individual components of the gas mixture in the goaf, the contribution of internal mass sources must be taken into account. As discussed earlier, these sources include methane inflow from the surrounding rock mass and carbon monoxide generation at the surface of solid particles as a product of reaction (A) (Eq. (12) and (13)). The total source terms ($S_i$) for N−1 components are defined separately for each component of the gas mixture:

✓ $i=1$, nitrogen (N₂) is treated as an inert component. Its variation results only from dilution due to the inflow of methane and carbon monoxide (a product of reaction (A)). Therefore, its total source term is determined solely based on Eq. (14) for ${\omega }_1^0=0.$

$$S_1=(s_3^0-s_6^0){\omega }_1$$

(19)

✓ $i=2$, oxygen (O₂) is a reactive component involved in all reactions (A), (B), and (C), leading to its consumption within the gas mixture. In addition, it is diluted in a manner analogous to nitrogen (${\omega }_2^0=0$). Its total source term is therefore given by:

$$S_2=M_2({\dot{R}}_A^{*}+\mathrm{0,5}{\dot{R}}_B+2{\dot{R}}_C)+(s_3^0-s_6^0){\omega }_2$$

(20)

✓ $i=3,$ carbon monoxide (CO) is generated as a product of reaction (A) at the surface of solid particles and is introduced into the gas phase in a pure form (${\omega }_3^0=1$). It is subsequently transported with the gas mixture and acts as a reactant in reaction (B). Taking into account Eq. (14), its total source term is defined as:

$$S_3=s_3^0(1-{\omega }_3)+M_3{\dot{R}}_B$$

(21)

✓ $i=4,$carbon dioxide (CO₂) and water vapour (H₂O) are products of volumetric gas-phase reactions (B) and (C). These components are also subject to dilution due to the inflow of methane and carbon monoxide so for them ${\omega }_4^0={\omega }_5^0=0$. Their total source terms are defined as:

$$S_4=-M_4({\dot{R}}_B+2{\dot{R}}_C)+(s_3^0-s_6^0){\omega }_4$$

(22)

$$S_5=-2M_5{\dot{R}}_C+(s_3^0-s_6^0){\omega }_5$$

(23)

These total source terms ensure consistency between chemical reactions and mass transport within the gas mixture.

The quantities ${\omega }_i$ appearing in Eqs. (14 - 18) represent the mass fractions of individual components in a fully mixed gas phase. These expressions ($S_i$) are further supplemented by reaction rate terms modified using Eq. (3), allowing the source terms to be expressed as functions of the mass fractions and the local gas density. The resulting expressions are not presented here due to their complexity; however, they are implemented in the numerical algorithm.

To determine the mass fraction fields (${\omega }_i$) for all components, a system of partial differential equations describing convective–diffusive transport is solved. These equations are derived from the mass balance of a given component within a differential control volume of the goaf, defined by horizontal dimensions dx$·$dy and a unit thickness (1 m for the differential time $dt$). The governing equation is derived using Eq. (4), which defines the mass flux densities of the i component across the boundaries of the control volume, combined with the corresponding source term ($S_i$). The transient character of the process requires inclusion of the accumulation term within the control volume. The general form of the governing equations is given as:

$${\displaystyle \frac{\partial }{\partial x}}\left(v_x\rho {\omega }_i\right)+{\displaystyle \frac{\partial }{\partial y}}\left(v_y\rho {\omega }_i\right)+{\displaystyle \frac{\partial \left(\rho {\omega }_i\right)}{\partial t}}=\epsilon D_i\left({\displaystyle \frac{{\partial }^2\left(\rho {\omega }_i\right)}{\partial x^2}}+{\displaystyle \frac{\begin{array}{c}\\ {\partial }^2\left(\rho {\omega }_i\right)\end{array}}{\partial y^2}}\right)+S_i$$

(24)

for i = 1, 2, …, N−1

The above equation is written in a simplified form based on the assumption of quasi-stationary density and Darcy-type flow in porous media. Under these conditions, the conservative form can be reduced to the presented formulation expressed in terms of mass fractions.

Methane (when i=N=6) is treated in the model as a complementary component of the gas mixture. Therefore, its mass fraction in both the goaf and the ventilation galleries is determined from Eq. (2). Accordingly, the following relationship is obtained:

$${\omega }_N=1-{\displaystyle {\sum }_{i=1}^{N-1}}{\omega }_i i=\mathrm{1,2},...,N$$

(25)

As mentioned earlier, the system of differential equations describing mass exchange and transport during a fire is supplemented by Eq. (9), which defines the rate of carbon consumption within spherical particles. However, this set of equations (N+1) is not sufficient to fully describe the fire initiation process. It is necessary to determine the temperature fields in both the solid phase (rock mass) and the filtering gas. During the self-heating process, temperature increases over time, which in turn accelerates the oxidation reactions.

2.3. Heat Exchange and Transport During A Fire

The increase in temperature within the goaf is caused by exothermic chemical reactions (A), (B), and (C). These processes are associated with heat release, the intensity of which depends on the rates of the individual reactions. Therefore, it is necessary to quantify these thermal effects based on appropriate thermochemical data. To this end, the local enthalpy change under isobaric conditions is defined for each reacting component. Nitrogen (N₂ - i=1) is excluded from this analysis, as it is chemically inert and does not participate in reactions (A), (B), or (C). The enthalpy changes of each gaseous component i = 2,3,…,N and carbon in the solid phase is expressed as:

$$\Delta H_{0,i}^T=\Delta H_i^{*}+{\int }_{T_0}^T{C}_{i}dT, J/mol$$

(26)

where molar heat ($C_i$) of $\mathit{i}$ component corresponds to $p=p_0=idem$, $T_0$ - initial temperature (T=298 K) and $T$ - is the local temperature at a given time during the self-heating process. These temperatures correspond to the initial and current thermodynamic states, respectively. The condition $\mathit{\Delta }{\mathit{H}}_{\mathit{i}}^{\mathit{*}}\ne 0$ applies only to chemical compounds (CO, CO₂, H₂O, CH₄), whereas for elements such as carbon and oxygen, $\mathit{\Delta }{\mathit{H}}_{\mathit{w}}^{\mathit{*}}=\mathit{\Delta }{\mathit{H}}_2^{\mathit{*}}=0$. The integral term in Eq. (26) accounts for the temperature (T) dependence of the molar heat capacity. This dependence is commonly approximated using the Kelley equation, which can be simplified to a linear function over a limited temperature range [25]. Under this assumption, the integral in Eq. (26) is approximated using the average value of the molar heat capacity ($\bar{\mathit{C}}(\mathit{T}-{\mathit{T}}_0)$) over the temperature interval. Based on the molar heat capacities of individual components, the local specific heat capacity of the gas mixture can be determined using the following relationship:

$${\mathit{c}}_{\mathit{g}}={\displaystyle {\sum }_{\mathit{i}=1}^{\mathit{N}}}{\displaystyle \frac{{\mathit{\omega }}_{\mathit{i}}{\mathit{C}}_{\mathit{i}}}{{\mathit{M}}_{\mathit{i}}}}-{\mathit{v}}_{\mathit{x}},\mathit{J}·{\mathit{k}\mathit{g}}^{-1}·{\mathit{K}}^{-1}$$

(27)

The average specific heat over a given temperature interval is denoted by an overbar. For an individual component, it is expressed as ${\bar{\mathit{c}}}_{\mathit{i}}={\displaystyle \frac{{\bar{\mathit{C}}}_{\mathit{i}}}{{\mathit{M}}_{\mathit{i}}}}$, while for the gas mixture it is denoted as ${\bar{\mathit{c}}}_{\mathit{g}}={\sum }_{\mathit{i}=1}^{\mathit{N}}{\mathit{\omega }}_{\mathit{i}}{\bar{\mathit{c}}}_{\mathit{i}}$. These quantities are used in the numerical algorithm.

The thermal effects associated with individual reactions (A), (B), and (C), conventionally expressed per mole of the reacting substrate, are determined using Kirchhoff’s law based on the quantities defined in Eq. (26). In the case considered, these effects are obtained from the following balance relations:

$$Q_A=2(H_{\mathrm{0,3}}^T-H_{0,w}^T)-H_{\mathrm{0,2}}^T,J·{mol}^{-1}{\mathrm{O}}_2$$

(28)

$$Q_B=H_{\mathrm{0,4}}^T-H_{\mathrm{0,3}}^T-\mathrm{0,5}{H}_{\mathrm{0,2}}^T, J·{mol}^{-1}\mathrm{C}\mathrm{O}$$

(29)

$$Q_C=H_{\mathrm{0,4}}^T+2H_{\mathrm{0,5}}^T-H_{\mathrm{0,6}}^T-2H_{\mathrm{0,2}}^{T}J·{mol}^{-1}{\mathrm{C}\mathrm{H}}_4$$

(30)

It should be noted that, according to Eq. (26), the thermochemical effects associated with possible phase transitions are generally negligible for most components involved in the oxidation process. In particular, no phase transitions occur for nitrogen and the majority of reacting species. An exception is methane. One of its combustion products, water vapour (H₂O), may undergo condensation at temperatures of approximately 373 K, releasing additional heat. However, this process occurs outside the fire zone, during the cooling of exhaust gases, and therefore does not affect the temperature field within the goaf. In this context, the thermal effect of combustion is interpreted as the calorific value. Based on the above considerations, the temperature field within the goaf is treated as being generated exclusively by internal heat sources associated with exothermic reactions (A), (B), and (C). Reaction (A) occurs at the surface of coal-containing particles; therefore, the corresponding heat source is also of a surface nature. Assuming local homogeneity of the gas phase in the vicinity of each particle, the temperature at the particle surface is also considered uniform. At the initial stage of self-heating, the temperatures of the gas and the solid surface are equal. However, as the process develops, temperature differences arise due to heat exchange between the gas and the solid phase. To summarise the assumptions regarding heat and mass transfer in the goaf, each differential control volume is assumed to be characterised by a homogeneous gas composition and a uniform temperature field. A similar assumption is adopted for the outer surface of the spherical particles. At the initial stage of the self-heating process, the temperatures of the gas and the solid surface are equal. As the process develops, temperature differences arise due to heat exchange between the gas phase and the solid particles. Due to the higher activation energy of reactions (B) and (C), the ignition of carbon monoxide and methane occurs only when the local gas temperature reaches sufficiently high values. At this stage, the gas temperature in certain regions of the goaf may exceed the temperature of the solid particles, resulting in a reversal of the direction of heat transfer between the two phases. In the initial phase of self-heating, when only reaction (A) takes place, heat is generated within the solid particles. This leads to an increase in their surface temperature and, consequently, to heating of the surrounding cooler gas. The heat exchange between the gas and the particle surface, occurring through convective heat transfer (so-called “gas penetration”), is described using Newton’s law of cooling and is introduced into the model as an internal heat source:

$$\dot{q_s}=\xi {\alpha }_{s-g}\left[T_s\right({\displaystyle \frac{d}{2}},x,y,t)-T_g(x,y,t\left)\right],W·m^{-3}$$

(31)

where T_s(r, x, y) denotes the temperature distribution within a spherical particle as a function of the radial coordinate x, y, and T_g(x, y) represents the temperature field of the gas phase in the goaf. The parameter a_s-g is the surface-averaged heat transfer coefficient in the gas–solid system. Its value is determined using the Froessling correlation [20], which relates the Nusselt number to the Reynolds and Prandtl numbers. Since the Reynolds number depends on the local gas velocity, the value of the heat transfer coefficient must be evaluated at each point of the flow domain. The detailed procedure for calculating local values of a_s-g is not presented here due to its complexity. In the gas phase, additional internal heat sources arise from the exothermic reactions (B) and (C). The intensity of these sources is proportional to the rates of carbon monoxide and methane combustion. Using Eqs. (10) and (11), the total heat source term can be expressed as:

$${\dot{q}}_{BC}={\dot{R}}_B{Q}_B+{\dot{R}}_C{Q}_C,W·m^{-3}$$

(32)

where the values of the thermal effects $Q_B$ and $Q_C$ are determined based on Eqs. (29 - 30).

The energy balance of the gas phase is additionally influenced by mass sources associated with methane inflow and carbon monoxide generation as a product of reaction (A). These processes introduce additional enthalpy fluxes into the gas mixture, which modify both its thermal state and mass flux.

The total mass flux consists of two components: methane (CH₄) originating from the surrounding rock mass and carbon monoxide (CO) generated by the heterogeneous reaction (A) at the surface of solid particles. Methane entering the goaf from relatively cold rock layers is assumed to have a constant initial temperature (T₀). In contrast, the temperature of carbon monoxide corresponds to the local surface temperature of the particles, which varies in both space and time due to the progression of reaction (A). The resulting internal heat source associated with the inflow of these components is expressed as:

$${\dot{q}}_{\mathrm{3,6}}=s_3^0{\bar{c}}_3\left[T_s({\displaystyle \frac{d}{2}},x,y,t)-T_g(x,y,t)\right]+s_6^0{\bar{c}}_6\left[T_0-T_g(x,y,t)\right],W·m^{-3}$$

(33)

where the molar heat capacities of carbon monoxide and methane are taken as average values over the temperature intervals specified in square brackets. It should be noted that the second term in Eq. (33) remains non-positive throughout the entire process, since the corresponding temperature difference is negative (T_g $\ge$T₀). As a result, the inflow of methane leads to a cooling effect on the gas mixture.

The total heat source within the gas phase in the goaf is obtained as the sum of the contributions defined by Eqs. (31 - 33). Therefore, the overall heat source term can be written as:

$$\dot{Q}={\dot{q}}_s+{\dot{q}}_{BC}+{\dot{q}}_{\mathrm{3,6}}$$

(34)

In order to derive the governing differential equation for the transient temperature field in the gas phase (T_g(x,y,t)), it is necessary to define the components of the heat flux density vector in the domain under consideration. Similarly to the mass transport model (Eq. (18)), the heat flux is expressed as the sum of convective and conductive contributions, in accordance with Fourier’s law. The corresponding expressions are given as follows:

$${\dot{q}}_x=v_x\rho c_g{T}_g-\epsilon {\lambda }_g{\displaystyle \frac{\partial T_g}{\partial x}}\mathrm{and} {\dot{q}}_y=v_y\rho c_g{T}_g-\epsilon {\lambda }_g{\displaystyle \frac{\partial T_g}{\partial y}},W·m^{-3}$$

(35)

where${\lambda }_g$ - thermal conductivity coefficient of the gas mixture.

The heat transport equation is derived analogously to the mass transport equation, based on the energy balance for a differential control volume (H$·$dx$·$dy). The resulting equation is given as follows:

$${\displaystyle \frac{\partial }{\partial x}}\left(v_x\rho c_g{T}_g\right)+{\displaystyle \frac{\partial }{\partial y}}\left(v_y\rho c_g{T}_g\right)+\rho c_g{\displaystyle \frac{\partial T_g}{\partial t}}=\epsilon {\lambda }_g\left({\displaystyle \frac{{\partial }^2{T}_g}{\partial x^2}}+{\displaystyle \frac{{\partial }^2{T}_g}{\partial y^2}}\right)+\dot{Q}$$

(36)

where (T_g(x,y,t) ) is the transient temperature field in the gas phase within the porous space of the goaf, (0 $\le$x $\le$X) and (0$\le$y$\le$ Y) are spatial coordinates, and (t $\ge$0) denotes the time elapsed since the onset of the self-heating process. The total heat source term ($\dot{Q}$) is defined according to Eq. (34). Within the spherical solid particles, heat transfer occurs exclusively by conduction. Based on the assumption of local homogeneity of the gas temperature in the vicinity of the particles, as well as uniform thermophysical properties within the solid phase, the temperature field inside each sphere can be considered radially symmetric. Therefore, it is appropriate to formulate the heat conduction equation in a spherical coordinate system. The governing equation is given as follows:

$${\displaystyle \frac{\chi }{r^2}}{\displaystyle \frac{\partial }{\partial r}}\left(r^2{\displaystyle \frac{\partial T_s}{\partial r}}\right)=(1-f_0+f){\displaystyle \frac{\partial T_s}{\partial t}}, \mathrm{for} 0<r\le {\displaystyle \frac{d}{2}},t\ge 0$$

(37)

where (r) is the radial coordinate, (T_s(r, x, y, t)) denotes the temperature distribution within a spherical particle at a given time and spatial location in the goaf, and (f₀=f(x, y, 0)) is the mass fraction of carbon in the solid phase at the initial stage of the process (i.e. at (t = 0 )). The factor (1 - f₀ + f) introduced in Eq. (37) accounts for the change in the effective heat capacity of the solid phase resulting from the progressive consumption of carbon during the self-heating process. As the carbon content decreases, the thermophysical properties of the material change, which is reflected in the modified transient term. The thermal diffusivity ($\chi$) is defined as: $\chi ={\displaystyle \frac{{\lambda }_s}{{\rho }_s{c}_s}}$ where ( ${\lambda }_s$ ) is the thermal conductivity of the solid, (${\rho }_s$) is its density, and ($c_s$) is the specific heat capacity. The term appearing in parentheses on the right-hand side of Eq. (37) accounts for the variation of the heat capacity of the solid phase during the self-heating process and at later stages of the fire. In contrast, changes in thermal conductivity due to carbon consumption are neglected. The existence of a minimum temperature at the centre of the sphere is expressed by the following boundary condition:

$${\displaystyle \frac{\partial T_s}{\partial r}}=0, \mathrm{for} r=0$$

(38)

Using Eqs. (1), (6) - (8), and (21), a boundary condition can be formulated that defines the surface heat source at the outer surface of a spherical particle. This condition accounts simultaneously for heat generation due to the surface reaction, heat conduction within the particle, and heat exchange with the surrounding gas phase. The corresponding boundary condition is expressed as follows:

$$-{\lambda }_s{\displaystyle \frac{\partial T_s}{\partial r}}+{\alpha }_{s-g}(T_g-T_s)+{\displaystyle \frac{{\dot{R}}_A}{\xi }}{Q}_A=0,\mathrm{for} r={\displaystyle \frac{d}{2}}$$

(39)

The last term on the left-hand side of Eq. (39) represents the rate of heat release per unit surface area of a spherical particle associated with the combustion of carbon within the solid phase (reaction (A)), i.e. the surface heat source. The remaining terms describe the distribution of this heat between the interior of the particle (heat conduction according to Fourier’s law) and the surrounding gas phase (heat transfer according to Newton’s law). It is assumed that the temperature field is homogeneous over the total surface area of particles within a unit volume assigned to a given point. A similar assumption is adopted for the gas phase. As a result, a temperature difference arises between the solid and gas phases during the process, which drives heat exchange between them. Locally, this heat transfer may be very intense. The system of governing equations, including Eqs. (9), (15), and (34), together with Eqs. (19) and (24) defining the methane mass fraction, describes a coupled set of physical fields within the goaf: $p(x,y),{\omega }_i(x,y)\hspace{0.33em}\mathrm{f}\mathrm{o}\mathrm{r}\hspace{0.33em}i=\mathrm{1,2},...,N,\hspace{0.33em}f(x,y),\hspace{0.33em}{T}_s(x,y)$and $T_g(x,y)$.

From a mathematical standpoint, the solution of the formulated problem is not unique unless additional conditions are specified. In the absence of such conditions, the solution constitutes an infinite set of possible fields, from which the physically relevant one must be selected. In the non-stationary case considered here, the uniqueness of the solution is ensured by two elements: the initial condition and the boundary conditions. The initial condition has already been discussed in Section 1. The boundary conditions are defined in the regions of airflow along the galleries adjacent to the goaf and are therefore relatively complex. Their formulation requires the specification of pressure, mass fractions of gas components, and temperature distributions along the coordinate axis parallel to the galleries. For this reason, the boundary conditions are presented in detail in Section 3.

3. Mass and heat transport in the gallery space, boundary conditions, and summary of the mathematical model

This study presents a numerical simulation of the self-heating process of coal residues induced by air penetration within a rectangular region of a mine goaf. The analysed domain is bounded by four sides, forming a rectangular area. Some of these boundaries are adjacent to ventilation galleries through which air flows, while the remaining boundaries are in contact with the surrounding, unmined rock mass. Simulation calculations were performed for two ventilation system configurations (Figure 1 and Figure 2). In the first configuration, the analysed region is separated from the longwall face by a single gallery of width (d_y) (Figure 1). The right boundary of this region, corresponding to the longwall front, coincides with the (Oy) axis of the Cartesian coordinate system. The domain extends along this boundary over the interval (x=0 $\le$ y $\le Y$), where (Y) is the width of the goaf equal to the longwall length. The airflow enters at the point (y=0) the gallery uniformly along its entire width at the inlet boundary. The outlet is located at the opposite end of the gallery (point (x=0, y=Y)). In this configuration, only the boundary parallel to the longwall is in direct contact with the flowing air.

In the second configuration (Figure 2), two galleries intersect at right angles at a junction point (x=0, y=Y). The width of the second gallery is denoted by (d_x - usually d_x$\ne$d_y). At the intersection, two air streams merge, one of which originates from the longwall face. The volumetric flow rates of these streams, which may differ, combine at the junction to form a single airflow. The resulting air stream, containing methane released from the rock mass as well as gaseous products of reactions (A), (B), and (C), is directed towards the outlet gallery. According to the model assumptions, the outlet is located at the corner of the analysed domain, with coordinates (x = X, y = 0), where (X) denotes the length of the goaf region.

Airflow along the galleries is fully forced, as it is generated by the main ventilation system. These flows induce air filtration into the goaf, which in turn promotes the self-heating of coal residues and increases the associated fire hazard. The airflow in the galleries forms the basis for defining the boundary conditions of the governing equations (10), (15), (19), (24), (30), and (34), which describe mass and heat transport within the goaf. To determine these boundary conditions, additional equations are introduced to describe the distributions of pressure, mass fractions of gas components, and temperature in the airflow within the galleries. Owing to the predominantly one-dimensional nature of these flows, a new numerical axis (O$\zeta$) is introduced along the direction of airflow. In the first configuration, the numerical axis coincides with the axis (Oy) aligned with the gallery adjacent to the goaf and varies within the interval (0 $\le$ $\zeta \le$Y). In the second configuration, due to the right-angle intersection of the galleries, the axis (O$\zeta$) is extended in a piecewise manner. Its direction changes at the junction point (x=0, y=Y), where the galleries intersect, forming a continuous coordinate along the airflow path. After the change in direction of (O axis, the variable $\left(\zeta \right)$ spans the interval (Y $\le \zeta \le$X + Y).

In contrast to filtration within the goaf, the physical fields describing airflow in the galleries can be considered predominantly one-dimensional. However, gas flow along the galleries is accompanied by turbulence, resulting in fluctuations of velocity and temperature. These effects lead to turbulent diffusion, which constitutes a significant mechanism of momentum and heat transport along the airflow path. Similarly to the goaf region, the solution of the governing equations for momentum, mass, and heat transport in the galleries allows determination of the transient distributions of the following quantities (denoted by the subscript ($\zeta$), referring to the gallery atmosphere):

√

air pressure differential ($p_{\zeta }^{*}$),

√

mass fractions of individual gas components (${\omega }_{\zeta i}$),

Gas Temperature ($T_{g\zeta }$)

The pressure field ( $p_{\zeta }^{*}$($\zeta$,t)) is described by a differential equation analogous to Eq. (10), formulated for a single spatial variable ($\zeta$ ), and takes the following form:

$${\displaystyle \frac{\partial }{\partial \zeta }}\left(E{\rho }_{\zeta }{\displaystyle \frac{\partial p_{\zeta }^{*}}{\partial \zeta }}\right)+S_{m\zeta }=0$$

(40)

where: E- equivalent to K_x, and K_y, (see equation (15)), ${\rho }_{\zeta }$- local density of the gas mixture in the gallery atmosphere.

Mass and heat transport associated with airflow along the galleries, in which forced convection is the dominant transport mechanism, are described by a system of ordinary differential equations for the dependent variables as functions of the variable ($\zeta$ ), and time (t).

$${\displaystyle \frac{\partial }{\partial \zeta }}\left(v_{\zeta }{\rho }_{\zeta }{\omega }_{\zeta i}\right)+{\displaystyle \frac{\partial \left({\rho }_{\zeta }{\omega }_{\zeta i}\right)}{\partial t}}=S_{\zeta i}$$

(41)

$${\displaystyle \frac{\partial }{\partial x}}\left(v_{\zeta }{\rho }_{\zeta }{c}_{g\zeta }{T}_{g\zeta }\right)+{\rho }_{\zeta }{c}_{g\zeta }{\displaystyle \frac{\partial T_{g\zeta }}{\partial t}}=\dot{Q_{\zeta }}$$

(42)

where, according to the assumptions of the flow model, the gas velocity represents the cross-sectional average (superficial) linear velocity of the flowing medium.

$$v_{\zeta }=-E{\displaystyle \frac{\partial p_{\zeta }^{*}}{\partial \zeta }} m·s^{-1}$$

(43)

The source terms appearing in Eqs. (41) - (43) provide the coupling between the processes occurring within the goaf and the airflow in the galleries. They represent the exchange of mass and heat between these two regions and are defined at the interface between the goaf and the galleries. To formulate these terms, an additional coordinate ($n_{\zeta }$), normal to ($\zeta$ ), is introduced. For both configurations, the coordinate ($n_{\zeta }$) is defined as ($n_{\zeta }=X-x$) in the direction perpendicular to the gallery boundary (for 0$\le$y$\le$Y). In the second configuration, due to the geometry of intersecting galleries, the definition of ($n_{\zeta }=y$) is extended locally to account for the change in orientation of the boundary (for 0$\le$x$\le$X)

$$S_{m\zeta }=-{\displaystyle \frac{K_{\zeta }\rho }{{\delta }_{\zeta }}}{\displaystyle \frac{\partial p^{*}}{\partial n_{\zeta }}},kg/m^{3}s$$

(44)

$${\dot{Q}}_{\zeta }={\displaystyle \frac{1}{{\delta }_{\zeta }}}[v_{n\zeta }\rho {\omega }_i-\epsilon D_i{\displaystyle \frac{\partial \left(\rho {\omega }_i\right)}{\partial n_{\zeta }}}], kg/m^{3}s$$

(45)

$$S_{\zeta i}={\displaystyle \frac{1}{{\delta }_{\zeta }}}[v_{n\zeta }\rho c_g{T}_g-\epsilon {\lambda }_g{\displaystyle \frac{\partial T_g}{\partial n_{\zeta }}}],W·m^{-3}{s}^{-1}$$

(46)

The quantities ${(K}_{\zeta })$, (r), (D_i), (c_g) and (l_g) describe the properties of the gas phase flowing in the galleries in the immediate vicinity of the goaf. The apparent linear velocity of the gas, directed normal to the gallery boundary and towards the gallery, is denoted as ($v_{n\zeta }$). This velocity is perpendicular to the gallery axis and represents the filtration flow between the goaf and the gallery. For both configurations, ($v_{n\zeta }$) is equal to (-v_x) and defined consistently along the interface (for 0$\le$y$\le$Y). In the second configuration, its definition is extended locally to account for the change in geometry associated with the intersection of galleries (v_y for 0$\le$x$\le$X). The width of the galleries is denoted by (${\delta }_{\zeta }$) with the subscript indicating the specific gallery segment (x or y) depending on the value of the coordinate ($\zeta$).

The boundary conditions (Eqs. (44 – 46)) describe the exchange of mass and heat between the porous goaf region and the adjacent ventilation galleries. They are formulated based on Darcy’s law and the continuity of mass and energy fluxes across the interface between the two regions. The normal component of the velocity is derived from the pressure gradient in accordance with the filtration model, while the source terms represent the net flux of mass and energy exchanged between the goaf and the gallery flow. This formulation ensures consistency between the porous medium flow and the one-dimensional flow in the galleries.

In both configurations, the boundaries of the goaf that are in direct contact with the unmined rock mass are also taken into account. It is assumed that no mass or heat transfer occurs across these boundaries in the direction normal to their surfaces (O$\zeta$). This condition is expressed mathematically by the following boundary conditions:

$${\displaystyle \frac{\partial \mathsf{\Phi }}{\partial \zeta }}=0$$

(47)

The generalised source term ($\Phi )$is defined as a set of component functions determined within the model, i.e. ($\Phi \in \left\{p^{*},{\omega }_i,T_g\right.$) for (i=1,2,…,N-1). The conditions given in Eq. (46) apply to all boundaries directly adjacent to the unmined rock mass.

For the system of ordinary differential equations (40 - 42), the boundary conditions at the inlet point ($\zeta =0$) for (x=0, y=0) are of particular importance. At this location, clean air enters the gallery with a specified mass flow rate ($\dot{w_0}$) and atmospheric composition. The inlet air is assumed to consist of nitrogen (${\gamma }_1^0=0.79)$and oxygen (${\gamma }_2^0=0.21)$, respectively. The total pressure (p), and temperature ($T_g^0$), are known from operational measurements. Based on these parameters, the density of the incoming air and the corresponding mass fractions (${\omega }_i^0)$of its components can be determined. Under the assumption of dry air, the remaining components are neglected, i.e. their mass fractions are equal to zero. The boundary conditions at the inlet cross-section ($\zeta =0$) are therefore defined as follows:

• ✓ In the case of equation (40), there is a boundary condition of the second kind of the following form:

$$-E\rho {\displaystyle \frac{\partial p_{\zeta }^{*}}{\partial \zeta }}=\dot{w_0}$$

(48)

• ✓ In the case of equations (41) and (42), only first-order conditions occur at this point:

$${\omega }_{\zeta i}={\omega }_i^0$$

(49)

$$T_{g\zeta }=T_0$$

(50)

Since Eq. (40) is a second-order differential equation with respect to the spatial coordinate ($\zeta$), an additional boundary condition is required to ensure the uniqueness of its solution. This condition is specified at the outlet section of the gallery, i.e. at ($\zeta$ =Y) for the first configuration and at the corresponding terminal point of the flow path for the second configuration ($\zeta$ =Y+X). In both cases, the following boundary condition applies:

$$p_{\zeta }^{*}=0$$

(51)

In the second configuration, special attention must be given to the junction point where two air streams merge ($\zeta$ =Y). As previously described, the combined stream flows towards the outlet section ($\zeta$ =Y+X). In the proposed model, this junction is treated as a point source of mass and heat, representing the inflow of an additional air stream. At this location, the standard source terms defined by Eqs. (44) - (46) are replaced by modified expressions accounting for the additional inflow. In particular, the local mass source term in Eq. (40) increases due to the contribution of the incoming air stream. For the component balance (Eq. (41)), the general formulation given by Eq. (14) is applied, where the inlet composition ($\dot{s_i^0})$ corresponds to that of air supplied from the longwall. The mass source terms for individual components are determined based on the known mass fractions of the incoming air (${\omega }_i^0$ for i=1,2,…,N). For brevity, the explicit expressions are not presented here.

To summarise the proposed mathematical model, the numerical simulation provides approximate solutions for the following scalar fields within the goaf and the adjacent gallery regions:

Pressure Differential Field in the Goaf${(p}^{*}(x,y)$) and Along the Galleries${(p}_{\zeta }^{*}\left(\zeta \right)$), Pa

• Mass fraction fields of individual gas components in the goaf ${(\omega }_i)$ and along the galleries ${(\omega }_{\zeta i})$ for (i = 1,2,…,N), kg·kg⁻¹,
• Temperature field of the gas phase in the goaf ($T_g)$ and along the galleries ($T_{g\zeta })$, K,
• Radially symmetric temperature field within spherical solids in the goaf ($T_s)$, K,
• Carbon mass fraction in the solid phase within the goaf ($f)$, kg·kg⁻¹.

In addition, the model determines the velocity field of the gas mixture within the goaf ($v_x, v_y$), as well as the airflow velocity along the galleries ($v_{\zeta })$.

All of the above-mentioned fields are non-stationary, meaning that their local values vary over time during the development of the self-heating process and potential fire.

4. Numerical Solution of the Equations of the Proposed Model

The initial–boundary problem formulated in the previous sections is strongly non-linear and cannot be solved analytically. Therefore, a numerical approach based on the Control Volume Method (CVM), also known as the Finite Volume Method, was applied. This method was further developed by Patankar (1980) for solving transport problems involving combined convection and diffusion processes [26]. The method is based on the discretisation of independent variables. For steady-state problems, only spatial discretisation is required, whereas for non-stationary processes, discretisation is also performed in time. The numerical scheme used in this study follows Patankar’s approach, which ensures stability and convergence of the solution, particularly in convection–diffusion problems. The convective terms are approximated using the POWER-LAW scheme, which is a modification of the classical upwind method. The spatial discretisation is performed by dividing the computational domain into a finite number of control volumes associated with numerical nodes representing the geometry of the excavations and adjacent roadways. The computational domain was discretised using a structured grid consisting of 3 × 5 control volumes. This grid resolution was adopted as a compromise between computational efficiency and the ability to capture the main features of the analysed processes. Due to the preliminary character of the study, a grid independence analysis was not performed and will be considered in future work. The governing differential equations are transformed into algebraic form by approximating derivatives using finite differences based on the distances between neighbouring nodes (dx, dy, $d\zeta$). For non-stationary processes, time discretisation is introduced by dividing the simulation time into finite increments, and the solution is obtained sequentially at discrete time levels. The resulting solution is discrete, as the values of the physical fields are determined only at the computational nodes. Each algebraic equation represents a balance of the transported quantity over a control volume surrounding a central node and its neighbouring nodes, forming a computational stencil. In a Cartesian coordinate system, control volumes take the form of rectangular elements. To properly account for convective transport, a staggered grid arrangement is employed, in which velocity components (v_x and v_y) are defined at the faces of control volumes. The discrete form of the governing equations for a representative control volume, corresponding to a non-stationary two-dimensional case, is presented below:

$$a_w\left({\mathsf{\Phi }}_W-{\mathsf{\Phi }}_P\right)+a_e\left({\mathsf{\Phi }}_E-{\mathsf{\Phi }}_P\right)+a_s\left({\mathsf{\Phi }}_S-{\mathsf{\Phi }}_P\right)+a_n\left({\mathsf{\Phi }}_N-{\mathsf{\Phi }}_P\right)+\dot{S_P}=b_P{\displaystyle \frac{{\mathsf{\Phi }}_P-{\mathsf{\Phi }}_P^0}{∆t}}$$

(52)

where ($\mathsf{\Phi }$) denotes a generalised physical field, whose spatial and temporal variations describe the transport of an extensive quantity, such as momentum, mass, or thermal energy. The subscripts denoted by capital letters refer to the central node (P) and its neighbouring nodes (W, E, S, N), corresponding to the west, east, south, and north directions, respectively. Lowercase indices denote the locations of the cell faces between neighbouring nodes, where the velocity components are defined and where the corresponding transport coefficients are evaluated. The coefficient (b_P) in Eq. (52) is associated with the control volume surrounding node (P) and represents the storage term: mass in the case of mass conservation equations (Eq. (24)) and heat capacity in the case of the energy equation (Eq. (36)). As can be seen from equation (52), the calculation uses a discrete solution corresponding to the previous time level. The quantity ${(\mathsf{\Phi }}_P^0)$ represents the numerical value of this solution at the central point. The term on the right-hand side of equation (52) represents the temporal accumulation and is evaluated using a backward difference approximation of the time derivative, which ensures numerical stability for non-stationary simulations. At the initial time level, the values of ${(\mathsf{\Phi }}_P^0)$ are defined by the initial conditions. The solution is then advanced in time by successively updating the field values from the previous time level to the current one. This sequential procedure involves solving the system of algebraic equations at each time step. The calculations may be continued for a prescribed number of time steps or until a steady-state solution is reached, if such a state exists. The numerical solution at each time step is obtained according to the following computational sequence:

• Solution of the filtration equation (15) to determine the pressure differential field ($p^{*})$ in the goaf, including adjacent galleries.
• Calculation of the gas velocity field in the goaf based on Darcy’s law (Eq. (17)), using the previously determined pressure field.
• Determination of the carbon mass fraction in the solid phase (f) by solving Eq. (9) within the goaf region.
• Solution of the system of equations (24) to obtain the mass fractions (${\omega }_i)$ of gas components (excluding methane) for i=1,2,…,N-1 in the goaf and adjacent galleries.
• Determination of the methane mass fraction field for working areas including adjacent galleries using Eq. (25).
• Solution of the energy equation to obtain the temperature field (T_g) of the gas phase in the goaf and adjacent galleries.
• Numerical solution of Eq. (37), together with boundary conditions Eq. (38) and (39), to determine the temperature distribution within spherical solids.

Due to the strong non-linearity of the governing equations, the above computational steps are performed iteratively at each time level. In addition, non-linear source terms in Eq. (24) are linearised using standard procedures described by Patankar (1980) [26].

To ensure numerical stability and accuracy, a variable time step is applied. Its value is determined adaptively based on the local rate of combustion reactions. Specifically, the time step is selected such that the decrease in oxygen mass fraction within the most reactive control volume does not exceed a prescribed threshold (Δω_min = 0.02). This value was chosen based on numerical experiments as a compromise between computational efficiency and convergence of the solution. Smaller values lead to significantly increased computation time, while larger values reduce the accuracy of capturing local oxygen depletion.

It should be noted that detailed analysis of grid independence, convergence criteria, and solver parameters was not included in the present study, as the primary objective was the formulation of the mathematical model and its physical consistency. These aspects will be addressed in future work.

5. Example Calculations and Discussion of Results

Based on the developed mathematical model and the implemented simulation program, exemplary results are presented to demonstrate the capabilities of the proposed approach. The calculations were performed for the longwall ventilation systems shown in Figure 1 and Figure 2. Direct validation of the obtained results against in-situ measurements within the goaf is not feasible due to technical limitations. The goaf constitutes an inaccessible region, where installation of appropriate measurement equipment is not possible. Therefore, only indirect validation can be carried out based on measurements at the boundaries of the analysed area. It should be emphasised that the primary objective of this study was not quantitative validation of the model, but rather the development and demonstration of a physically consistent simulation framework. Detailed validation studies are planned for future work. The numerical results were obtained in discrete form and subsequently visualised using the Surfer software, which enables representation of results in the form of contour maps. In the following figures, only selected fields are presented, namely methane (CH₄), oxygen (O₂), carbon monoxide (CO), and gas temperature distributions. Figure 3, Figure 4, Figure 5 and Figure 6 show contour maps obtained for the ventilation system in the U configuration (Figure 1), including:

–

methane concentration (Figure 3),

–

oxygen concentration (Figure 4),

–

carbon monoxide concentration (Figure 5),

–

gas temperature (Figure 6).

The coordinate system used in these figures corresponds directly to that shown in Figure The vertical axis (0–250 m) represents the longwall face length, while the horizontal axis corresponds to the goaf area formed as a result of mining operations. At the location corresponding to 250 m, fresh air enters the longwall panel, initiating airflow penetration into the goaf. This results in a decrease in methane concentration near the inlet region (Figure 3). In contrast, at the outlet (0 m), an increase in methane concentration is observed due to gas accumulation and outflow. A similar trend is observed in (Figure 4), where its concentration increases near the inlet and penetrates approximately 75 m into the goaf. This oxygen participates in oxidation reactions of residual coal, leading to self-heating processes. As a consequence, carbon monoxide is generated (Figure 5), and an increase in gas temperature is observed (Figure 6), indicating the development of self-heating phenomena.

Due to the lack of direct access to the goaf, the presented distributions cannot be validated directly. However, comparison with operational mine measurements at the inlet (250 m) and outlet (0 m) indicates that the qualitative trends of gas concentration distributions are consistent with practical observations.

Figure 7, Figure 8, Figure 9 and Figure 10 present analogous contour maps for the Y ventilation system (Figure 2), including:

–

methane concentration (Figure 7),

–

oxygen concentration (Figure 8),

–

carbon monoxide concentration (Figure 9),

–

gas temperature (Figure 10).

The coordinate system corresponds to that shown in Figure 2. The vertical axis (0–250 m) represents the longwall face length, while the horizontal axis corresponds to the goaf region adjacent to the return airway. In this configuration, fresh air enters at 250 m, while additional air mixing occurs near the outlet (0 m). As a result, the airflow penetration into the goaf is significantly more intensive than in the U system. Consequently, higher oxygen availability within the goaf leads to intensified oxidation processes. This results in different spatial distributions of methane, oxygen, and carbon monoxide compared to the U system. In particular, the Y system is more prone to self-heating phenomena, especially in areas with residual coal and geological disturbances.

For the Y-type ventilation configuration, the presence of an additional return airway provides more boundary locations at which future indirect validation of the numerical results may be performed using operational measurements.

Although direct validation, especially in first case scenerio, is not feasible due to the inaccessibility of the goaf, the obtained results are consistent with qualitative observations reported in industrial practice and previous studies. Above that, the presented simulations demonstrate that the proposed model can be a valuable tool for designing and optimising fire prevention strategies in goaf areas, both at the planning stage and during ongoing mining operations. Future work will focus on quantitative validation using indirect measurement data.

6. Conclusion

The present study developed a mathematical and numerical framework for analysing transient coupled heat and mass transfer processes associated with coal self-heating in a reactive porous medium. The model simultaneously accounts for gas filtration, species transport, heat transfer between the gas and solid phases, heat generation due to heterogeneous and homogeneous chemical reactions, and the resulting evolution of temperature and gas composition. A key feature of the proposed approach is the incorporation of spatially variable permeability based on in-situ mining data, enabling heterogeneous flow and transport conditions within the longwall goaf to be represented.

Numerical simulations performed for U-type and Y-type ventilation configurations demonstrated that the ventilation layout significantly affects gas flow, oxygen penetration, temperature development, and the spatial distribution of zones susceptible to coal self-heating. The Y-type ventilation system promoted deeper oxygen penetration into the goaf, thereby creating conditions that may favour more extensive coal oxidation and self-heating compared with the U-type configuration. These results demonstrate the strong coupling between ventilation-induced gas flow, species transport, chemical reactions, and thermal processes in the porous goaf medium.

From a thermal engineering perspective, the results show that reliable analysis of coal self-heating requires the simultaneous consideration of fluid flow, species transport, reaction heat generation, and heat exchange between the gas and solid phases. Spatial variations in permeability affect local transport conditions and therefore influence oxygen availability, heat accumulation, and temperature development within the porous medium. The proposed modelling framework provides a tool for analysing these coupled phenomena under transient operating conditions and for evaluating the influence of ventilation configurations on thermal hazard development.

The present model is subject to several limitations. The computational domain was discretised using a relatively coarse structured grid, and a systematic grid-independence analysis was not performed. In addition, the predicted temperature and gas concentration fields have not yet been validated against dedicated experimental or field measurements. The model does not explicitly account for moisture transport and related evaporation and condensation processes, which may influence heat transfer and coal oxidation behaviour.

Further research should therefore focus on experimental and field validation of the predicted temperature and gas concentration fields, systematic grid-independence and sensitivity analyses, and extension of the mathematical formulation to include moisture transport and phase-change phenomena. Application of the modelling framework to other reactive coal-bearing porous systems would additionally require appropriate modification of the geometry, transport properties, initial conditions, and boundary conditions.