Modified CPD Model for Coal Devolatilization at UCTT Conditions

To study coal pyrolysis behavior at underground coal thermal treatment (UCTT) conditions, a modified CPD (M-CPD) model was developed and evaluated using two scales of experiments as well as two different coals, Utah Sufco and Illinois #6. Compared with the original CPD model, three major aspects were changed, including (1) rewriting the algorithm using MATLAB built-in functions; (2) developing a correlation method between coal elemental composition and the initial functional group population to calculate the effect of cross-linking on the lattice structure; and (3) using a two-activation energy approach for gas generation kinetics modeling. Results predicted by the M-CPD model were compared with those from the CPD model and experiments, and reasons for observed differences between the two models and the experimental data are discussed. The instantaneous devolatilization rates and yields of total volatiles and tars predicted by the M-CPD model were in good agreement with experimental data obtained at UCTT conditions.

Modified CPD Model for Coal Devolatilization at UCTT Conditions Ding Wang,1 Thomas H. Fletcher,2 Swomitra Mohanty,1 Haoquan Hu,3 and Eric G. Eddings1,* 1 Department of Chemical Engineering, the University of Utah, Salt Lake City, Utah, 84112, Unite States 2 Department of Chemical Engineering, Brigham Young University, Provo, Utah 84602, United States 3 State Key Laboratory of Fine Chemicals, Institute of Coal Chemical Engineering, School of Chemical Engineering, Dalian University of Technology, Dalian 116024, China Corresponding Author Email: [e-mail] 1 ABSTRACT To study coal pyrolysis behavior at underground coal thermal treatment (UCTT) conditions, a modified CPD (M-CPD) model was developed and evaluated using two scales of experiments, as well as two different coals, Utah Sufco and Illinois #6. Compared with the original CPD model, three major aspects were changed, including: 1) rewriting the algorithm using Matlab built-in functions; 2) developing a correlation method between coal elemental composition and the initial functional group population to calculate the effect of crosslinking on the lattice structure, and 3) using a two activation energy approach for gas generation kinetics modeling. Results predicted by the M-CPD model were compared with those from the CPD model and experiments, and reasons for observed differences between the two models and the experimental data are discussed. The instantaneous devolatilization rates and yields of total volatiles and tars predicted by the M-CPD model were in good agreement with experimental data obtained at UCTT conditions. 1. INTRODUCTION Cheap natural gas as a result of the shale gas revolution has radically changed the electricity generation market in the U.S., which was historically dominated by coal but is now surpassed by natural gas1. The shift from coal to natural gas is also driven by environmental concerns over climate change caused by high CO2 emissions from coal combustion. The development of clean coal technologies involving low-carbon emissions are necessary from both environmental and economic perspectives. Underground coal thermal treatment (UCTT) is a proposed method for producing light hydrocarbon compounds from underground coal seams (Figure 1). By in-situ heating of coal in the absence of an oxidant, the coal would be thermally decomposed to a product mixture liquids and gas with a higher hydrogen/carbon ratio than the original coal, and a solid 2 product of char, the latter being left underground along with the mineral matter inherent in coal. The UCTT method would reduce the overall CO2 footprint for coal use at the source, since even the residual porous carbon matrix could be utilized for later CO2 sequestration. Other environmental advantages of UCTT include the absence of surface depressions, as the subterranean cavity would be filled with the char matrix, and also the overall reduction in the release of sulfur and heavy metals. A similar in-situ coal utilization process, underground coal gasification (UCG), which produces syngas by controlled coal seam oxidation, has been investigated by many countries at pilot scale2. In contrast with UCG, UCTT heats the coal seam at lower temperatures (e.g., 573-873 K), and operations can be discontinued readily, which improves the overall safety of the process. In addition, UCTT avoids the subsidence issues often encountered with UCG. On the other hand, UCTT requires an external power input to drive the pyrolysis reaction, whereas UCG does not. Figure 1. Example of a UCTT process. The process is not drawn to scale. 3 Our group has been focusing on the UCTT process for several years: Kelly et al.3 evaluated the life-cycle energy and greenhouse gas (GHG) impact of UCTT based on two scales of experiments. The results showed that the feasibility of UCTT is determined by heater temperature, number of heaters and moisture content of the coal seam. Gneshin et al.4 compared particle morphology changes of a Utah bituminous coal at conditions of atmospheric pressure and heating rates of 10 and 0.1 K/min, and found that no plastic deformation occurred at lower heating rates, as compared to raw coal. Zhang et al.5 recently published a review of underground in situ processes relevant to UCTT, and as pointed out in this paper, the models that have been developed for conventional coal pyrolysis may not necessarily be suitable for UCTT due to the vast differences in temperatures and heating rates used. The modeling of a UCTT process is thus still in its infancy. To date, several network models have been developed for conventional pulverized coal devolatilization. Solomon et al. combined the FG (functional group) model with the DVC (depolymerization, vaporization, and crosslinking) model, using a Monte Carlo method to produce the general FG-DVC model6-10. Later, another version of the two sigma FG-DVC model involving percolation theory was developed to speed up the calculation11. Fletcher et al.12-14 made percolation statistics applicable to a Bethe lattice, which is employed to represent the coal molecular structure, allowing a mathematical description of the bridge-breaking process in a closed form. Niksa et al. also adopted the percolation theory of a Bethe lattice in a model called FLASHCHAIN15,16, which when compared to its previous model DISCHAIN17,18, shared the same qualitative features. Besides network model forms, one-step and multi-step kinetic models are increasing in popularity. Ranzi et al.19,20 assumed a devolatilization mechanism in which all coals were considered as a linear combination of three reference coals. The devolatilization of a chosen coal can be expressed as a straightforward weighted combination of the pyrolysis of the reference 4 coals. Among the models mentioned above, the CPD model is widely used to describe the devolatilization processes during combustion and gasification, both of which feature high temperatures and high heating rates. Nevertheless, UCTT is typically a low heating rate and low temperature process, which lies outside of the current CPD model's capability. Therefore, some modifications to the CPD model are necessary for this particular application. This paper will focus on the modification procedure and the analysis of the results obtained from the CPD and the modified CPD (M-CPD) models, as compared to experimental results. UCTT also features conditions such as high pressure and physical confinement that can lead to additional mass transfer considerations and secondary reactions; however, these issues will not be discussed in the current paper. 2. MODIFIED CPD MODEL 2.1. limitations of the current CPD model The CPD model has been proven to predict coal devolatilization behavior (e.g., total volatile and tar evolution) correctly, except for cases of very low heating rate and/or low temperature. Figure 2a illustrates a total volatile yield comparison between TGA data for a Utah Sufco coal and CPD model results at 1 K/min. Like most solid fuel TGA devolatilization data, a typical sigmoid evolution trend can be found in both curves. Below 650 K, a slight increase in apparent yield happens due to moisture and CO2 release (verified by a micro-GC, as shown later), which is possible due to low temperature crosslinking and physical desorption. As the temperature rises, total volatile evolution speeds up. Between 650 and 750 K, the majority of devolatilization reactions happen, and fluidic products such as tar and H2O are generated. Above 750 K, fewer reactions but with high activation energy continue, such as H2 formation due to aromatic 5 condensation. Therefore, the total volatile yield gradually reaches a plateau. In contrast, the CPD model predicts zero yields below 600 K, since the temperature is not high enough for noteworthy numbers of labile bridges to break. As the temperature ramps up, an approximate linear evolution trend is predicted all the way to the onset of the final temperature, and then the slope decreases. Overall, the total volatile yield predicted by the CPD model is close to the TGA data, and the starting point for the non-zero yield is reasonable since products generated at such low temperatures are not of interest. However, the total volatile yields predicted by the CPD model are much lower than the TGA results in the main devolatilization stage (650-750 K). 6 Figure 2. Comparison between the CPD model predictions (red lines) and experiments. (a) Utah Sufco coal TGA result (black lines, daf basis). The experiment was conducted with a heating rate of 1 K/min and a final temperature of 873 K; (b) Utah Sufco coal core pyrolysis results (black dots, daf basis). Experiments were conducted with a heating rate of 10 K/min and a final temperature of 683 K. The comparison of tar evolution (Figure 2b) at 10 K/min indicates an even greater difference. Because of the heterogeneous nature of the core samples, the tar yields exhibit significant data scatter over a large range (0.075-0.18). Even so, it is clear that a significant amount of tar would 7 be generated for these temperature and heating rate conditions. Nevertheless, the tar yield calculated by the CPD model is close to zero (0.003) during the whole experiment. 2.2. Logic flowchart of the modified CPD model To address deficiencies in CPD model predictions at low heating rate or low temperature conditions (UCTT conditions), some modifications are necessary for the CPD model. Figure 3 shows the logic flow charts of the original CPD model (Figure 3a) and the M-CPD model (Figure 3b). As shown in Figure 3a, the CPD model needs input parameters from the interpretation of NMR data of the sample coal, to calculate the structure coefficients such as side chain populations (δ). For cases of the initial fraction of intact bridges (p0) larger than the threshold of a corresponding Bethe lattice, the temperature and rate constants for each time step will be calculated. In this module, a simple lookup table is designed to calculate the inverse of the area under the normal curve. To solve for the new L (labile bridge population), δ, and C (char bridge population) equation set, a hand-written predictor-corrector method is used. In the next module, the total fraction of finite fragments (tar plus metaplastic) and total fraction of gas can be obtained from a Lattice statistics calculation. However, not all components can be released from the mixture due to their different saturated vapor pressures. Therefore, the weight fraction of total fragments with different fragment sizes (for n-mer, ff(n)), and the molecular weight of the one fragment with different size (for nmer, mf(n)), both need to be calculated as well. What's left in the matrix (metaplastic) is consumed by crosslinking reactions with char; however, this crosslinking calculation does not change the structure parameters such as L and C. In the end, flash distillation determines how much volatiles are released from the matrix. The fraction of tar, together with the fraction of gas, combine to 8 provide the total volatiles. In this step, a hand-written secant method is used for the Beta (the fraction of feed that is vaporized) calculation in the Rachford-Rice formulation. Figure 3. Logic flowcharts for the CPD model (a) and the M-CPD model (b). Modules that have been modified are marked in red color. The modifications made in the M-CPD model are shown in Figure 3b, and in general, three main aspects of the CPD model have been changed: 1. Some algorithms were rewritten using Matlab built-in functions (block 1, 2, 3) 9 2. Part of Ranzi's model19,20 was introduced to calculate the crosslinking effect to the lattice structure (block 4,5,6) 3. Two activation energy sets instead of one were applied for computing the gas generation kinetics (block 1) The purpose of rewriting the algorithm was to improve the speed and accuracy of the CPD model. The detailed procedures of the modification will be discussed in the following. The original CPD model's crosslinking module is not related to the change of the bridge population as mentioned earlier. In reality, crosslinking does change the bridge population, which ultimately changes the lattice structure. Two sources of crosslinking are discussed in the current paper; namely, lowtemperature crosslinking (LT crosslinking) and mid-temperature crosslinking (MT crosslinking). For the third modification, a normal distribution of the activation energy for gas release reaction kinetics is always a good assumption in high heating rate pyrolysis conditions. However, for low heating rate conditions, this assumption needs to be reconsidered. All the arguments and the methods of execution will be discussed in the following. 2.3. Rewriting the algorithm using Matlab built-in functions Three hand-written methods in the original algorithm have been changed to Matlab built-in functions including: 1. The Norminv function instead of a lookup table to calculate the inverse of the area under the normal curve (block 1) 2. The Ode45 function instead of a hand-written Predictor-Corrector method to solve the ordinary differential equations (ODE), (block 2) 10 3. The Fzero function instead of a hand-written Secant method to calculate Rachford-Rice formulation in flash distillation (block 3) 11 Figure 4. Variation of the model parameter Beta versus holding time, calculated from the MCPD model (blue lines) and the CPD model (red lines) at different conditions (in the order of heating rate, final temperature): (a) 10 K/min, 673 K. (b) 10 K/min, 1073 K. (c) 60000 K/min, 1073 K. 12 The first two modifications have little effect on the results of L, δ, C, ff(n) and mf(n) when compared to the CPD model, as shown in S1 and S2 (Supplementary information). However, dramatic differences are found in the Beta calculation due to the Fzero function substitution. In the original CPD model, the flash distillation module simulates the separation of the volatiles that will "leave the matrix" and the metaplastic that will stay in the matrix, both of which are calculated in the preceding modules. From a Bethe lattice perspective, the feed stream includes the newly generated finite clusters, together with the remaining finite clusters (the fraction after the crosslinking reaction from the original CPD model, not the LT and MT crosslinking described in this paper). Beta represents the instantaneous gas phase fraction (relatively small clusters) of the feed stream before separation. By solving the Rachford-Rice formulation numerically, the value of Beta and the final tar fraction (the finite clusters that "leave the matrix") can be obtained. Figure 4 shows the comparison for 3 cases; namely, (a) low temperature, low heating rate; (b) high temperature, low heating rate; and (c) high temperature, high heating rate. For case (a), at the early stage below 550 K, the results from both the CPD model and M-CPD model are constant. Although the values are distinct, they are the same as the initial guess (defined by the user), 0.3 and 0, respectively. As the temperature continues to ramp up, some small fragments in the lattice start to leave the matrix (corresponding to real world volatiles vaporization). Note that mass transfer resistances are not considered by either the CPD or the M-CPD model. As the labile bridge breaking eventually initiates, values predicted by the CPD model are close to zero across the whole devolatilization process. On the other hand, a high-low trend is obtained by the M-CPD model, which seems more reasonable. According to experimental observations, gases (small fragments) are the primary products below 650 K, and later liquid (larger fragments) starts to be generated due to labile bridge breaking, which corresponds to the downward trend of Beta as predicted by 13 the M-CPD model, and these values should be greater than zero, since gas generation is always happening at temperatures above 650 K. Figure 4b shows the Beta comparison at a higher final temperature while maintaining the same heating rate as Figure 4a. Besides the discovered difference at a low-mid temperature range (below 900 K), both the CPD and the M-CPD models jump to a value close to 1 after this range, then another distinction appears. A downward trend with fluctuation is observed from the result of the CPD model, while the M-CPD model holds a constant value of close to 1 above 900 K. The latter result is consistent with the fact that, at a relatively high temperature, e.g., >900 K, most labile bridges have already been broken. Thus, less heavy finite clusters can be generated, and the crosslinking reaction from the original CPD model will also continue consuming the "surviving" heavier clusters (metaplastic). With the resulting smaller liquid fraction in the feed stream, and the continued generation of gaseous product, the value of Beta will increase to a value close to 1. To further verify the heating rate and temperature effects on the CPD model, a comparison at high temperature and high heating rate conditions was also conducted. As demonstrated by Figure 4c, besides the existence of the initial gap discussed above, both the CPD and the M-CPD models predict values close to 1 above 900 K, which can reflect the reality as explained above. It can be concluded that the algorithm in the flash vaporization calculation of the M-CPD model works better at low temperature and low heating rate conditions compared to that in the CPD model. 2.4. Modeling the crosslinking effect The crosslinking effect has been studied by many groups. Yang et al.21 built a model for the evolution of the pore structure in a lignite particle during pyrolysis, based on the CPD model. For describing the effects of crosslinking reactions, a correlation between crosslinks and cleaving of side chains was established to modify the number of bridges calculated by the CPD model. 14 However, this model has only been verified for high heating rate and high temperature experiments. Furthermore, many researchers have systematically studied the parameters that control crosslinking, including Solomon et al.22, Suuberg et al.23,24, and Ibarra et al.25,26. Their studies show that two crosslinking events exist: LT crosslinking (below 673 K) occurs prior to tar evolution, with CO2 and H2O evolving simultaneously (in low-rank coals only); MT crosslinking (673-773 K) occurs slightly after tar evolution, hence the tar yield is less affected and it correlates best with methane formation. Their observations indicate that the crosslinking effect in a Bethe lattice system can be calculated from the carboxyl/hydroxyl and methyl group distributions in unreacted coal, since these functional groups are the source of the gases mentioned above, according to the FG model7,8. For different ranks of coal, the initial functional group distribution is vastly different. Zhao et al.27 published a method to correlate these parameters to the elemental composition of the sample coal. By using a set of library coals (with functional group compositions known) to form a triangular mesh in the Van Krevelen diagram, an unknown coal can be interpolated from those of the three library coals based solely on its elemental composition. However, the elemental composition of the Utah Sufco coal is outside of the triangular mesh range. Therefore, a new method needs to be developed. As mentioned in the introduction, Ranzi et al.19 published a correlation that uses three reference coal (molecules) to represent a real coal. If the monomer in the Bethe lattice is represented by these reference molecules, the initial functional group distribution can be calculated based on the functional group numbers on each reference molecule and its weight percentage. Since a Bethe lattice instead of the real coal is considered as the matrix for the modeling, the crosslinking reaction regulations need to be specified as well. 15 2.4.1. Use of Ranzi's model In Ranzi's model19, the first step is coal characterization. The elemental analysis of the coal is normalized to the C, H and O contents on a dry ash-free basis to formulate the sample coal molecule. The second step is to determine the weighting factors of the reference molecules to represent the sample coal. In this step, the high hydrogen content reference molecules COAL1 (C12H11-, 50:50 mol mixture of -C12H10- and -C12H12-), and the high oxygen content reference molecule COAL3 (-C12H12O5-), together with the pure carbon CHARC (-C-) form a triangle in a van Krevelen diagram. The bituminous reference molecule COAL2 (-C14H10O-) locates in the middle of this triangle, so that the sample coal molecule can fall into one of the resulting three zones, and then the weighting factors and the associated initial functional group population can be calculated. The last step, which is the core of Ranzi's model, is to form the complete reaction scheme of the sample coal molecule by using the reference coal molecules reaction table. Since we have adopted the devolatilization mechanism of the CPD model, the last step is not employed. 2.4.2. Modeling the low-temperature crosslinking effect To establish these crosslinking reaction regulations, knowledge relating to the reaction mechanism is necessary. Recent investigations are reviewed below. Hayashi et al.28 found a relationship between the loss of aliphatic carbon and the yield of H2O, and that the CO2 yield did not correlate with the observed heating rate effect on the maximum tar yield. The investigation conducted by Eskay et al.29 suggests that the CO2 evolution from thermal decarboxylation of free carboxylic acids may not be directly related to crosslinking reactions. Mae et al.30 proposed six kinds of crosslinking reactions based on the types of hydrogen bonding formed 16 between the carboxyl and hydroxyl functional groups, and the results show that the number of crosslinks is equal to the amount of H2O formed. To summarize the observations mentioned above, H2O instead of CO2 seems to be a better correlating product for LT crosslinking, which means that decarboxylation may not be necessary for crosslink formation. Therefore, it is highly possible that the crosslink formed is a labile bridge instead of char bridge, due to the difficulty level of the corresponding reactions. In addition, the labile bridge formed will convert to either side chains or a char bridge later, which follows the 0.9 rule from the CPD model relating the side chain rate constant to that of char. After putting all these observations/assumptions into Bethe lattice frame, and taking consideration of the limited variety of the functional groups from the three reference molecules, a model crosslinking process involving only carboxyl and hydroxyl groups can be established. Figure 5. Example of crosslinking formation in a Bethe lattice (coordination number = 3, initial intact bridge population = 0.5). (a) An example configuration of functional groups for LT crosslinking (left circle) and MT crosslinking (right circle). (b) The newly formed labile bridge is represented by the green line, while the newly formed char bridge is represented by a red line. 17 Figure 5 illustrates the idea presented. An example of how LT crosslinks form between clusters is shown in the left yellow circles in Figures 5a and 5b, noting that the OH group inside can be replaced by COOH. This model greatly reduces the complexity of the possible reactions, by considering only carboxyl and hydroxyl groups. Due to the properties of Bethe lattices, no loops can be formed in the matrix, which means that this labile bridge can be formed only when the OH and COOH group are in the same ‘channel' (channel indicates the possible path that exists in a Bethe lattice). Based on these rules, a formula for the initial LT crosslinking population can be established, as shown in Equation 1. w3 w3 w2 w3 CRLT0 = (1 − p0) ∗ (σ+1 ∗ σ+1 + 2 ∗ σ+1 ∗ σ+1) ∗ (σ + 1) = (1 − p0) ∗ w32 +2∗w2∗w3 σ+1 Eq. 1 From left to right, explanations are provided for each term. Firstly, the (1 − p0) refers to the possibility of no bridge existing between two random adjacent clusters (double bridges are not w3 considered); secondly, σ+1 refers to the possibility for a COOH group of a chosen cluster in one w2 channel (one cluster has (σ +1) channels), and σ+1 refers to the possibility for an OH group of a neighbor cluster in one channel. For a chosen cluster, it has (σ + 1) neighbor clusters (which is (σ + 1) possibilities, the last item). Thus, Equation 1 above expresses the coexistence possibility of either a COOH-COOH or COOH-OH pair in one channel (OH-OH pair is not considered since the corresponding reaction happens above 673 K, which is beyond the scope of LT crosslinking30). The coefficient 2 represents the interconvertible nature of a COOH-OH pair between two neighboring sites. The hidden assumption in this formula is: if COOH/OH groups are in one channel, a crosslinking reaction will happen since the hydrogen bond is likely present in between the two groups. 18 The next step is to determine the kinetic parameters. Solomon's work is presented here in Equation 211: E LT k LT = ALT exp (R∗T ) Eq. 2 Since the author assumed that LT crosslinking is a first-order reaction, it is easy to derive the following equations: CRLTnew = CRLTold ∗ exp (- k LT ∗ dt) Eq. 3 Lnew = Lold + k LT ∗ CRLT ∗ dt Eq. 4 A new L for each time step is obtained; therefore, the Bethe lattice structure changes accordingly. 2.4.3. Modeling the mid-temperature crosslinking effect To establish a set of MT crosslinking reaction guidelines similar to that of LT crosslinking, a relevant literature review is necessary. Few papers have investigated the detailed reaction mechanism for MT crosslinking of coal. However, publications in asphaltenes and thermal cracking/coking provide insight to this complicated process. Ranzi et al.31 indicates that the addition reactions between polyaromatic radicals through dealkylation are more favored at a position where methyl side chains are present (larger side chains are avoided by β-scission), and the methyl radicals abstract hydrogen on the substrate to produce methane. Alshareef et al.32,33 found that alkyl-alkyl and alkyl-aryl addition are more favorable reactions than aryl-aryl addition, due to the higher reactivity of alkane side chains. Also, even broader molecular weight distributions than the parent compounds can be obtained in this liquid phase reaction, which is analogous to the parallel reaction mechanism used in the CPD model. 19 After putting all these observations above into the Bethe lattice frame, a model crosslinking process involving methyl groups can be established. An example similar to LT crosslinking is provided in the right yellow circles of Figure 5. In short, one random cluster with a methyl group will form a char bridge with its neighbor cluster when they are not linked in this channel. Because the MT crosslinking happens after the labile bridge breaking stage, the formation of a char bridge rather than a labile bridge seems more reasonable from an energy perspective. The released methyl free radical abstracts hydrogen from the neighboring cluster or its side chain to form a methane molecule. Owing to the free radical's short lifetime and high activity, acquiring hydrogen from the nearest position to stabilize itself is a reasonable assumption (which means they are in the same channel). It is certain that not every methyl group dissociation leads to crosslinking. Therefore, a coefficient is necessary in the calculation, considering the possibility that the cluster is stabilized by other long-lifetime floating free radicals. The last guideline is that the population of newly-formed methyl groups from labile bridge breaking is not taken into account in this calculation, since the chances of char bridge formation are already counted in the original CPD model. Based on the discussion above, the formula for the initial MT crosslinking population can be established, as shown in Equation 5. CRMT0 = (1 − p0) ∗ 1.5∗w1 σ+1 ∗ 1 σ+1 ∗ (σ + 1) ∗ 2 1.9 = (1−p0)∗1.5∗w1∗2 (σ+1)∗1.9 Eq. 5 Once again, the possibility of no bridge existing between two random adjacent clusters is (1 − p0). The term 1.5∗w1 σ+1 refers to the population of CH3 of a chosen cluster in one channel (1.5 comes from the reference molecule COAL119). Since hydrogen exists in almost any cluster/side 20 chain, the population of 1 σ+1 is used; for a chosen cluster, it has (σ + 1) choices to link a neighbor. For the last term, the coefficient 2 represents the interconvertible nature of a CH3-H pair, and the coefficient 1.9 is divided to consider the chances of floating free radical stabilization, which is consistent with the value in the original CPD model. Similar to the LT crosslinking, the kinetic parameters are taken directly from Solomon's work11, and are shown below: E MT k MT = AMT exp ( R∗T ) Eq. 6 CRMTnew = CRMTold ∗ exp (- k MT ∗ dt) Eq. 7 Cnew = Cold + k MT ∗ CRMT ∗ dt Eq. 8 A new C is obtained, therefore the Bethe lattice structure changes accordingly again. 21 2.5. Two activation energy sets for gas generation Figure 6. Light gaseous products evolution over holding time from the Sufco coal core pyrolysis experiment. The experiment was conducted with a heating rate of 1 K/min and a final temperature of 873 K. In addition to the introduction of crosslinking modules, another modification was made to the CPD model to simulate the gas evolution process; namely, the application of two activation energy sets for gas generation. The simplification of using one activation energy for gaseous product prediction in the CPD model is appropriate for high heating rate applications, since all the components are released almost simultaneously during flash devolatilization. With the aid of the standard deviation (normal distribution), the CPD model can predict total volatiles evolution accurately at as low as 10 K/min. However, at lower heating rates, this assumption no longer holds because different gas components are produced in sequential order. An example of evolution from coal at a heating rate of 1 K/min is shown in Figure 6. Note that even the same species may come 22 from different sources (notice the multiple peaks for CO2 and CO), which causes the integrated peak (total gas generation) to become more and more ‘flat'. The reason for this phenomenon is clear, since every individual gas component has its own one or more activation energies, and all the gas generation reactions are independent, according to the FG model7,8. Then why use two activation energies instead of three or more? The idea was born out of Figure 6: most of the peaks in the figure can be grouped into two categories, those associated with the CH4 peak, and those with the H2 peak (notice the two black markers). Solomon et al.8 made a list of activation energies for major gases, in which two categories of 60000 and 80000 cal/mol can be found as well, corresponding to CH4 and H2, respectively. Therefore, to not over-complicate the modeling process, no more than two activation energy sets are used, and the standard deviations for these two gases in the list are adopted in the calculation. In addition, a ratio of the quantity for these two groups is specified as 1:2, and to maintain the accuracy of the high heating rate predictions, the M-CPD model was designed to switch back to one activation energy set when the heating rate is higher than 1 K/second. 3. EXPERIMENTAL SECTION 3.1. Coal sample preparation Powder preparation: Ground Utah Sufco Bituminous coal and Illinois #6 coal were sieved to give particle sizes below 88 μm, and then dried at 383 K for 24 hours to remove the moisture before use. Table 1 provides a standard analysis of the two types of coal. 23 Table 1. Ultimate and Proximate Analyses of the Coal Samples Ultimate analysis (wt % , daf) Proximate analysis (wt %) Coal C H N O Adry Vdaf Fdaf M Sufco 77.85 5.84 1.63 14.06 10.92 45.25 54.75 3.21 Illinois #6 75.6 5.8 1.51 12.8 15.5 44.9 55.1 3.56 Coal core preparation: Coal cores (only for Utah Sufco coal) of approximately 2 cm in diameter were prepared from larger blocks of coal using a 7/8 in. diamond-grit hole saw. Cores were drilled with minimum rotation speed (400 rpm) incrementally while cleaning away dust and debris at intervals to prevent overheating of the core surface. Once a core had reached the depth of the hole saw, a screwdriver was used to pry the core off the block, or it would come off by itself due to a natural cleat having been reached. Cores only above 9 g were accepted for experimentation. 3.2. Two scales of experiments Studies at two experimental scales were conducted including thermogravimetric analysis (TGA) and fixed-bed core pyrolysis experiments. TGA: A Q600 (TA Instruments) TGA was used to obtain the pyrolysis kinetics for the pulverized coal samples (approximate 20 mg for each test). Two heating rates of 1 and 10 K/min were applied, with a nitrogen purge to maintain an inert environment. Sample weight and temperature data were recorded every 1.5 seconds for the entire duration of the experiment. All the TGA data presented in current paper were converted from weight loss to volatiles yields to facilitate comparison with the results calculated from both the CPD and M-CPD models. 24 Fixed-bed core pyrolysis experiment: Details of the equipment used are given in Figure 7. The mass flow controller provided a constant flow (500 ml/min) of nitrogen carrier gas through the apparatus. Coal core samples with a stainless-steel (SS) sleeve were placed horizontally into a 316 SS fixed-bed reactor. The sleeve was used to provide better heat conduction. A thermocouple was touching the surface center on one end of the coal core and monitored the sample temperature. Figure 7. The coal core pyrolysis reactor system, including in-line tar sampling and gas analysis apparatus. Two linear heating rates were applied, 1 and 10 K/min, up to three different final temperatures, 683, 773, and 873 K, respectively. The samples were held at each final temperature for different durations to have an accurate tar devolatilization curve for later comparison with model results. In order to minimize the effect of inhomogeneity of the core samples, replicates of each experiment were carried out. Liquid pyrolysis products were collected by a modified bubbler train system, comparable in design to the Tar Guideline34, which is a tar sampling standard for biomass gasification. The main difference was that two U-shaped SS tubes with glass wool inside replaced 25 the two upstream bubblers in the cooling bath, because of the lower solubility of Utah Sufco coal tar in isopropanol, relative to that of the biomass tar. The temperature of the heating and the cooling bathes were set to 313 K and 258 K, respectively. An activated carbon column was placed downstream to recover the vaporized isopropanol, so that a simple weight difference determination method could be used to calculate the liquid yield (tar + water). Tar yield was obtained after the water concentration was measured by a Karl Fischer titration unit (Mettler Toledo V30). Gas composition was determined by an Agilent 490 Micro-GC equipped with two columns (PoraPlot Q and Molsieve 5 Å). 4. MODELING For the modeling work, the commercial software Matlab was used for programming of the M-CPD model. The CPD model source code was provided by Dr. Thomas H. Fletcher, Brigham Young University. Both the CPD and M-CPD model require coal-specific chemical structure parameters found using 13 C NMR spectroscopy. To reduce the expense of this process, Genetti et al. 35 developed a correlation to determine these NMR parameters, using the proximate and ultimate analysis data for each coal type. In this study, the chemical structure parameters for the two coals used are calculated based on this method, and the results are shown in Table 2 26 Table 2. Chemical Structure Parameters of the Coal Samples Coal p0 c0 Mcluster Mδ (g/mol) (g/mol) σ+1 Sufco 0.483 0.0827 457.8 45.7 4.78 Illinois #6 0.61 0 267 29 4.6 Besides the NMR parameters, both models also require the coal-nonspecific rate parameters for different reactions to be specified, which are listed in Table 3. The additional parameters needed in the M-CPD model are included in the table as well. Most values that are used in both models are kept the same, except for the first parameter Eb. In the M-CPD model, a different value is adopted, which comes from Yan's work36, since a slightly better prediction can be obtained by using the new value. As mentioned above, the M-CPD model uses two activation energy sets for the gas generation calculation, and the parameters associated in this step are Eg1, σg1, Eg2, and σg2. The values for these four parameters are directly cited from the study by Solomon et al.8, and they have also considered the applicability of the parameters for low temperature and low heating rate conditions in their previous study7. Similarly, another of Solomon's papers11 is used for the values of ELT, ALT, EMT, and AMT, which are used in the LT and MT crosslinking modules of the M-CPD model. 27 Table 3. Rate Parameters Used in the CPD and M-CPD Model Parameters Value (CPD) Value (M-CPD) Description Eb 55400 cal/mol 57620 cal/mol Bridge scission activation energy Ab 2.6 * 10 σb 1800 cal/mol 1800 cal/mol Standard deviation for distributed Eb Eg 69000 cal/mol 69000 cal/mol Gas release activation energy Ag 15 15 3 * 10 s s -1 -1 15 2.6 * 10 3 * 10 15 s s -1 -1 Bridge scission frequency factor Gas release frequency factor σg 8100 cal/mol 8100 cal/mol Standard deviation for distributed Eg Eg1 NA 60000 cal/mol New Gas release activation energy σg1 NA 4000 cal/mol Standard deviation for distributed Eg1 Eg2 NA 80000 cal/mol New Gas release activation energy σg2 NA 12000 cal/mol Standard deviation for distributed Eg2 Ecr 65000 cal/mol 65000 cal/mol Crosslinking activation energy Acr 15 3 * 10 s -1 3 * 10 15 s -1 Crosslinking frequency factor ELT NA 45000 cal/mol ALT NA 0.8 * 10 EMT NA 60000 cal/mol AMT NA 0.75 * 10 13 14 s -1 s -1 New LT crosslinking activation energy New LT crosslinking frequency factor New MT crosslinking activation energy New MT crosslinking frequency factor 28 5. RESULTS AND DISCUSSION 5.1. Model modified with Matlab built-in functions Figure 8. Comparison of the M-CPD-1 model predictions (blue lines) with the Sufco coal experimental data (black lines and dots, daf basis), and with the CPD model predictions (red lines). For (a) and (c), a heating rate of 10 K/min and a final temperature of 873 K were used. For (b) and (d), a heating rate of 60000 K/min and a final temperature of 1073 K were used. The CPD model that is only modified with Matlab built-in functions (refers to M-CPD-1) is compared with experiments and the CPD model for two specialized cases: 1) a mid-range temperature with a low heating rate (Figure 8a and Figure 8c) and 2) a high temperature with a high heating rate (Figure 8b and Figure 8d). Due to experimental limitations of the TGA and core 29 pyrolysis reactor, no data could be obtained for the high heating rate case and only the CPD model is used as the standard for this comparison, since it has been previously validated for such conditions of high temperature with high heating rate13,36. As indicated in Figure 8, the total volatiles and tar yields predicted by the M-CPD-1 model are higher than those predicted by the CPD model in both cases. These results indicate that the other two modifications previously discussed, namely, introducing modules to calculate the crosslinking effect to the lattice structure and using two activation energy sets for the gas generation kinetics calculation, may be necessary. Each of these modifications were thus included in the M-CPD model, and while each individual modification provided better predictions, both were required to provide the best agreement with the low heating rate and low temperature experimental data. After all these modifications were made, the predictive capability of the final M-CPD model was validated using experiments that investigated three key factors: evolution kinetics, final temperature and heating rate. 5.2. Prediction of instantaneous devolatilization rates A comparison of the volatiles evolution kinetics for different final temperatures is shown in Figure 9 (total volatiles) for the 1 K/min experiments and Figure 10 (tars) for the 10 K/min experiments. In all three cases of Figure 9, both models underestimate the total volatiles yield before about 650 K, due to the same reason in the above discussion of Figure 2a. From 650 K to the end of each ramping stage (transition point of each temperature line), the trend predicted by the CPD model is essentially a straight line. Although the final yield is close to the TGA data in the 873 K case (Figure 9a), the details of the evolution rate (slope) in the experimental curve are not well represented. 30 On the other hand, the results from the M-CPD model match the experimental data more closely, owing to the more accurate calculation provided by the two activation energy modification. After the final temperature is reached, all the curves tend to level off. Yields from the M-CPD model are quite similar to the experimental data, while the CPD model predicts much lower values at lower temperatures (Figures 9b and 9c), which can be explained mainly by the Beta calculation difference mentioned earlier. 31 Figure 9. Comparison of the M-CPD model predictions (blue lines) with the CPD model predictions (red lines) and Sufco coal TGA results (black lines, daf basis). Experiments were conducted with a heating rate of 1 K/min and different final temperatures of 873 K (a), 773 K (b) and 683 K (c). 32 Figure 10. Comparison of the M-CPD model predictions (blue lines) with the CPD model predictions (red lines) and Sufco coal tar yields data from core pyrolysis experiments (black dots, daf basis). Experiments were conducted with a heating rate of 10 K/min and different final temperatures of 873 K (a), 773 K (b) and 683 K (c). 33 Comparisons for the tar yield are shown in Figure 10. Although the experimental data points show some scatter, the final average yield is about the same (approximately 0.15). For the two relatively higher temperature cases (Figure 10a and b), a clear asymptote can be found in both experimental data and model results, with the same corresponding temperature of 773 K or so. However, no such asymptote can be found in the low temperature case (Figure 10c), as the yield is still increasing after reaching the final temperature of 683 K. Furthermore, there is a gap between the curve of the M-CPD model and the experimental data before the holding time of 10000 s. All these phenomena are related to mass transfer effects, which is not considered in either model. For most softening bituminous coal, the lowest coal/char viscosity shows up around 773 K, at which temperature (and 10 K/min of heating rate) most volatiles can be transported by the bubbles generated instead of pure diffusion. Such a mechanism reduces the effect of mass transfer resistance. In the case of 683 K, mass transfer effects cannot be ignored due to the relatively higher viscosity37,38, even though the swelling ratio of chars measured after the experiments indicates that vast numbers of bubbles are created during the reaction. In general, the M-CPD model can successfully predict the values and tendency of the curve. 5.3. Final temperature effect on the ultimate yield Besides the volatiles evolution kinetics, another useful validation for the M-CPD model is to study the effect of final temperature on the ultimate yields. As shown in Figure 11, the experimental total volatile yields increase with temperature, while the tar yields plateau above 723 K. The reason for this different behavior is that the gaseous products are included in the total volatile yield (in addition to tar), and it is the primary volatile component after most tar has been released. The slope change at 683 K in Figure 11a can be explained by this fact as well, since the rate (slope) of total volatile release decreases without the contribution of the tar release. In terms of the model 34 predictions, the CPD model result develops a more linear trend, which fits well on both ends of the studied temperature range but loses accuracy in between, and this issue is addressed by the MCPD model. Figure 11. Comparison of the M-CPD model predictions (blue lines) with the CPD model predictions (red lines) and the yields data from Sufco coal core pyrolysis experiments (black dots, daf basis). Experiments were conducted with a heating rate of 1 K/min and the samples held at different final temperatures for 2 hours and 40 mins. 35 For the tar yields comparison (Figure 11b), the CPD model result fits well at low temperatures (e.g., 570 K and 605 K), and then gradually deviates from the experimental data at higher final temperatures. In contrast, the M-CPD model result follows the trend of the experimental data, yet overestimates the values below 683 K. Besides the mass transfer effect mentioned above, another possible explanation for this overestimation can be found in Burnham's publication37, in which a correction factor is used to consider the lower volatility of liquid product from oil shale at low heating rate and/or low temperatures conditions. It would be reasonable to extend this observation to coal owing to their similar hydrocarbon nature. For the comparison beyond 683K, the M-CPD model result captures most experimental values, whereas the CPD model does not. The differences in total volatiles and tar yield between the two models can be explained mainly by the modification of the Beta calculation. Since the original CPD model's flash distillation module predicted a smaller Beta parameter at both a relatively lower temperature region (<900 K) and lower heating rate such as 10 K/min (as shown in Figure 4b), the fraction of the finite clusters that separated into the vapor phase was reduced. This reduction indicated a lower final tar yield and also a lower total volatiles yield consisting of the final tar and gaseous products. Although the heating rate of 1 K/min used in these cases is not the same as the one used in the Beta calculation, the observation still holds true. 5.4. Heating rate effect on the ultimate yield The third useful factor that needs to be investigated is the heating rate effect. Due to heating rate limitations of the experimental setup, only the CPD model results are used for comparison. A relatively high final temperature of 1073 K is used as the input in all calculations to ensure the only factor being varied is the heating rate. The unit of time for these comparisons was changed to 36 seconds, since it is more standard for high heating rate applications. The comparison for the total volatile yield is shown in Figure 12a. The CPD model predicts that the yields increase linearly over the heating rate, while the trend predicted by the M-CPD model has more of a ‘hockey stick' shape. As mentioned above, real coal volatiles release through two mechanisms: bubble transport (convection) and diffusion. Higher heating rates facilitate bubble transport, resulting in less residence time in the ramping stage and less secondary reactions of the volatiles, and thus a higher total volatiles yield is expected. For cases with a lower heating rate, although the total volatiles yield contributed by the bubble transport is significantly lower, the entire heating time increases dramatically (e.g., 5 days for the 0.001667 K/second case), which also has positive influence on the total volatile yield. Therefore, it is reasonable to expect the occurrence of a minimum when the heating rate shifts from high to low. 37 Figure 12. Comparison of the M-CPD model predictions (blue dots) with the CPD model predictions (red dots). Input conditions: a final temperature of 1073 K, and the samples held at this final temperatures for 15 seconds. For the tar yield comparison (Figure 12b), both models predict essentially linear trends over the range of heating rates. The CPD model predicts a steeper slope, with the yield of the lowest heating rate close to zero, whereas the corresponding yield from the M-CPD model is about 0.12. As shown in Figure 11, which represents a heating rate of 0.01667 K/second, the M-CPD model predicted 38 an asymptotic tar yield of approximately 0.13 (with a much longer holding time at the final temperature), consistent with the experimental data. Therefore, in the lower heating rate regions of Figure 12b, it is clear that the M-CPD predictions are closer to the experimental value of 0.13 (c.f., second point from the left, 0.01667 K/second). Please note that mass transfer and volatility effects are not considered in these predictions, but may become more significant at increasingly lower heating rates. On the other hand, the values predicted by both models are very close at heating rates higher than 100 K/second, which proves that the M-CPD model preserves the predictive capability of the CPD model at high heating rate conditions. The differences in Beta calculation as discussed above are primarily responsible for the differences in the predictions from two models, as shown in Figure 12. Specifically, the relatively high temperature region needs to be considered as well (c.f. Figure 4b and 4c). The Beta value calculated by the M-CPD model is kept constant (close to 1) above 900 K, while the value from the CPD model decreases with some fluctuations. This distinction leads to both the lower tar yield and total volatiles yield predicted by the CPD model at relatively low heating rate regions, as shown in Figure 12. 5.5. Extension to other coals We also investigated the pyrolysis behavior of Illinois #6 coal to validate the universality of the M-CPD model. Figure 13 shows the total volatiles yield comparisons between the TGA data and the M-CPD model with different final temperatures, 873 K, 773 K, and 683 K, respectively. The input parameters from Table 1 and Table 2 for Illinois #6 coal were used for these calculations, with the same rate coefficients. It clearly demonstrates that the M-CPD model can predict the overall devolatilization trend and the total volatiles yield values for Illinois #6 coal. Therefore, it 39 is anticipated that the M-CPD model would be readily applicable to other types of coal and provide fundamental insights to experimental research. 40 Figure 13. Total volatiles yield comparison of the M-CPD model predictions (blue lines) with the Illinois #6 coal TGA results (black lines, daf basis). Experiments were conducted with a heating rate of 1 K/min and different final temperatures of (a) 873 K, (b) 773 K and (c) 683 K. 41 6. CONCLUSION In summary, a modified CPD model has been developed to predict the coal devolatilization behavior at low temperature and lower heating rate conditions typical of UCTT. The modified CPD model was compared to the original CPD model in terms of their respective ability to predict total volatile yield, tar and instantaneous devolatilization rates, for different heating rates, final temperatures, and coal types. Several concluding comments can be made regarding this investigation: 1. A more reasonable result for the flash evaporation calculation is obtained using the Matlab built-in Fzero function than the existing algorithm in the original CPD model. 2. A correlation method between coal elemental composition and the initial functional group population for LT and MT crosslinking was established by introducing part of Ranzi's model19,20 . 3. A two activation energy kinetic calculation method was developed to model the gas generation behavior at low heating rate conditions. 4. Both total volatiles and tar yields predicted by the M-CPD model agree well with those observed experimentally for Utah Sufco and Illinois #6 coals under UCTT conditions for heating rates of 1 and 10 K/min. 5. Mass transfer and volatility limitations may need to be considered at low temperature (e.g., < 683 K) and/or low heating rate (e.g., <1 K/min) devolatilization conditions. 42 AUTHOR INFORMATION Corresponding Author *E-mail: [e-mail] ORCID Ding Wang: 0000-0001-7520-7768 Thomas H. Fletcher: 0000-0002-9999-4492 Haoquan Hu: 0000-0002-5288-2186 Eric G. Eddings: 0000-0002-0212-9188 Notes Eric Eddings has financial interest in the company Amaron Energy, Inc., which has developed mobile pyrolysis technology for biomass processing. ACKNOWLEDGMENTS This work was partially supported by the United States Department of Energy under Award Numbers DE-FE0001243 and DE-NT0005015. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. 43 REFERENCES (1) U.S. Energy Information Administration. Annual Energy Outlook 2018 with Projections to 2050. 2018. (2) Bhutto, A. W.; Bazmi, A. A.; Zahedi, G. Underground Coal Gasification: From Fundamentals to Applications. Progress in Energy and Combustion Science 2013, 39 (1), 189-214. (3) Kelly, K. E.; Wang, D.; Hradisky, M.; Silcox, G. D.; Smith, P. J.; Eddings, E. G.; Pershing, D. W. Underground Coal Thermal Treatment as a Potential Low-Carbon Energy Source. Fuel Processing Technology 2016, 144, 8-19. (4) Gneshin, K. W.; Krumm, R. L.; Eddings, E. G. Porosity and Structure Evolution during Coal Pyrolysis in Large Particles at Very Slow Heating Rates. Energy and Fuels 2015, 29 (3), 1574- 1589. (5) Zhang, H. R.; Li, S.; Kelly, K. E.; Eddings, E. G. Underground in Situ Coal Thermal Treatment for Synthetic Fuels Production. Progress in Energy and Combustion Science 2017, 62, 1-32. (6) Solomon, P. R.; Hamblen, D. G.; Carangelo, R. M.; Serio, M. A.; Deshpande, G. V. General Model of Coal Devolatilization. Energy & Fuels 1988, 2 (4), 405-422. (7) Solomon, P. R.; Hamblen, D. G. Finding Order in Coal Pyrolysis Kinetics. Progress in Energy and Combustion Science 1983, 9, 323-361. (8) Serio, M. A.; Hamblen, D. G.; Markham, J. R.; Solomon, P. R. Kinetics of Volatile Product Evolution in Coal Pyrolysis: Experiment and Theory. Energy and Fuels 1987, 1 (2), 138-152. (9) Solomon, P. R.; Fletcher, T. H.; Pugmire, R. J. Progress in Coal Pyrolysis. Fuel 1993, 72 (5), 587-597. (10) Solomon P, R.; Hamblen D, G.; Serio M, A.; Yu, Z.; Charpenay, S. A Characterization Method and Model for Predicting Coal Conversion Behaviour. Fuel 1993, 72 (4), 469-488. (11) Solomon, P. R.; Hamblen, D. G.; Yu, Z.; Serio, M. A. Network Models of Coal Thermal Decomposition. Fuel 1990, 69, 754-763. (12) Grant, D. M.; Pugmire, R. J.; Fletcher, T. H.; Kerstein, A. R. Chemical Model of Coal Devolatilization Using Percolation Lattice Statistics. Energy and Fuels 1989, 3 (2), 175-186. (13) Fletcher, T. H.; Kerstein, A. R.; Pugmire, R. J.; Grant, D. M. Chemical Percolation Model for Devolatilization. 2. Temperature and Heating Rate Effects on Product Yields. Energy & Fuels 1990, 4 (1), 54-60. 44 (14) Fletcher, T. H.; Kerstein, A. R.; Pugmire, R. J.; Solum, M. S.; Grant, D. M. Chemical Percolation Model for Devolatilization. 3. Direct Use of 13C NMR Data To Predict Effects of Coal Type. Energy & Fuels 1992, 6 (10), 414-431. (15) Niksa, S.; Kerstein, A. R. FLASHCHAIN Theory for Rapid Coal Devolatilization Kinetics, 1. Formulation. Energy & Fuels 1991, 5 (5), 647-665. (16) Niksa, S. FLASHCHAIN Theory for Rapid Coal Devolatilization Kinetics. 2. Impact of Operating Conditions. 1991, 5 (5), 665-673. (17) Niksa, S.; Kerstein, A. R. The Distributed-Energy Chain Model for Rapid Coal Devolatilization Kinetics. Part I: Formulation. Combustion and Flame 1986, 66 (2), 95-109. (18) Niksa, S. The Distributed-Energy Chain Model for Rapid Coal Devolatilization Kinetics. Part II: Transient Weight Loss Correlations. Combustion and Flame 1986, 66 (2), 111-119. (19) Sommariva, S.; Maffei, T.; Migliavacca, G.; Faravelli, T.; Ranzi, E. A Predictive Multi-Step Kinetic Model of Coal Devolatilization. Fuel 2010, 89 (2), 318-328. (20) Maffei, T.; Frassoldati, A.; Cuoci, A.; Ranzi, E.; Faravelli, T. Predictive One Step Kinetic Model of Coal Pyrolysis for CFD Applications. Proceedings of the Combustion Institute 2013, 34 (2), 2401-2410. (21) Yang, H.; Fletcher, T. H.; Li, S.; Hu, H.; Jin, L.; Li, Y. Model for the Evolution of Pore Structure in a Lignite Particle during Pyrolysis. 2. Influence of Cross-Linking Reactions, Molten Metaplast, and Molten Ash on Particle Surface Area. Energy & Fuels 2017, 31 (8), 8036-8044. (22) Solomon, P. R.; Serio, M. A.; Despande, G. V.; Kroo, E. Cross-Linking Reactions during Coal Conversion. Energy and Fuels 1990, 4 (1), 42-54. (23) Suuberg, E. M.; Lee, D.; Larsen, J. W. Temperature Dependence of Crosslinking Processes in Pyrolysing Coals. Fuel. 1985, pp 1668-1671. (24) Suuberg, E. M.; Unger, P. E.; Larsen, J. W. Relation Between Tar and Extractables Formation and Crosslinking During Coal Pyrolysis. Energy and Fuels 1987, 1 (3), 305-308. (25) Ibarra, J. V.; Moliner, R.; Gavilan, M. P. Functional Group Dependence of Cross-Linking Reactions during Pyrolysis of Coal. Fuel 1991, 70 (3), 408-413. (26) Ibarra, J. V.; Cervero, I.; García, M.; Moliner, R. Influence of Cross-Linking on Tar Formation during Pyrolysis of Low-Rank Coals. Fuel Processing Technology 1990, 24, 19-25. 45 (27) Zhao, Y.; Serio, M. A.; Bassilakis, R.; Solomon, P. R. A Method of Predicting Coal Devolatilization Behaviour Based on the Elemental Composition. Symposium(International) on Combustion 1994, 25 (1), 553-560. (28) Hayashi, J. I.; Takahashi, H.; Doi, S.; Kumagai, H.; Chiba, T.; Yoshida, T.; Tsutsumi, A. Reactions in Brown Coal Pyrolysis Responsible for Heating Rate Effect on Tar Yield. Energy and Fuels 2000, 14 (2), 400-408. (29) Eskay, T. P.; Britt, P. F.; Buchanan, A. C. Does Decarboxylation Lead to Cross-Linking in Low-Rank Coals? Energy and Fuels 1996, 10 (6), 1257-1261. (30) Mae, K.; Maki, T.; Miura, K. A New Method for Estimating the Cross-Linking Reaction during the Pyrolysis of Brown Coal. Journal of Chemical Engineering of Japan 2002, 35 (8), 778-785. (31) Ranzi, E.; Dente, M.; Goldaniga, A.; Bozzano, G.; Faravelli, T. Lumping Procedures in Detailed Kinetic Modeling of Gasification, Pyrolysis, Partial Oxidation and Combustion of Hydrocarbon Mixtures. Progress in Energy and Combustion Science 2001, 27 (1), 99-139. (32) Alshareef, A. H.; Scherer, A.; Tan, X.; Azyat, K.; Stryker, J. M.; Tykwinski, R. R.; Gray, M. R. Formation of Archipelago Structures during Thermal Cracking Implicates a Chemical Mechanism for the Formation of Petroleum Asphaltenes. Energy and Fuels 2011, 25 (5), 2130- 2136. (33) Alshareef, A. H.; Scherer, A.; Tan, X.; Azyat, K.; Stryker, J. M.; Tykwinski, R. R.; Gray, M. R. Effect of Chemical Structure on the Cracking and Coking of Archipelago Model Compounds Representative of Asphaltenes. Energy and Fuels 2012, 26 (3), 1828-1843. (34) Good, J.; Ventress, L.; Knoef, H.; Zielke, U.; Lyck Hansen, P. Sampling and Analysis of Tar and Particles in Biomass Producer Gases - Technical Report. 2005, 1-44. (35) Genetti, D.; Fletcher, T. H.; Pugmire, R. J. Development and Application of a Correlation of 13 C NMR Chemical Structural Analyses of Coal Based on Elemental Composition and Volatile Matter Content. Energy & Fuels 1999, 13 (1), 60-68. (36) Yan, B.; Cheng, Y.; Xu, P.; Cao, C.; Cheng, Y. Generalized Model of Heat Transfer and Volatiles Evolution Inside Particles for Coal Devolatilization. AIChE Journal 2014, 60 (8), 2893-2906. (37) Oh, M. S.; Peters, W. A.; Howard, J. B. An Experimental and Modeling Study of Softening Coal Pyrolysis. AIChE Journal 1989, 35 (5), 775-792. 46 (38) Yu, J.; Lucas, J. A.; Wall, T. F. Formation of the Structure of Chars during Devolatilization of Pulverized Coal and Its Thermoproperties: A Review. Progress in Energy and Combustion Science 2007, 33 (2), 135-170. 47 Supplementary information Modified CPD Model for Coal Devolatilization at UCTT Conditions Ding Wang,1 Thomas H. Fletcher,2 Swomitra Mohanty,1 Haoquan Hu,3 and Eric G. Eddings1,* 1Department of Chemical Engineering, the University of Utah, Salt Lake City, Utah, 84112, Unite States 2Department of Chemical Engineering, Brigham Young University, Provo, Utah 84602, United States 3State Key Laboratory of Fine Chemicals, Institute of Coal Chemical Engineering, School of Chemical Engineering, Dalian University of Technology, Dalian 116024, China Figure S1. Bridge dynamic population parameters on a per-site basis as a function of holding time. (a) Labile bridges; (b) side chains; (c) char bridges. Calculations were conducted with a heating rate of 10K/min and a final temperature of 673K. Figure S1 shows the comparison results of the three kinds of bridge variations with time from both the CPD and M-CPD models as the pyrolysis proceeds. As shown in Figure a, the labile bridge (L) population maintains the initial value until 650 K then decreases rapidly, since this temperature is high enough to trigger some low activation energy reactions. The relatively small fraction of these reactants (labile bridges) quickly consume (break) to form either side chains or char bridges, resulting in a plateau when the final temperature is held for 2 hours. On the other hand, the population of the non-zero side chain (Figure b) or char bridge (Figure c) starts to increase at about 650 K, with a plateau reached due to the same reason as the labile bridge. Hence, the timing for the changes in bridge population (up and down) are perfectly matched. For comparison of the two algorithms, as indicated clearly by the three graphs, the explicit Runge-Kutta-based ODE45 algorithm and the Norminv Matlab built-in functions obtained almost the same results versus the hand-written Predictor-Corrector algorithm and the look-up table. Even though the lookup table greatly simplified the calculation process of the inverse of the area under the normal curve, the accuracy is maintained. Therefore, no significant differences for the parameter of the bridge population were found for this comparison. Figure S2. Dynamic weight fraction of total mass fragment of size n=1 (a), n=4 (b); dynamic molecular weight of one fragment of size n=1 (c), n=4 (d). 1 and 4 are randomly picking numbers. As demonstrated by Figure a and b, more fragments were generated as the labile bridges broke, leading to an increasing trend for the ff(n). For different n, the magnitude is different, since the possibility for smaller fragment formation is much higher. In contrast, the mf(n) evolution curves (Figure c and d) drop down at the corresponding temperature, since the increased mass of side chains (δ) does not compensate the mass loss of labile bridges (L) converting to char bridges (C). In addition, a bigger mer number (e.g., 1 to 4) causes a magnitude difference as well. As for the comparison of the two algorithms, no unusual distinction can be found. The two computation sample curves coincide almost completely.