Introduction

Cellulose pyrolysis has been a topic of scientific research since the early twentieth century (Lédé 2012). Much of the present-day interest stems from pyrolysis-based technologies that enable the conversion of lignocellulosic biomass into transportation fuels and commodity chemicals (Huber et al. 2006). Another reason for the interest is the fact that cellulose pyrolysis plays a central role in the combustion of wood, one of the oldest and most widely used building materials in the world. Despite decades of experimental research efforts and abundance of data, many (if not most) aspects of the decomposition process still remain debated or unknown (Lin et al. 2009; White et al. 2011; Lédé 2012; Mettler et al. 2012b; Burnham et al. 2015). This may be because sample preparation, sample size, temperature range, heating rate, gas/liquid environment, and experimental set-up all have been shown to have an effect on the process and, consequently, the interpretation of the results.

Universally it is agreed that cellulose does indeed chemically decompose due to an elevated temperature, that one of the major volatile decomposition products is levoglucosan (LGA), and that the process usually leaves a solid char residue whose chemical and physical composition is mostly unknown. A spectrum of lumped reaction kinetics models exist for describing this process, the Broido-Shafizadeh model (Bradbury et al. 1979) being one of the most widely used. It assumes that cellulose first converts to “active cellulose” and then further degrades via parallel pathways either to volatile organic compounds or to char and gases (CO, CO2, H2O). The lumped models rely on experimental parameters obtained from thermogravimetric analysis (TGA) or related techniques, usually being limited to global mass loss kinetics and yields of final gaseous products.

Describing the condensed-phase pyrolysis chemistry, i.e. the elementary reaction mechanisms underlying pyrolysis, is regarded as one of the fundamental challenges for a thorough understanding of cellulose (and in general polymer) pyrolysis (Mettler et al. 2012b). In the recent years, there have been attempts to complement the experimental data on cellulose pyrolysis through electronic structure calculations and atomistic simulations that are capable, one way or the other, of capturing individual chemical reactions that occur during the decomposition process (Liu et al. 2011b; Zhang et al. 2011b; Mayes and Broadbelt 2012; Agarwal et al. 2012; Hosoya and Sakaki 2013; Zhang et al. 2014; Murillo et al. 2015; Zhang et al. 2015a, b; Zheng et al. 2016). This is a formidable challenge, as the pyrolysis chemistry of cellulose is highly complex, involving thousands of chemical reactions and yielding hundreds of volatile species (Zhou et al. 2014). Clearly, any such computational approach must be carefully focused to address a specific detail.

While electronic structure calculations based on density-functional theory (DFT) can yield highly accurate information on pre-determined reaction paths, what is really desired is a technique capable of predicting these paths. Predictive molecular-level modelling of pyrolysis has become accessible through the development of reactive molecular dynamics (RMD) methods. Ab initio MD (AIMD) methods, which have been around longer, have the inherent capability of describing bond breaking and forming, along with other phenomena relevant for pyrolysis chemistry. However, AIMD is computationally expensive compared to the RMD approach, and thus more restricted when it comes to the size and time scales of the simulations. There are specialised methods for accelerating AIMD, but each has its own issues and subtleties (Voter 1997a, b, 1998; Sørensen and Voter 2000; Laio and Parrinello 2002). The RMD methods employ reactive interatomic potentials (Abell 1985; Tersoff 1988; Brenner 1990; Stuart et al. 2000; van Duin et al. 2001; Brenner et al. 2002; Su and Goddard 2007; Shan et al. 2010; Knippenberg et al. 2012; Liang et al. 2013), or are based on the use of switching functions (Nyden et al. 1992). Both approaches have been used for studying the thermal decomposition of polymers. ReaxFF (van Duin et al. 2001), a widely used reactive interatomic potential, has been utilised i.a. in RMD simulations of the thermal decomposition of polydimethylsiloxane (Chenoweth et al. 2005), 1-methylnaphthalene (Leininger et al. 2008), and polyethylene (Liu et al. 2014). Nyden et al. have developed a switching-function based RMD software that has been used to study the thermal decomposition of polyethylene, polypropylene, polyisobutylene (Nyden et al. 2004) and poly(methyl methacrylate) (Stoliarov et al. 2003).

Concerning cellulose, Agarwal et al. utilised metadynamics-accelerated AIMD to study the decomposition pathways of cellulose pyrolysis at 327 °C and 600 °C (Agarwal et al. 2012). Their simulations revealed a variety of possible initial reaction steps. Based on estimates of the reaction free-energy barriers, they concluded that a ring contraction mechanism was favourable at low temperatures, leading to a precursor for hydroxy-methylfurfural formation (via depolymerisation). At high temperatures, the most favourable pathway was glycosidic bond cleavage and subsequent formation of an epoxide precursor to LGA. Recently, Zheng et al. (2016) used ReaxFF-based RMD simulations to study the initial reaction mechanisms of cellulose pyrolysis, and the effect of temperature on the product distribution. Their model system consisted of an amorphous arrangement of six cellulose chains, each having 60 β-D-glucopyranose units. They considered temperatures ranging from 500 to 1400 K for the overall product distribution, and 1500–2800 K for the production of inorganic gases (H2O, CO, CO2). They proposed a detailed decomposition reaction scheme based on the analysis of the MD trajectories. Two major decomposition products were identified as glycolaldehyde and LGA (or its precursor). We note also that various non-reactive MD studies on the high-temperature behaviour of cellulose exist (Agarwal et al. 2011; Matthews et al. 2011; Bergenstråhle et al. 2007; Chen et al. 2012; Zhang et al. 2011a, b; Paavilainen et al. 2011).

In the present work, we used stochastic RMD simulations to identify and analyse the primary thermal decomposition reactions of an isolated cellulose molecule. We considered a range of models to capture the possible effects of molecular weight and initial conformation, and a range of temperatures to determine the kinetics of the dominant reactions. We adopted the approach from a study on the chain fragmentation of polyethylene (Smith et al. 2011). However, we used the ReaxFF reactive force field instead of the switching-function based approach, which was made possible by the availability of a ReaxFF parameterisation for general hydrocarbon and carbohydrate compounds (Chenoweth et al. 2008). The present approach corresponds to gas-phase conditions, i.e. a situation without the chemical and electrostatic environment of crystalline (or amorphous) cellulose. Solid-phase conditions are expected to affect the pyrolysis reactions (Mayes and Broadbelt 2012; Hosoya and Sakaki 2013), making the gas-phase approach more representative of sparse amorphous cellulose. Moreover, focusing on isolated molecules contributes to a bottom-up approach, in which we build towards solid-phase pyrolysis simulations in a controlled manner.

Methods

The outline of the computational analysis is the following. Repeated RMD simulations of an isolated cellulose molecule were performed at a range of temperatures. Each simulation was continued until the first decomposition reaction was observed. The time of the initial reaction was determined for each simulation, yielding a distribution of decomposition times for each temperature. The temperature-dependency of the average reaction rate constant was then calculated based on these distributions, and the corresponding kinetic parameters were determined by fitting against the Arrhenius equation. Details of the molecular models, the molecular dynamics set-up, the stochastic simulations and the post-processing techniques are given in what follows.

Molecular models

Cellulose molecules with degrees of polymerisation (DP) 8, 16, 32 and 64, both in elongated and coiled conformations, were used to consider the possible effects of molecular weight and initial conformation on the decomposition process. The cellulose models were based on crystallographic unit cell data for the cellulose Iβ polymorph (Nishiyama et al. 2002). The considered DPs are much lower than those encountered in natural cellulose. This choice was dictated by the requirement for feasible simulation times. Structural optimisation was carried out to adapt the molecules to the used interaction model. The optimisation was performed using a combination of the steepest descent, conjugate gradient and Hessian-free truncated Newton algorithms.

The coiled conformations were created as follows. The elongated cellulose molecule was simulated in the canonical ensemble (i.e. constant NVT) at a temperature of 300 K for a duration of 100 ps. During this time, the molecule would assume a compact globular conformation. A snapshot of the molecule was taken at 10 ps and at the end of the simulation, corresponding to a slightly coiled and a fully coiled conformation, respectively. The procedure was repeated five times, each time with a different initial velocity vector, yielding five variations of both cases. As an example, Fig. 1 shows an elongated, a slightly coiled and a fully coiled DP-32 cellulose molecule. It is noteworthy, that the used interatomic potential does not predict a linear configuration for an isolated cellulose molecule.

Fig. 1
figure 1

An elongated, a slightly coiled and a fully coiled DP-32 cellulose molecule. Colour code: (black) carbon, (grey) hydrogen, (red) oxygen

Full size image

Molecular dynamics set-up

All simulations were performed using the classical MD code LAMMPS (Plimpton 1995) (SVN revision 8489). The cellulose molecule was embedded in a cubic simulation domain with periodic boundary conditions. The side length of the domain was chosen to be slightly greater than the length of the elongated conformation to avoid self-interactions over the domain boundaries.

ReaxFF, a bond-order based reactive interatomic potential (van Duin et al. 2001), was used to describe the interatomic bonding. The parameterisation by Chenoweth et al. (2008), originally developed for the high-temperature gas-phase oxidation of hydrocarbons, was adopted for this work. It has been used in various studies for simulating the thermal decomposition of hydrocarbon and carbohydrate compounds (Chenoweth et al. 2009; Jiang et al. 2009; Desai et al. 2011; Salmon et al. 2009a, b; Wang et al. 2011; Liu et al. 2011a; Cheng et al. 2012). The uncertainties related to the ReaxFF approach in general, and the parameterisation by Chenoweth et al. in particular, are acknowledged. ReaxFF provides, by design, an approximate solution to the quantum mechanical problem at hand, and the propagation of uncertainty from the parameterisation to the predicted chemistry is non-trivial. It should be noted that ReaxFF treats electronic interactions in an implicit manner. Thus, it does not enable predictions with regard to electron transfer (cf. homolytic and heterolytic fission).

The Fortran implementation of ReaxFF integrated into LAMMPS was used. There are four user-given parameters that affect its behaviour: the hydrogen-bond cut-off, the hydrogen-bond function style switch, the triple bond stabilisation switch and the precision of charge equilibration. A hydrogen-bond cut-off of 6 Å was used. The new hydrogen bond function style was chosen over the old one. Stabilisation of all triple bonds was set off according to the instructions given by the authors of the ReaxFF parametrisation. The precision of charge equilibration was set to 10−6.

The equations of motion were solved using velocity-Verlet integration with a constant time step of 0.1 fs. The time step was chosen based on preliminary simulations and considerations presented in similar studies (Stoliarov et al. 2003; Smith et al. 2011; Knyazev 2007). This value is also recommended for high-temperature (~2500 K) simulations by the authors of the used ReaxFF parameterisation. Due to the nature of the simulations, neighbour list update was performed when at least one of the atoms had moved >1 Å since the previous update.

At the beginning of each simulation, the atoms were given random velocities that would produce the desired temperature in accordance with the classical equipartition theorem. The total linear and angular momentum of the system were adjusted to zero. The Berendsen thermostat (Berendsen et al. 1984) was used for temperature control, which results in an NVE-scaling scheme (an NVT-NVE hybrid). The NVE-scaling scheme with Berendsen thermostatting was chosen over true NVT Nosé–Hoover type approach (Tobias et al. 1993; Shinoda et al. 2004). This choice was based on preliminary simulations, where considerable temperature fluctuations were observed when using the Nosé–Hoover type approach. When using the Berendsen thermostat for temperature control, the fluctuations were an order of magnitude smaller. The same conclusion was reached by Stoliarov et al. in their study on the thermal decomposition of PMMA (Stoliarov et al. 2003). The damping parameter of the thermostat was set to 10 fs (i.e. 100 time steps), as suggested in the LAMMPS documentation.

Stochastic simulations

Repeated decomposition simulations were carried out at the temperatures 1400, 1600, 1800, 2000 and 2200 K. The simulations were realised in an iterative manner. The molecular dynamics were simulated for a predetermined time, after which the bond information output was parsed to find out if the cellulose molecule was still intact. If yes, the simulation would continue for another iteration. If not, the last time step with an intact molecule was determined, and the following time step was considered the decomposition time. The use of artificially high temperatures was dictated by the need for accelerated dynamics—a prerequisite for collecting a sufficient number of decomposition events with reasonable computational cost. It is noteworthy, that the use of high temperatures might open kinetic pathways that are not dominant at lower temperatures. This could potentially impair the correspondence between predicted reaction pathways and those encountered in practical pyrolysis conditions. One hundred decomposition simulations were performed for each combination of molecular weight (four), initial conformation (eleven) and temperature (five). There were, however, two exceptions: (1) for the DP-8 molecules, the lowest temperature was omitted, and (2) for the DP-64 molecules, variations of the coiled conformations were omitted. The omissions were made to reduce the total computational cost of the analysis. With these exceptions, the total number of simulations was 16,900.

Post-processing

Information on chemical bonding was output at 100 time step intervals for temperatures above 1400 K, and at 200 time step intervals at 1400 K. Thus, the time resolution for tracking chemical reactions was 10 fs in the former case and 20 fs in the latter. Lower time resolution was used for the lowest temperature to reduce the size of the bond information file. A bond-order cut-off value of 0.1 was used to distinguish between an intact and a broken bond.

In-house post-processing tools were used to parse the bond information files and to determine for each recorded time step of each simulation the following:

  • the number of molecules in the system,

  • the chemical composition of each molecule,

  • the number of pyranose, furanose and other ring structures in each molecule, and

  • the state of the β(1 → 4)-glycosidic bonds of the original cellulose molecule.

Taken together, this information was used to deduce whether the decomposition event was (1) a depolymerisation reaction (i.e. chain scission via glycosidic bond cleavage), (2) some other chain scission reaction, or (3) a release of a low molecular weight product (LMWP); and whether ring opening or contraction reactions were involved. For the purpose of this study, decomposition products with less than six carbon atoms (cf. β-D-glucopyranose) were considered LMWPs. For decomposition events involving glycosidic bond cleavage, it was further determined whether the glycosidic oxygen would remain bonded towards the reducing or the non-reducing end of the cellulose molecule. The smallest set of smallest rings (SSSR) algorithm, as implemented in the OpenBabel API (O’Boyle et al. 2011), was used for ring perception.

Specific algorithms were developed to search for LGA within the decomposition products, considering both individual molecules and LGA-end groups. A similar search was carried out for glycolaldehyde, formaldehyde, formic acid and the hydroxymethyl radical. Some smaller products such as water, carbon monoxide, carbon dioxide and the hydroxyl radical could be identified by their chemical composition alone.

Average reaction rate constants were determined for the depolymerisation reactions, considering separately each combination of molecular weight, initial conformation and temperature (Smith et al. 2011):

$$k = \frac{N}{{N_{\text{b}} }}\left( {\mathop \sum \limits_{{{\text{i}} = 1}}^{N} t_{\text{i}} } \right)^{ - 1}$$
(1)

where N is the total number of simulations, \(N_{\text{b}}\) is the number of β(1 → 4)-glycosidic bonds in the cellulose molecule and \(t_{\text{i}}\) is the decomposition time observed in the ith simulation. The Arrhenius equation relates the reaction rate constant with temperature through the frequency factor A and the activation energy \(E_{\text{a}}\). These were determined using linear least-squares fitting against the rate constant data.

Results and discussion

General observations

We performed altogether 16,900 simulations, recording the first decomposition event from each. Majority of the reactions (97.89%) resulted in two product molecules; with three products observed in 2.05%, and four products in only 0.06% of the simulations. The decomposition events could be divided into three main categories: (1) depolymerisation reactions, (2) other chain scission reactions, and (3) release of LMWPs. The initial decomposition event was identified as a chain scission reaction in 95.54% of the simulations. Of this figure, 90.75% was due to depolymerisation reactions (an example shown in Fig. 2), while other chain scission reactions contributed the remaining 4.79%. The initial decomposition event was identified as a release of an LMWP in 3.99% of the simulations. Finally, 0.47% of the events could not be identified by our algorithm. The proportions remained roughly the same when considering separately the different molecular weights, initial conformations and temperatures. Only the release of LMWPs showed an effect, becoming slightly more frequent towards the higher temperatures. The decomposition times were in accordance with the exponential distribution, as would be expected for a homogeneous Poisson process.

Fig. 2
figure 2

Example of an observed depolymerisation event: a the reducing end of an intact cellulose molecule, b the terminal β(1 → 4)-glycosidic bond breaks and c a hydroxymethyl radical is released in the process. The glycosidic oxygen remains bonded towards the non-reducing end of the cellulose molecule. Colour code: (dark grey) carbon, (light grey) hydrogen, (red) oxygen

Full size image

Depolymerisation reactions

Majority of the decomposition events were identified as depolymerisation reactions. This is in line with previous experimental and theoretical evidence suggesting that cellulose pyrolysis begins with depolymerisation reactions that ultimately lead to the oligo- and monosaccharide level (Pouwels et al. 1989; Várhegyi et al. 1994; Mamleev et al. 2006; Lin et al. 2009; Lédé 2012; Agarwal et al. 2012; Mayes and Broadbelt 2012; Hosoya and Sakaki 2013; Zhou et al. 2014; Zheng et al. 2016).

We determined the frequencies of end-chain and mid-chain depolymerisation reactions to find out if either of them was dominant. We found that the frequencies are in accordance with a probabilistic model that assumes uniform scission probability throughout the chain (Fig. 3). A further bond-wise analysis confirmed that each β(1 → 4)-glycosidic bond was an equally probable site for cleavage, implying that unzipping or other selective depolymerisation would not be favourable in gas-phase conditions. We also examined whether the glycosidic oxygen would prefer to remain bonded towards either end of the cellulose molecule. We found that bonding towards the non-reducing end was roughly twice as frequent as bonding towards the reducing end. Furthermore, we found that roughly half of the depolymerisation reactions involved a prior or concurrent opening of a pyranose ring. Based on the ring perception data, we could estimate that ring contraction was involved in 5–10% of the reactions.

Fig. 3
figure 3

Ratio of the frequency of end-chain depolymerisation events to the corresponding expectation value in the case of uniform scission probability of the β(1 → 4)-glycosidic bonds. The grey bars indicate the reference ratio of one to one. The error bars indicate standard deviation

Full size image

We determined average reaction rate constants for the depolymerisation reactions, considering separately each combination of molecular weight, initial conformation and temperature. The resulting Arrhenius plots are shown in Figs. 4, 5 and 6, and the corresponding kinetic parameters are given in Tables 1 and 2. The rate constants fall in straight lines in the ln (k) versus T −1 plots, indicating that the dominant reaction mechanism stays unchanged throughout the considered temperature range. The only exception is the case of the elongated DP-64 cellulose molecule, for which there is an obvious deviation from linearity (see Fig. 4). We were not able to attribute this effect to any evident change in the decomposition process. Due to this anomaly, we chose to omit the elongated conformation in the remaining discussion. Neither molecular weight nor initial conformation showed any significant effect on the kinetic parameters, as shown in Figs. 7 and 8. Moreover, the reactions in which the glycosidic oxygen remains bonded towards the non-reducing end, and the ones in which it remains bonded towards the reducing end, did not show any difference within the limits of uncertainty. Finally, taking an average over the kinetic parameters of the coiled conformations yields \(\left( {1.07 \pm 0.12} \right) \cdot 10^{15}\) s−1 for the frequency factor, and (171 ± 2) kJ mol−1 for the activation energy.

Fig. 4
figure 4

Arrhenius plots for depolymerisation reactions in the elongated cellulose molecules. Linear least-squares fits are shown as dashed lines

Full size image
Fig. 5
figure 5

Arrhenius plots for depolymerisation reactions in the 10 ps coiled cellulose molecules. Each data point is an average over the five variant conformations. Linear least-squares fits are shown as dashed lines

Full size image
Fig. 6
figure 6

Arrhenius plots for depolymerisation reactions in the 100 ps coiled cellulose molecules. Each data point is an average over the five variant conformations. Linear least-squares fits are shown as dashed lines

Full size image
Table 1 Arrhenius kinetic parameters for the depolymerisation reactions: DPs 8 and 16. The uncertainties are propagated from the standard deviations of the linear fitting parameters
Full size table
Table 2 Arrhenius kinetic parameters for the depolymerisation reactions: DPs 32 and 64. The uncertainties are propagated from the standard deviations of the linear fitting parameters
Full size table
Fig. 7
figure 7

Frequency factors for the depolymerisation reactions across the considered degrees of polymerisation. The error bars indicate uncertainty propagated from the standard deviations of the linear fitting parameters

Full size image
Fig. 8
figure 8

Activation energies for the depolymerisation reactions across the considered degrees of polymerisation. The error bars indicate uncertainty propagated from the standard deviations of the linear fitting parameters

Full size image

It is noteworthy, that the cellulose molecules undergo conformational change during the decomposition simulations. We observed a wide range of decomposition times in each set of simulations, which is characteristic of the exponential distribution. The average decomposition times ranged from 0.2 ps to 140 ps, depending on the temperature (longer times towards lower temperatures). At room temperature, the conformational change from elongated to fully coiled occurs in <100 ps for all considered molecular weights (see Methods chapter, “Molecular models” section). The behaviour is similar at the elevated temperatures. Thus, at the decomposition time, the initially elongated or slightly coiled molecules are somewhere between the initial and the fully coiled conformation.

We now compare our results against kinetic parameters extracted from cellulose pyrolysis experiments. These are typically based on simplified kinetic models that are parameterised against global mass loss measurements. Thus obtained kinetic parameters represent an average over the underlying elementary reaction mechanisms, and may also include contributions from, i.a., the vaporisation of volatile products. The comparison is, however, justifiable, since depolymerisation reactions have been associated with the rate-limiting step in various models (Bradbury et al. 1979; Várhegyi et al. 1994; Lédé et al. 2002; Mamleev et al. 2006). Lin et al. (2009) give a summary of kinetic studies of cellulose pyrolysis involving various kinetic models, sample materials and reaction conditions. They report activation energies ranging from 48 kJ mol−1 to 282 kJ mol−1 and frequency factors ranging from of 180 s−1 to 1.33 × 1023 s−1. Agarwal et al. (2012) point out that a much narrower range of activation energies, 190–200 kJ mol−1, corresponds to physically plausible frequency factors (1013–1014 s−1). This subset of apparent activation energies is remarkably close to our results. It is noteworthy, that considering the solid-phase electrostatic environment might lead to somewhat higher predictions for the activation energy, as demonstrated by Hosoya and Sakaki (2013).

We may also compare our results against kinetic parameters extracted from DFT-based electronic structure calculations. Agarwal et al. (2012) used metadynamics-accelerated AIMD simulations (at 600 and 873 K) to study nascent pyrolysis reactions in crystalline cellulose. They observed various reactions leading to depolymerisation, for which they reported free energy barriers ranging from 130 kJ mol−1 to 527 kJ mol−1. Roughly half of these fall within the narrower range from 151 kJ mol−1 to 180 kJ mol−1, which also encompasses our results for the activation energy (166–176 kJ mol−1). Pyranose ring opening, contraction or fragmentation was involved in roughly half of the reactions. Altogether, most of the proposed reaction mechanisms are compatible with our observations.

Hosoya and Sakaki (2013) used DFT calculations to study the depolymerisation of cellulose via a concerted mechanism for glycosidic bond cleavage and LGA formation. They determined activation energies for the reaction using one-, two- and three-chain models of glucose hexamers, thus considering both gas-phase and solid-phase electrostatic environments. For the gas-phase conditions, they reported values ranging from 160 kJ mol−1 to 198 kJ mol−1, depending on the site of cleavage. The range of activation energies is compatible with our results. The concerted mechanism was also studied by Mayes and Broadbelt (2012), who reported an activation energy of 231 kJ mol−1 for the reaction in the gas-phase.

Remarkably, we did not observe the formation of LGA in our simulations, neither as individual molecules nor as LGA-end structures, even though our activation energy for the depolymerisation reaction agrees with both experiments and DFT calculations. Due to the stochastic nature of our study, we could not record and analyse the trajectories with sufficient accuracy to reveal the details of the reaction path and the nature of the transition state. Therefore, we cannot provide a detailed account of the mechanism. We could however observe that the glycosidic oxygen preferred to remain bonded towards the non-reducing end of the cellulose molecule, while the concerted mechanism for LGA formation suggests the opposite.

One reason for the lack of LGA formation in our simulations may be the high temperatures, which could lead to radical chemistry dominating over the concerted routes to LGA. Another factor may be the absence of the chemical and electrostatic environment of crystalline (or amorphous) cellulose, especially since the presence of hydroxyl groups on neighbouring chains has been proposed to catalyse concerted reactions (Assary and Curtiss 2012; Seshadri and Westmoreland 2012).

Release of low molecular weight products

The initial decomposition event was identified as a release of an LMWP in 3.99% of the simulations. The species most commonly observed were glycolaldehyde (32.8%), water (17.5%), the hydroxymethyl radical (9.6%), the hydroxyl radical (6.1%), formaldehyde (6.1%) and formic acid (5.5%). These account for more than three quarters of all the observed LMWPs. The formation of glycolaldehyde and formic acid was always accompanied by an opening of a pyranose ring. Ring opening reactions were also involved, with a few exceptions, in the release of formaldehyde. Dehydration reactions, on the other hand, occurred mostly without prior or concurrent ring opening. The observed LMWPs were formed directly from the cellulose molecule, and no reaction intermediates were involved. Based on the ring perception data, we could estimate that ring contraction was involved in 2–9% of the reactions. The stable species, namely glycolaldehyde, water, formaldehyde and formic acid, are recognised products of cellulose pyrolysis (Vinu and Broadbelt 2012; Mettler et al. 2012a, b; Zhou et al. 2014).

The ReaxFF-MD study by Zheng et al. (2016) provides a relevant point of comparison to our results. Firstly, they report glycolaldehyde and water among the major products, which is in agreement with our results. The remainder of their list, namely hydroxyacetone, 2-hydroxypropionaldehyde, LGA,Footnote 1 CO and CO2 are missing from our simulations (or only observed in trace amounts). This is reasonable, as they observed the species in question being formed through reaction intermediates. Secondly, Zheng et al. report observing 5-hydroxymethyl furfural (HMF), again formed through an intermediate. We did not observe HMF in our simulations, as would be expected. Moreover, we observed furanose ring structures, which might indicate the presence of HMF precursors (Agarwal et al. 2012), in <1% of our simulations. Lastly, Zheng et al. report the water formation rate to be an increasing function of temperature at temperatures below 2000 K. In our simulations, the rate is roughly constant for all considered temperatures. The distinct observations suggest that the primary water formation mechanisms in cellulose pyrolysis involve either reaction intermediates or inter-chain interactions, or both. Zheng et al. do not propose a mechanism for the water formation, and further comparison is thus not possible. Altogether, the observations of Zheng et al. are in agreement with ours.

Conclusions

We have shown that stochastic ReaxFF-MD simulations can be used to study the thermal decomposition of the cellulose molecule, and to predict the kinetics of the primary decomposition reactions. Our simulations suggest that, in gas-phase conditions, the decomposition begins with random cleavage of the β(1 → 4)-glycosidic bonds. The predicted kinetic parameters, \(A = \left( {1.07 \pm 0.12} \right) \cdot 10^{15}\) s−1 and \(E_{\text{a}} = \left( {171 \pm 2} \right)\) kJ mol−1, are in line with previous ab initio predictions for an initial decomposition step, and within the range of values reported for the global mass loss kinetics of cellulose pyrolysis. Our results imply that the depolymerisation mechanism stays unchanged throughout the considered temperature range from 1400 to 2200 K. We did not observe dependency of the kinetic parameters on molecular weight or initial conformation. Interestingly, our simulations suggest that the glycosidic oxygen prefers to remain bonded towards the non-reducing end of the cellulose molecule. Besides the depolymerisation reactions, we also observed the release of LMWPs directly from the cellulose molecule. Excluding radicals, the most commonly observed products were glycolaldehyde, water, formaldehyde and formic acid, which are all recognised products of cellulose pyrolysis. We did not observe LGA during the simulations, neither as individual molecules nor as LGA-end structures. The general compatibility of our results with the existing knowledge supports further use of ReaxFF-MD in the study of cellulose pyrolysis.

The present work can be extended in a number of directions. Perhaps most importantly, a similar approach could be developed for the solid-phase pyrolysis reactions in crystalline and amorphous cellulose. Moreover, there are certain technical improvements that could be realised. Firstly, it would be feasible to use higher sampling frequencies for the chemical bonding information, and to analyse the resulting data in bond-specific detail. With suitable algorithms, it would then be possible to track individual decomposition reactions and to determine their associated kinetics. Secondly, the stochastic approach could be directly applied to secondary and tertiary reactions. This could be realised by running longer simulations with more reactive events, or by using the products of primary reactions as a starting point. Thirdly, the reliability of the predictions could be improved by shifting the simulation temperatures towards practical pyrolysis conditions. This could be feasible with recently introduced software and hardware acceleration methods, such as the accelerated ReaxFF reactive dynamics (aARRDyn) (Cheng et al. 2014), reactive parallel replica dynamics (RPRD) (Joshi et al. 2013) and the GPU-accelerated version of ReaxFF (Zheng et al. 2013). Together, the suggested extensions would constitute a ReaxFF-MD framework for predicting the primary reaction pathways of cellulose pyrolysis, the associated kinetics and the main features of the product spectrum.