A computational model to simulate self-heating ignition across scales, con ﬁ gurations, and coal origins

Self-heating of fuel layers can trigger ignition when the temperature of the surroundings is su ﬃ ciently high. Self- heating ignition has been a hazard and safety concern in raw materials production, transportation, and storage facilities for centuries. Hot plate and oven-basket experiments are the two most used lab-scale experiments to assess the hazard of self-heating ignition. While extensive experiments have been done to study this phenomenon, modelling of the experiments is substantially lagging behind. A computational model that can accurately simulate self-heating ignition under the two experimental con ﬁ gurations has not been developed yet. In this study, we build such a model by coupling heat transfer, mass transfer, and chemistry using the open-source code Gpyro. Due to the accessibility of large amount of experimental data, coal is chosen as the material for model validation. A literature review of the kinetic parameters for coal samples from di ﬀ erent origins reveals that there is a compensation e ﬀ ect between the activation energy and exponential factor. Combining the compensation e ﬀ ect with our model, we simulate 6 di ﬀ erent experimental studies covering the two experimental con ﬁ gura- tions, a wide range of sample sizes (heights ranging from 5mm to 126mm), and various coal origins (6 coun-tries). The model accurately predicts critical ignition temperature ( T ig ) for all 24 experiments with an error below 7°C. This computational model uni ﬁ es for the ﬁ rst time the two most used self-heating ignition experiments and provides theoretical insights to understand self-ignition for di ﬀ erent fuels under di ﬀ erent conditions.


Introduction
Self-heating is the temperature rise tendency of a material due to heat generated by exothermic processes taking place within the body of material [1]. Self-heating can lead to ignition when the rate of heat generation is greater than the rate at which the heat dissipated to the surrounding. In many industrial processes, fuel dust is often produced as main or by-product. Fuel dust can unintentionally accumulate or be intentionally stored in layers. A fuel layer has self-heating propensity, because exothermic low-temperature oxidation can occur between solid phase fuel dust and the oxygen diffusing from ambient atmosphere inside the fuel layer. Self-heating of a fuel layer can lead to ignition when the surrounding that the layer is exposed to is sufficiently hot. Self-heating ignition of coal layers has been a safety concern in fuel production, transportation, and storage facilities for centuries.
Hot plate experiment and oven-basket experiment are the two most commonly used experimental techniques to assess the self-heating and ignition hazard of fuel dusts. The hot plate experiment is particularly relevant for the situation where dust can accumulate on hot surfaces [2]. In this experiment, a uniform thin layer of dust is placed on a hot plate that is maintained at a constant temperature, while the top face of the layer is exposed to the ambient air with side faces framed by a metal ring [3,4]. The hot plate temperature is increased between different experimental runs until ignition is observed. The experiment aims to measure the critical temperature of the hot plate that triggers the ignition of a dust layer. This critical temperature is referred as minimum ignition temperature [5].
The oven-basket experiment is used to classify self-heating substances for the transport of dangerous goods [6]. In the oven-basket experiment [4,7], the bulk dust is packed inside a basket which is then hung in an isothermal oven with constant temperature. The oven temperature is increased between runs until ignition occurs. This experiment aims to measure the critical oven temperature that triggers the ignition of the dust bulk. This critical temperature is referred as selfignition temperature [8].
It is seen that both experiments aim to measure the critical surrounding temperature (either hot plate or oven temperature) that triggers ignition. This temperature is one of the most important indicators to assess the self-heating and ignition hazard for fuel layers. For the convenience of notation, the minimum ignition temperature in hot plate experiment and the self-ignition temperature in oven-basket experiment are denoted uniformly as T ig in this paper .
Both types of experiments have been conducted extensively to evaluate the self-heating and ignition hazard of various fuels, such as coals [7,[9][10][11][12][13], carbon-rich soils [14], shale rocks [15,55], and biomass [16][17][18][19]. Janes et al. [2] tested 14 different materials and reported a correlation between the two experiments. Based on the correlation, a procedure is proposed to estimate the minimum ignition temperature of a dust layer using the oven-basket experiment when it is not possible to carry out hot plate experiment.
To gain more insights into these two experiments, theoretical studies have also been conducted. Frank-Kamenetskii [20] established the fundamental self-heating ignition theory for porous solid materials. He proposed a dimensionless number, the Frank-Kamenetskii (F-K) parameter, to determine the critical ignition condition. Bowes and Thomas [1,21] advanced Frank-Kamenetskii's theory by taking the convection on the boundary into consideration and obtained a corrected F-K parameter which includes the Biot number. From the perspective of these fundamental theories, the self-heating ignition that occurs in these two experimental configurations is the same phenomenon and the mechanisms (heat transfer, mass transfer, and chemistry) that control self-heating ignition do not change with experimental configurations. The main difference between these two experiments is the boundary condition: Oven-basket experiments have a symmetrical boundary condition whereas hot plate experiments have an unsymmetrical one. By modifying the boundary condition, the ignition criterion for two different experimental configurations can be analytically determined within the same theoretical frame [1]. However, these theories assume steady-state and mainly focus on the criticality. The transient behavior of self-heating ignition, such as the evolution of temperature, the time to ignition, and transient temperature profile cannot be obtained. To better understand this phenomenon, we need computational models that can conduct transient simulation.
So far a number of computational models have been built for various types of materials, covering coal [9,[22][23][24][25][26][27], metal [28,29], and biomass [30,31]. However, all of these computational models were developed for a specific experimental configuration, either for hot plate configuration [22][23][24]29,31] or for oven-basket configuration [9,25,26,30]. This split implies that the physics of self-heating ignition depend on the experimental configuration. However, the fundamental mechanisms of self-heating ignition should be the same for two experimental configurations as explained by Frank-Kamenetskii theory. Thus it is possible to develop a computational model that can simulate self-heating ignition of dust layers under different configurations. The objective of this study is to build such as a model. As one of the most important solid fuels, coal has received the most attention in the past. Self-heating and ignition hazard of coal has been evaluated with various sources and under different configurations [7,[9][10][11][12]. Due to the availability of a large amount of experimental data, coal is chosen as the material to validate the model. Our model is first validated against two high-fidelity experiments (one hot plate experiment and one ovenbasket experiment). Then six different experimental studies (24 experiments in total) covering both experimental configurations, a wide range of sample sizes, and various coal origins are simulated using this validated model.

1D governing equations
Self-heating ignition of fuel dust is driven by the low-temperature heterogeneous oxidation occurring between solid phase fuel and gas phase oxygen in a porous media system. Oxygen diffusion and heat transfer inside the porous media system are the two key processes that control this phenomenon [2]. To simulate these two processes, the open-source code Gpyro is used here [32].
Gpyro is a generalized code that can be applied to simulate pyrolysis and smouldering of reactive solids, especially porous media. Uniquely, Gpyro has a convective-diffusive solver to simulate gas diffusion from the ambient into the porous media, providing the capability to calculate the transient composition of gaseous species (reactants and products) at different locations [32]. Combining this solver with the conservation equations of mass, energy, and species in both the condensed phase and gas phase, the distributions of temperature and species inside a thermally stimulated porous media can be calculated. Fig. 1 displays the computational domain and boundary conditions for the two experimental configurations. The two configurations can be simulated by the identical set of governing equations (Eqs. (1)-(6)), which impose conservations of (1) mass, (2) species and (3) energy in condensed phase as well as (4) mass, (5) species, and (6) momentum (Darcy's law) in gas phase. Eqs. (3) and (6) are two key equations in the simulation. They are responsible for the calculation of heat transfer and oxygen concentration. Subscript i and j refer to solid species and gas species respectively. In this work, i can be coal/ash and j can be oxygen/nitrogen/gaseous products. All the other symbols are explained in the nomenclature.
In the model, the averaged properties in each cell are calculated by weighting mass or volume fractions. All gaseous species have unit Schmidt number, and equal diffusion coefficient and specific heat. In this study, the properties of coal are assumed to be independent of temperature. The coal sample is assumed to be dry and the effect of water adsorption and desorption is not considered in the simulation. The details of the mathematical formulation of Gpyro can be found in [32].
The two different experimental configurations can be simulated using the same governing equations, but they do have different boundary conditions.
In the hot plate configuration, the sample is a circular-shaped thin layer, as shown in Fig. 1(a). Since the thickness of the layer is significantly smaller than its diameter, the radial heat transfer from the side can be neglected [23]. Therefore, hot plate experiments can be simulated with a 1D model.
The boundary conditions for hot plate experiments are as follows: the top face (z = L) is open to the air atmosphere; the bottom face (z = 0) is impermeable to mass flow, which means no mass flux and diffusion of gas species though the bottom face. The boundary conditions can be defined by Eqs. (7)- (12). Eqs. (7) and (8) are the boundary conditions for the energy conservation equation. The temperature of the hot plate is a fixed value in the experiments [3], and this constant temperature is, therefore, set as the boundary condition for the bottom (Eq. (8)). The heat transfer between the sample and hot plate surface can be calculated through the temperature gradient at z = 0. Both convective and radiative heat transfer are considered on the top face (Eq. (9)) and the emissivity (ε) of coal is 0.8 [33]. Eqs. (9)-(12) are the boundary conditions for the conservation equations of gas momentum and gas species. h m is mass transfer coefficient, which can be calculated by h m = h c /C g according to mass-heat transport analogy [32]. C g, the specific heat capacity for gas species, is assumed to be 1100 J/kg-K as a constant in the simulation [32].
At t = 0, the entire coal sample is unreacted (Y c = 1). The initial temperature for both coal and air inside the sample is the same as the ambient temperature (T ∞ = 20°C). The whole computational domain, including the pores inside the dust, has the same initial gas composition as ambient air (Y O = 0.22/Y N = 0.78).
In the oven-basket configuration ( Fig. 1(b)), the coal dust sample has a cubic or equi-distant cylindrical (L = diameter of base) geometry with all sides open to air atmosphere. Since the height of the sample is equal to the length of its base, the heat transfer through the side face should be taken into account. Thus this is a multi-dimensional problem. However, it can still be reduced to a 1D problem as long as the heat losses through the side faces are considered in the model. Here, the volumetric heat loss coefficient (h vl ) is incorporated into the energy conservation equation to account for these heat losses as shown in Eq.
(3). This coefficient can be used to describe the heat transfer process between sample center and ambient air [34]. It can be regarded as the reciprocal of effective heat resistance, which consists of two parts: natural convection (1/h face ) at side faces and conduction (L/2k c ) from the sample center to the side face. Then the effective heat transfer coefficient is converted into the volumetric term h vl by multiplying the ratio of side area to volume. Since oven-basket experiments use simple basket geometries, either cubic or equal-distant cylinder, the ratio of side area to volume is 4/L. Thus h vl can be calculated by Eq. (13) [35], where the coefficient of natural convection is calculated using Eq. (14) [36]. For the hot plate configuration, h vl is equal to 0 since there are no heat losses through side faces In the oven-basket experiment, boundary conditions for top and bottom faces are the same. Both faces are open to the air atmosphere and have the same boundary conditions as the top face in the hot plate experiment. They can be described by Eq. (15)- (19).
At t = 0, the entire sample is unreacted coal (Y c = 1). The initial temperature for coal and the air inside the sample is 20°C, while ambient temperature is set to T oven. The whole computational domain, including the pores inside the dust, has the same initial gas composition as the ambient air (Y O = 0.22/Y N = 0.78).

Chemical kinetics of coal self-heating ignition
The exothermic reaction during coal self-heating ignition can be described by a global one-step kinetic scheme [23,37]. According to [23], 1 kg of reactive component in coal reacts with approximately 2 kg of oxygen to generate 3 kg of gases with a certain amount of ash left. The stoichiometry of this scheme can be determined accordingly as shown in Eq. (20). It should be pointed out that for heterogeneous reaction, the stoichiometry is usually given on the mass basis rather than molar basis. This is because the solid reactant in the heterogeneous reaction is usually a mixture of different components and the exact molecule composition is difficult to specify experimentally. It is a typical way to use one or several global reactions along with a few lumped species to describe the heterogeneous chemistry process [38]. Thus, the mass-based rather than molar-based stoichiometry is used to quantify these lumped species in the heterogeneous reactions as shown in previous works [32,38].
The reaction rate is assumed to have an Arrhenius form [23,37]. The reaction rate is assumed [1,39] to depend on both the mass fraction of coal (Y c ) and oxygen (Y O ) as shown in Eq. (21). n c and n O represent the reaction order of coal and oxygen respectively.
a Only E is optimized in the validating simulation to account for the uncertainty.

. Parametrization
A base case is first simulated following the hot plate experiment conducted by Park et al. [3]. In this experimental work, four different layer heights (6.4 mm, 12.7 mm, 19.1 mm, and 25.4 mm) are experimented and the hot plate temperature was increased by 5-10°C for each run until ignition was observed. In the experiments, almost all thermo-physical parameters required in the simulation were measured with accuracy and high quality. Specifically, the experiment conducted with the layer height of 12.7 mm is chosen here for model validation. This is because the minimum ignition temperature measured in the case of 12.7 mm has a much smaller uncertainty (5°C) than the other three cases (10°C), therefore providing the highest accuracy. Table 1 lists the essential parameters used in the simulation. Most of the physical parameters were measured or calculated in the original experiment [3], except for heat capacity (c), which is obtained from [40]. The original experiment estimated chemical parameters (E and HA Δ ) as well. However, these chemical parameters were obtained based on a simple kinetic scheme without considering the influence of oxygen concentration (n O ). In addition, H A Δ and were measured as a couple ( HA Δ ). The coupled value is not suitable to be input in the simulation. Thus, here we use the kinetic parameters obtained from [37], which reported a more complete set of kinetic parameters for the same Pittsburgh seam coal. To account for the influence of mass fraction of coal (as shown in Eq. (21)), n c is added to the scheme and assumed to be 1 [23] (As discussed latter in the sensitivity study this assumption does not affect simulation results). It should be noted that the value of E (activation energy) finally used in the validating simulation has been optimized, since the uncertainty of E can lead to a significant difference in the simulation results. A detailed discussion on the optimization process is given in the Section 3.1.2.

Sensitivity analysis
Simulation requires the values of various parameters as inputs. Most parameters can be measured from experiments, but uncertainty always exists. Because of the measurement uncertainty, it is necessary to quantitatively analyze the influence of different parameters. For the parameters that have significant effects on simulation results, their uncertainties should be considered in the simulation.
For the base case simulation, we first conducted a tentative calculation using literature values for all the parameters. The tentative simulation over-predicted the T ig by around 40°C (the green star in Fig. 2(a)). Thus, a sensitivity analysis is conducted to find the most sensitive parameter(s), whose uncertainty could be the source of the over-prediction.
Two analysis approaches are used here: scatter-plot analysis [41] and one-at-a-time (OAT) analysis [41]. Scatter plot is a graphical sensitivity analysis approach, which can demonstrate how the output changes with the parameter of interest while the other parameters fixed at their base values. Fig. 2(a) is the result of scatter plot analysis. Every single line displays how the simulated T ig changes with the variation of a single parameter while the other parameters are kept at their base values (i.e the literature values in Table 1). For the convenience of comparison, the variation of different parameters has been normalized by their base values. It can be observed that E has the dominant influence on T ig . Varying it by only 7%, from 89.9 kJ/mol to 83.6 kJ/mol, the experimental T ig can be accurately predicted. For all other parameters, the experimental T ig cannot be obtained even with an approximately 100% variation (The variations of k and ϕ are less than 100% to avoid unphysical values). Therefore, we choose E as the parameter to be optimized as a consideration of measurement uncertainty and the 7% variation is within ∼10% uncertainty range in the measurement of E [42].
One-at-a-time analysis is the technique to quantify the sensitivity levels of different parameters using sensitivity coefficient. Sensitivity coefficient is the derivative of output (y) with respect to the parameter of interest (X). Since the parameters compared here have different units, the dimensionless sensitivity coefficient s is used. It can be calculated by Eq. (22) [41].
In this study, T ig is chosen as y and the 10 parameters listed in Table 1 are treated as X i . The OAT analysis results are summarized in Fig. 2(b) and the absolute value of s i is presented for the convenience of comparison. Agreeing with the scatter-plot analysis results, E is the determining parameter with a sensitivity coefficient 10 times larger than any other parameter. A, ΔH, and n O play similar secondary roles. The only exception in chemical parameters is n c, which has almost no influence on T ig. This is because at the critical point of self-heating ignition, the depletion of coal is not significant yet and Y c is still close to 1 [1]. In this circumstance, the influence of n c is negligible as the calculation of reaction rate (Eq. (21)) is almost independent of n c when Y c ≈ 1. In terms of physical parameters, ρ and k have the secondary influences similar to A, ΔH, and n O , while the effects of h c , ϕ and c are negligible.
An intuitive explanation for the sensitivity analysis result can be given based on the energy balance equation at criticality. The situation at the critical point can be considered as a steady state [1]. A quasisteady energy conservation equation, Eq. (23), is therefore derived from Eq. (3) by cancelling the time-variant term and substituting Eq. (21) into it.
Eq. (23) shows that at the critical condition, the heat generated from the oxidation reaction (left side) balances with the heat dissipated through diffusion (right side). The ignition threshold (T ig ) would reduce as the heat generation increases or heat dissipation decreases, which can be achieved by the variation of parameters. Among these parameters, E should have the largest influence. This is because a slight variation of E would lead to an exponential change in heat generation and hence largely influence the temperature (T ig ) needed to trigger the ignition. The effect of E is therefore much larger than other parameters.
The significant impact of E on the T ig was reported in both experimental [6] and numerical [43] studies before. As discussed above, this is due to the exponential effect of E in Arrehnius reaction rate. Since the Arrhenius equation is most commonly used to describe reaction rate for self-heating ignition in the numerical study, this implies the need to improve the measurement accuracy of E. The present measurement uncertainty of E is 5-10% [42], a relatively large value that can lead to significant differences in the prediction of T ig as shown in Fig. 2. Therefore, E is optimized here as a consideration of measurement uncertainty.

Base case validation
In this study, a grid/time step independence study is first conducted to ensure the simulation results are not dependent on grid size or time step. For the heat conduction driven problem investigated here, the refining of grid size and time step should obey Eq. (24) to avoid unphysical solutions [44]. This means that grid size and time step should be refined simultaneously: Whenever grid space is reduced by a factor of n, the time step should be reduced by a factor of n 2 . As shown in Fig. 3, the solution tends to converge at = z Δ 0.1mm ( = t Δ 0.01s), whereas the computing time keeps increasing exponentially after that point. Thus, = z Δ 0.1mm and = t Δ 0.01s are the independence points and the following simulations in this study are conducted at this grid size and time step.
The validating simulation is conducted with the optimized parameters listed in Table 1. The comparison between simulation and experiment is displayed in Fig. 4. Ignition is considered to occur when the temperature of sample is 50°C higher than the hot plate temperature [45]. The time to ignition is also defined at this point, i.e. the time when the temperature of the sample is 50°C over the hot plate temperature. In this work, simulations stop when the sample temperature reaches the ignition criterion. As the focus of the paper is the onset of ignition, the post-ignition scenario has not been considered in the present model.
As shown in Fig. 4(a) and (b), the temperature curves flatten out when hot plate temperature (T hp ) is 210°C. When T hp is increased to 215°C, temperature curves significantly increase after 3000 s and finally exceed 265°C, the ignition criterion in this specific case. The model accurately predicts T ig = 215°C. Besides, the model also properly simulates the evolution of temperature at different heights. The predicted temperature histories agree well with the experimental date with a ∼10% predicting error. In the experiment, ignition occurs at around 3000 s, which is also captured by the simulation. Fig. 4(c) and (d) compare the predicted transient profiles of temperature and mass fraction of oxygen (Y O ) with the experimental data. In the subcritical condition (T hp = 210°C), the profile of temperature demonstrates a parabolic-like shape, where the highest point is close to the bottom. The profile of Y O approximately mirrors the temperature profile with the lowest point near the bottom. This is because at a lower location sample dust has a higher temperature, leading to a higher reaction rate and hence more consumption of oxygen. After 3000 s, the profiles of both quantities don't change much anymore with time. At this point, the heat generation rate and heat dissipation rate start to balance with each other. The diffusion of oxygen also balance with the consumption of oxygen. It can be regarded as a steady state. In the supercritical condition (T hp = 215°C), both profiles change  significantly after 3000 s. This is because at that time ignition has been triggered and the exothermic reaction starts to accelerate. Due to the increasing reaction rate, the consumption rate of oxygen starts to exceed the diffusion rate, leading to the significant decrease of Y O . At the same time, heat generation rate increases significantly and exceeds heat dissipation rate, resulting in temperature increase after ignition. The predicted temperature profile shows good agreement with experimental data. The largest predicting error is around 30°C. In addition, it is seen that the ignition tends to occur near the bottom of the sample and the hot spot (the highest temperature point) is located around z = 3-4 mm, also agreeing with classical theory [46].

Validation against oven-basket experiment
The model is also validated against an oven-basket experiment. Here the experiment conducted by Nugroho et al. [11] is used for model validation. In the experiment, the sample is a cubic with a length of L = 50 mm. The oven temperature was increased by 1°C for each run until ignition was observed.
As discussed in Section 2, the model for the oven-basket experiment is built on the model for the hot plate experiment, only with the volumetric heat loss coefficient introduced to account for heat losses through side faces. The underlying physical mechanisms of the two models are the same. Thus, sensitivity analysis result obtained and the parametrization methodology used in the base case simulation of the hot plate experiment can be applied here. Table 2 lists the essential parameters used in the simulation. Most parameters can be obtained from the original experiment. For the same reason discussed in Section 3.1.2, E is optimized to account for the measurement uncertainty. For the unreported parameters, their base case values in Table 1 are used here. The volumetric heat loss coefficient is calculated using the method mentioned in Section 2.
The evolution of temperature was only monitored at the center of the sample in the experiment. Therefore, the comparison between simulation and experiment is conducted for the temperature history at z = 25 mm as shown in Fig. 5. In oven-basket experiments, ignition is considered to occur when the sample temperature rises 60°C above the oven temperature [8,47].
The sample ignites when T oven is increased to 127°C. The model accurately predicts T ig , the most important index to assess the selfheating and ignition hazard. Since the experiment [11] only reported the temperature trend after the center temperature exceeds the oven temperature (exceeding point), the temperature evolution before that point cannot be observed in the Fig. 5. Thus the comparison between simulation and experiment is conducted based on the data after the exceeding point. The model predicts that ignition occurs at approximately 65 min (i.e the time when the center temperature is 60°C over T oven ). This prediction is only 5 min later than the experimental observation. Considering the significantly longer induction delay in ovenbasket experiment than hot plate experiment, this prediction error is satisfactory. In terms of the temperature trend, the agreement between simulation and experiment is reasonable. The temperature curves in both super-critical and under-critical conditions are slightly underpredicted. The under-prediction is within 15°C in the case of T oven = 127°C and 30°C in the case of T oven = 126°C. Considering the simplification made to reduce the 3D problem to a 1D model, the reasonable agreement can be regarded as a validation of our model.

The compensation effect between E and A across coals
As shown in Section 3.1.2, kinetic parameters have a large influence on the simulation results. In previous computational studies, an experiment to obtain kinetic parameter was always conducted before the simulation. Fig. 6 summarizes the E and A obtained from experiments for different coal samples. Most of the experimental data (except for the experimental data from Reddy et al. [12]) are directly obtained from oven-basket experiments, which are the most commonly used technique to obtain E and A for self-heating ignition at present. We can find that there exists a linear correlation between E and ln(A) across different coals. This linear correlation is referred as the compensation effect [48].
A similar compensation effect has been reported by Nugroho et al. [42] previously, but in their study the two correlated parameters were E and ln( H |Δ |A). In the correlation presented here, H |Δ | is removed to demonstrate a purely chemical compensation effect between only E and A. In this way, a higher correlation coefficient can be obtained with even more experiments included.
The compensation effect is a characteristic pattern of rate behavior that can be found across similar reactants [49]. The compensation effect has been found in the pyrolysis of pitch [50] and cellulosic materials [51,52] through TGA experiments, where cummings kinetic model [53] is used to obtain E and A. Eqs. (25) and (26) can be used to interpret the compensation effect between E and ln(A) [48]. η is the reaction rate constant in Arrhenius equation (Eq.(21)). The occurrence of compensation effect here might suggest the existence of some characteristic temperature T iso , around which the exothermic reactions occurring during self-heating ignition of different coals have a same reaction rate constant, which can be referred as η iso .
In practice, the compensation effect between E and ln(A) means the two parameters are correlated and independent. The reactive variation among different coal samples can be described by the variation of a single kinetic parameter (the activation energy).
It is interesting to note that the E and A used in the base case validation are in the middle range of the correlation line. They can be regarded as an average representation of kinetic parameters for different coal samples investigated here.

Prediction of T ig for different coal self-heating ignition experiments
Combining the findings from Sections 3 and 4, we build a model to predict the critical ignition temperature (T ig ) across experimental configurations, scales, and coal origins. Six different experimental studies (24 experiments in total) covering two experimental configurations, a wide range of sample sizes (heights ranging from 5 mm to 126 mm) and various coal origins (6 different countries) are simulated to test the generality of the model. Such an experimental dataset can be regarded as a representative of most lab-scale coal self-heating ignition experiments.
The parametrization in the simulations is largely based on the base  Fig. 6, the coal from different origins tend to have a compensation effect between E and A. The base case values (the optimized value in Table 1) are used here, as they are the average representation for different coals studied here. In addition, n O and n c, the kinetic parameters with secondary influence on the simulation (as shown in the sensitivity analysis of Section 3.1) are usually not measured or reported in the experiments. The information for the precise parametrization for these two parameters is not enough. Taking all these factors into consideration, here we only vary E to account for the kinetic difference between different coals, while the values of A, n O , and n c are kept fixed as the base case values. This allows to reduce the number of variables and therefore largely simplify the parametrization of kinetics. For the physical parameters, some of them were reported in the corresponding studies as listed in Table 3. For the unreported ones (c, ϕ, h c ), their base case values (listed in Table 1) are used since they do not vary much with experiments and have secondary or negligible influence on simulation results as shown in the Section 3.1.2. For oven-basket experiments, h vl is calculated to estimate the heat losses through side faces using the method mentioned in Section 2. Since its value varies Fig. 5. Comparison of predicted (solid) and experimental [11] (dash) temperature histories at the center of the sample (z = 25 mm) for oven-basket experiment. Fig. 6. The compensation effect between E and ln(A) across coal origins; The legend refers to the source of experimental data [7,[9][10][11][12]27] and the origin countries of coal sample.  with the sample size, the value range for each study is given in Table 3. The comparison between simulations and experiments is displayed in Fig. 7. In the simulations, the hot plate temperature and oven temperature are increased by 2°C between two consecutive runs until ignition occurs. Error bars are presented for the studies that reported them. It is seen that simulations give a satisfactory prediction of T ig for all 24 experiments coming from six different studies, which cover a wide range of sample sizes (from L = 5 mm to L = 126 mm). The predicting error is less than 7°C. T ig decreases as the simple size increases for both hot plate and oven-basket experiments. The decrease of T ig is more than 45°C as the height of the sample (L) increases from 20 mm to 110 mm in the oven-basket experiments (the red curve in Fig. 7). Physically, this is because as the sample size increases the heat dissipation rate becomes smaller and therefore less heat generation from oxidation is required to trigger ignition, leading to the decrease of T ig .
As shown in Figs. 4 and 7, the present model is able to simulate oxygen diffusion inside the sample and successfully captures the effect of height on heat transfer. Oxygen diffusion and heat transfer are two key processes that determine self-heating ignition. This fact does not change with scale. The authors think that as long as the two key processes are correctly modelled (as the present model does), real-size scenario can be simulated. Due to the lack of real-size experimental data in the literature, the model is validated against the largest size experiment we could find (L = 126 mm). The validation against realsized experiments will be the interest of future work.
It should be noted that two studies (Schmidt et al. [25] and Wilén et al. [54]) simulated here didn't report kinetic parameters in their original work. Thus, prior to the simulations, we didn't know whether the kinetics of the coals used in these two studies comply to the compensation effect shown in Fig. 6. However, the experimental data in the two studies are still well predicted. This can be regarded as a blind test of the compensation effect across different coals and indicates the compensation effect might not be just limited to the coals summarized in Fig. 5.
The good agreement with the wide-ranging experimental dataset tests the generality of our model and demonstrates the model is able to predict T ig for coal self-heating ignition across configurations, scales, and coal origins.
As discussed in Section 3.1.2, activation energy has a significant influence on T ig . On the other hand, at present the measurement uncertainty of E is relatively large. To improve the predicting accuracy of the computational model, here we propose a possible procedure to use a simple experiment for model calibration: Use the base values (listed in Table.1) of kinetic parameters (E, A, n c , n O ) to conduct tentative simulation. Then adjust E, with the other three parameters fixed, to match the experimental T ig obtained in an experiment of a single size. Then use the calibrated model to estimate T ig for the specific coal at other sizes. In fact, that is how the simulations in Fig. 7 are conducted.
At present, evaluating self-heating ignition hazard for a specific coal has an intensive requirement on the quantity and/or quality of experiments, which are used to obtain E and A: The traditional Frank-Kamenetskii method requires a number of experiments conducted at several different sizes; The crossing point method may only need one experiment, but has a strict requirement on the positions of thermocouples [11].
In the procedure proposed here, the experiment needed is less resource consuming due to the merit of the model. It would improve the efficiency to evaluate self-heating and ignition hazard for a specific coal. To prove the universal validity of this procedure, further study are required.

Conclusions
In this study, we build a 1D computational model that couples heat transfer, mass transfer, and chemistry to simulate the self-heating ignition of a fuel layer under two different experimental configurations. The model is validated against two high-fidelity experiments (one hot plate experiment and one oven-basket experiment). The good agreement between simulation and experimental data demonstrates that our model can reasonably predict T ig , the time to ignition, and the location of the hot spot (for hot plate experiment). A literature review of the kinetic parameters for coal samples from different origins reveals that there is a compensation effect between the activation energy and exponential factor. This indicates that these two kinetic parameters are correlated and dependent. Combing the compensation effect with the validated model, six different experimental studies covering the two experimental configurations, a wide range of sample sizes (heights ranging from 5 mm to 126 mm), and various coal origins (6 countries) are simulated. The model well predicts T ig for all 24 experiments in the six studies with a predicting error less than 7°C. This is the first time that a computational model is developed to simulate self-heating ignition across experimental configurations (hot plate and oven-basket), scales, and coal origins.
Due to abundant experimental data in the literature for coal, at present our model is validated against coal. The validation shows that the model is capable of simulating two key processes, oxygen diffusion and heat transfer, during self-heating ignition of a reactive porous media. Thus, it is possible to investigate self-heating ignition for different porous fuels based on the present model.