Multiscale Informatics for Low-Temperature ... - ACS Publications

May 6, 2015 - Department of Mechanical Engineering, Department of Chemical Engineering, and Data Sciences Institute, Columbia University,. New York, N...
1 downloads 10 Views 3MB Size
Article pubs.acs.org/JPCA

Multiscale Informatics for Low-Temperature Propane Oxidation: Further Complexities in Studies of Complex Reactions Michael P. Burke,*,†,‡ C. Franklin Goldsmith,‡,§ Stephen J. Klippenstein,‡ Oliver Welz,∥,⊥ Haifeng Huang,∥ Ivan O. Antonov,∥ John D. Savee,∥ David L. Osborn,∥ Judit Zádor,∥ Craig A. Taatjes,∥ and Leonid Sheps∥ †

Department of Mechanical Engineering, Department of Chemical Engineering, and Data Sciences Institute, Columbia University, New York, New York, United States ‡ Chemical Sciences and Engineering Division, Argonne National Laboratory, Argonne, Illinois, United States § School of Engineering, Brown University, Providence, Rhode Island, United States ∥ Combustion Research Facility, Sandia National Laboratories, Livermore, California, United States S Supporting Information *

ABSTRACT: The present paper describes further development of the multiscale informatics approach to kinetic model formulation of Burke et al. (Burke, M. P.; Klippenstein, S. J.; Harding, L. B. Proc. Combust. Inst. 2013, 34, 547−555) that directly incorporates elementary kinetic theories as a means to provide reliable, physics-based extrapolation of kinetic models to unexplored conditions. Here, we extend and generalize the multiscale informatics strategy to treat systems of considerable complexity involving multiwell reactions, potentially missing reactions, nonstatistical product branching ratios, and non-Boltzmann (i.e., nonthermal) reactant distributions. The methodology is demonstrated here for a subsystem of low-temperature propane oxidation, as a representative system for low-temperature fuel oxidation. A multiscale model is assembled and informed by a wide variety of targets that include ab initio calculations of molecular properties, rate constant measurements of isolated reactions, and complex systems measurements. Active model parameters are chosen to accommodate both “parametric” and “structural” uncertainties. Theoretical parameters (e.g., barrier heights) are included as active model parameters to account for parametric uncertainties in the theoretical treatment; experimental parameters (e.g., initial temperatures) are included to account for parametric uncertainties in the physical models of the experiments. RMG software is used to assess potential structural uncertainties due to missing reactions. Additionally, branching ratios among product channels are included as active model parameters to account for structural uncertainties related to difficulties in modeling sequences of multiple chemically activated steps. The approach is demonstrated here for interpreting time-resolved measurements of OH, HO2, n-propyl, i-propyl, propene, oxetane, and methyloxirane from photolysis-initiated low-temperature oxidation of propane at pressures from 4 to 60 Torr and temperatures from 300 to 700 K. In particular, the multiscale informed model provides a consistent quantitative explanation of both ab initio calculations and time-resolved species measurements. The present results show that interpretations of OH measurements are significantly more complicated than previously thoughtin addition to barrier heights for key transition states considered previously, OH profiles also depend on additional theoretical parameters for R + O2 reactions, secondary reactions, QOOH + O2 reactions, and treatment of non-Boltzmann reaction sequences. Extraction of physically rigorous information from those measurements may require more sophisticated treatment of all of those model aspects, as well as additional experimental data under more conditions, to discriminate among possible interpretations and ensure model reliability.

1. INTRODUCTION

fuel and oxygen to carbon dioxide and water in fact proceeds through numerous (sometimes reaching thousands of) intermediate species, all of which are undergoing simultaneous sequential and parallel reactions, each of which occur at rates that can be temperature, pressure, and mixture-composition (T/P/M)

Predictive simulations employing computational fluid dynamics with detailed descriptions of the governing fluid mechanics and chemical kinetics offer promise to accelerate advanced engine development.1−4 Naturally, the utility of these computational design tools ultimately hinges on the accuracy of the submodels used to describe the relevant fluid mechanics and chemical kinetics. With regard to the chemical kinetics, many advanced engine concepts in particular tend to be far more sensitive to fuel chemistry compared to conventional engines. These chemical processes that occur in combustion are both highly complex and highly condition-dependent. Namely, the overall conversion of © XXXX American Chemical Society

Special Issue: 100 Years of Combustion Kinetics at Argonne: A Festschrift for Lawrence B. Harding, Joe V. Michael, and Albert F. Wagner Received: February 9, 2015 Revised: May 3, 2015

A

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

empirical rate parametrization expressions. Theoretical kinetics calculations have in fact often been used for extrapolation of limited experimental measurements in this manner for individual reactions.23−29 Furthermore, the accuracy of ab initio electronic structure calculations has advanced to a sufficient extent that the calculations (in addition to experimental measurements) impose meaningful constraints on the model parameters. A particularly important feature of the constraints imposed by theoretical kinetics calculations is that they span all T/P/M conditions essentially “filling in the gaps” in vast parameter space and allowing reliable extrapolation of limited data to previously unexplored conditions. These ideas largely motivated the multiscale informatics approach.30 In contrast to conventional kinetic models represented by a set of rate parameters within a model structure comprised of parametrized rate expressions, multiscale informed models are instead represented by a set of theoretical kinetics parameters within a model structure comprised of elementary kinetics theories. This set of theoretical kinetics parameters is then informed by theoretical and experimental data across a wide range of scales, including ab initio electronic structure calculations of molecular properties, rate constant determinations for individual reactions, and measured global observables in multireaction systems. The resulting kinetic model consists of a set of theoretical parameters (with constrained uncertainties), which can be related through elementary kinetics models to T/P/M-dependent rate constants (with propagated uncertainties), which in turn can be related through physical models to global observable behavior (with propagated uncertainties). Within this multiscale framework, data from ab initio calculations and experimental measurements are considered simultaneously. In this way, theory guides the interpretation of global data in a manner restricted to be physically meaningful, while experimental measurements impose tight constraints on the uncertain aspects of the theoretical treatment to improve their accuracy both facilitating experimental interpretations and enabling the constraints imposed by limited experimental measurements to be propagated to predictions at different T/P/M conditions. Here, we extend and generalize the multiscale informatics strategy30 to treat systems of substantial complexity, where both “parametric” and “structural” uncertainties are present, as discussed below (cf. section 2 below). This overall modeling strategy was developed to address complex reaction systems with multiwell reactions, potentially missing reactions in the starting mechanism, and potential non-Boltzmann reaction sequences31−40 (cf. section 3.3 below). The methodology is demonstrated for a subsystem of lowtemperature propane oxidation, as a representative system for low-temperature fuel oxidation. Low-temperature propane oxidation has been the subject of considerable attention,26,39−61 especially as a prototype for low-temperature alkane oxidation for which high-level theory is computationally tractable. We focus in particular on oxidation under dilute conditions intended to accentuate R + O2 (and perhaps QOOH + O2) reactions (where R denotes n/i-propyl radicals and QOOH refers to one of the three isomers of hydroperoxy propyl radicals, Figure 1). Thermokinetic parameters for these reactions in particular have been demonstrated to be essential to highaccuracy predictions of 0-D,62−66 nonpremixed counterflow,65 and cool flame55 systems for a variety of fuels. In many earlier studies,26,50−54,56 experimental measurements of time-dependent OH or HO2 concentrations from photolysis-initiated propane oxidation were used to adjust barrier heights of selected

dependent. Given that the reaction pathways and the accentuated portion of their associated rates within T/P/M space vary strongly with the exact local thermodynamic conditions, kinetic models are seldom applied to kinetic situations similar to those probed in limited experimental validation data a difficulty that is only exacerbated by the fact that advanced engines operate at pressure and temperature conditions that are especially challenging to study experimentally. Consequently, extrapolation of kinetic models, particularly to the higher pressures and lower temperatures relevant to advanced engines, is often problematicas illustrated by recent examples for even relatively simple combustion systems where the reaction proceeds through only ∼20−30 reactions.5−7 Such requirements for reliable extrapolation combined with observed deficiencies of model predictions involving extrapolation emphasize the need for improved strategies for kinetic model development. Historically, conventional kinetic models have been comprised of rate parameters that represent the T/P/M-dependent rate constants of the underlying reactions within parametrized rate expressions (e.g., modified Arrhenius, Troe, and PLOG expressions). Nearly all strategies for kinetic model development to date are based on identifying a set of rate parameters informed by rate constant determinations for individual reactions and measurements of global observables in multireaction systems. Most rigorous model development strategies contain elements of the comprehensive modeling philosophy of Westbrook and Dryer8 where models are tested at as wide a range of conditions as possible. The roots of mathematical approaches to combustion model development can in large part be traced to the early work of Frenklach9 and co-workers.10 Therein, solution maps (or “surrogate models”), representing the response of model predictions to variations in active parameters, were employed to optimize A-factors in rate constant expressions to combustion systems measurements.9−14 Later, Frenklach, Wang, and their co-workers15,16 generalized these techniques to enable uncertainty quantification of the complex constraints imposed by global observable targets on the active model parameters. Within the past few years, Turányi and co-workers17−19 and Sheen et al.20 have further advanced these rate-parameter-based approaches to include temperature-dependent rate parameters among the active parameters (in addition to A-factors), where Turányi and co-workers17−19 also consider the constraints imposed by rate constant determinations among the targets (in addition to global observables in combustion systems). Additionally, Prager et al.21 and Cai and Pitsch22 have extended A-factor uncertainty quantification approaches to include consideration of the cross-correlations among rate parameters introduced by the use of rate estimation rules. However, given that the rate constant expressions that form the model structure in rate-parameter-based strategies are simply fitting formulas without direct physical meaning in many cases, models constructed in this manner are essentially empirical. As such, massive amounts of experimental measurements are required in order to sufficiently constrain conventional (rate-parameterbased) kinetic models over all T/P/M conditions of potential interest. This amount of data is available at most for only a handful of combustion systems, such that the reliability of model extrapolation to new conditions in most cases is at best uncertain. To address the difficulties associated with extrapolation, it appears that integrating theoretical kinetics directly into the model development process provides a unique advantage. In such a way, the underlying model structure can be comprised of physically meaningful elementary kinetics theories in lieu of B

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

measurements beyond previous work.26 While model predictions using the rate constants of Huang et al.26 (derived from their OH measurements at 10 Torr) are consistent with these higher pressure measurements at intermediate temperatures, the predictions suggest a much stronger temperature dependence than observed in these new experimentswith model predictions yielding an order of magnitude lower OH at 300 K and 50% higher OH at 733 K (cf. Figure 13 below). In our companion paper, Welz et al.61 also investigated photolysis-initiated propane oxidation at low-pressure, lowtemperature conditions (4 Torr and 530−670 K) but performed measurements employing time-resolved multiplexed photoionization mass spectrometry (MPIMS). This approach provides the advantage of simultaneous time-resolved detection of multiple species, including the stable coproducts linked to formation of OH (oxetane, acetone, propanal, and methyloxirane) and HO2 (propene) from the R + O2 reaction surfaces. These MPIMS measurements serve as powerful complements to measurements of OH and HO2 time profiles by providing more direct constraints on individual sources of OH and HO2 production. Most notably, model predictions using the rate constants of Huang et al.26 (derived from their OH measurements) predict 4 times more methyloxirane (the presumed OH coproduct) than observed in the MPIMS measurements (e.g., see Figure 8 below). Clearly, the OH measurements of the present study and methyloxirane measurements of Welz et al.61 both imply deficiencies in the earlier interpretations.26 Recent calculations from Goldsmith et al.,59 which employed higher level theory for R + O2, do not resolve this disagreement, though they presented the first high-level theoretical kinetics calculations of QOOH + O2. Interestingly, as discussed below, QOOH + O2 emerged as an unexpected key factor in the multiscale informed model results but only after the constraints from the MPIMS measurements are imposed. Below, we use the multiscale informatics approach (1) to derive a consistent interpretation of both earlier and recent data and (2) to identify all relevant parameters that influence predictions and, consequently, interpretations of the experiments. First, we describe the overall methodology, implementation for the low-temperature propane oxidation system, and targets considered for the optimization and uncertainty quantification. Next, we highlight the most important parameters in predictions of select experimental targets and we present a multiscale informed model consistent with ab initio calculations and time-resolved measurements of OH, HO2, n-R, i-R, propene, oxetane, and methyloxirane at pressures from 4 to 60 Torr and temperatures from 300 to 700 K. Interestingly, we find that the interpretations of many earlier measurements, including the OH time profiles at 10 Torr,26 are affected by uncertainties in many more parameters than considered previously. In particular, the multiscale informed model results identify an unexpected role of QOOH + O2 reactions and potential non-Boltzmann reaction sequences in interpretations of these earlier OH measurements.26

Figure 1. Potential energy surfaces for the (a) n-R + O2 (R1), (b) i-R + O2 (R2), and (c) QOOH-1 + O2 (R3) reaction surfaces.59 Symbols are used for propyl (R), propylperoxy (RO2), 3-hydroperoxy-1-propyl (QOOH-1), 3-hydroperoxy-2-propyl (QOOH-2), 2-hydroperoxy-2-propyl (QOOH-3), 3-hydroperoxy-n-propylperoxy (O2QOOH), 2-formyl-ethylhydoperoxide (OQ′OOH), and allylhydroperoxide (Q−HOOH).

transition states in the R + O2 reaction surfaces found to be important under the given experimental conditions, under the assumption that uncertainty in the saddle point energies was the dominant source of uncertainty in the model. Although OH and HO2 are central to radical chain branching, propagation, and termination during low-temperature oxidation,67 measurements of these species alone provide relatively indirect kinetic information, particularly regarding the exact reactions responsible for their production. Previously, Huang et al.26 concluded that R + O2 ↔ OH + methyloxirane was responsible for OH production at their experimental conditions of 10 Torr and 426−698 K. They adjusted the barrier heights for rate-limiting transition states to form OH + methyloxirane (TS3 and TS11, following the nomenclature of Goldsmith et al.,59 cf. Figure 1 and Table 1) to match their measured OH time profiles. No new interpretations of the OH measurements of Huang et al.26 have been provided since their study in 2011. However, new experimental results included in the present paper and in our companion paper61 now provide a means to evaluate the validity of these previous interpretations. As part of the present study, we performed laser-spectroscopic measurements of OH time profiles at 30 Torr and 300−733 K extending the range of temperatures and pressures of OH

2. MULTISCALE INFORMATICS APPROACH The present modeling technique identifies optimized values and quantified uncertainties for a set of theoretical kinetics, rate constant, and physical model parameters that describe the overall kinetic model and physical models for the prediction of global observables. These active parameters are informed by available data from ab initio calculations, rate constant measurements, and C

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 1. List of Reactants, Wells, and Transition States; Energies; and Active Theoretical Kinetics Parametersa energiesb (kcal/mol) reactants, wells, and transition states

a priori

MS-informed

theoretical kinetics parameters

R1 r1 w1 w2 w3 TS1 TS2 TS3 TS4 TS5 TS6 TS7 TS-R3 TS8 TS9

n-R + O2 n-RO2 QOOH-1 QOOH-2 n-R + O2 ↔ n-RO2 n-RO2 ↔ QOOH-1 n-RO2 ↔ QOOH-2 n-RO2 ↔ HO2 + propene n-RO2 ↔ OH + propanal QOOH-1 ↔ OH + oxetane QOOH-1 ↔ OH + CH2O + C2H4 QOOH-1 ↔ R3c QOOH-2 ↔ HO2 + propene QOOH-2 ↔ OH + methyloxirane

r2 w4 w5 TS10 TS11 TS12 TS13 TS14 TS15

i-R + O2 i-RO2 QOOH-3 i-R + O2 ↔ i-RO2 i-RO2 ↔ QOOH-3 i-RO2 ↔ HO2 + propene i-RO2 ↔ OH + acetone QOOH-3 ↔ HO2 + propene QOOH-3 ↔ OH + methyloxirane

r3 w6 TS16 TS18 TS19

QOOH-1 + O2 O2QOOH QOOH-1 + O2 ↔ O2QOOH O2QOOH ↔ HO2 + Q−HOOH O2QOOH ↔ OH + OQ′OOH*

−33.1 −17.7 −20.0

−33.6 −17.6 −20.4

−8.1 0.2 −2.8 8.6 3.6 9.0

−7.3 −0.4 −4.1 8.5 2.1 8.9

v′low,r1 (0.3) Ew1 (1 kcal/mol), vlow,w1 ′ (0.3), αw1 ′ (0.4), βw1 (0.3) Ew2 (1 kcal/mol), vlow,w2 ′ (0.3), αw2 ′ (0.4), βw2 (0.3) Ew3 (1 kcal/mol), v′low,w3 (0.3), α′w3 (0.4), βw3 (0.3) f ′VRC‑TST,TS1 (0.4) E†TS2 (3 kcal/mol), v†low,TS2 ′ (0.3), vimag,TS2 ′ (0.3) E†TS3 (3 kcal/mol), v†low,TS3 ′ (0.3), vimag,TS3 ′ (0.3) E†TS4 (3 kcal/mol), v†low,TS4 ′ (0.3), v′imag,TS4 (0.3) E†TS5 (3 kcal/mol), v†low,TS5 ′ (0.3), v′imag,TS5 (0.3) E†TS6 (3 kcal/mol), v†low,TS6 ′ (0.3), vimag,TS6 ′ (0.3) E†TS7 (3 kcal/mol), v†low,TS7 ′ (0.3), vimag,TS7 ′ (0.3)

−2.4 −7.5

−2.5 −7.4

E†TS8 (3 kcal/mol), v†low,TS8 ′ (0.3), vimag,TS8 ′ (0.3) E†TS9 (3 kcal/mol), v†low,TS9 ′ (0.3), vimag,TS9 ′ (0.3)

R2 −34.6 −17.6

−35.1 −17.5

1.6 −4.3 6.0 −1.3 −5.0

1.0 −5.1 5.9 −1.2 −5.0

v′low,r2 (0.3) Ew4 (1 kcal/mol), vlow,w4 ′ (0.3), αw4 ′ (0.4), βw4 (0.3) Ew5 (1 kcal/mol), vlow,w5 ′ (0.3), αw5 ′ (0.4), βw5 (0.3) f ′VRC‑TST,TS10 (0.4) E†TS11 (3 kcal/mol), v†low,TS11 ′ (0.3), v′imag,TS11 (0.3) E†TS12 (3 kcal/mol), v†low,TS12 ′ (0.3), vimag,TS12 ′ (0.3) E†TS13 (3 kcal/mol), v†low,TS13 ′ (0.3), vimag,TS13 ′ (0.3) E†TS14 (3 kcal/mol), v†low,TS14 ′ (0.3), v′imag,TS14 (0.3) E†TS15 (3 kcal/mol), v†low,TS15 ′ (0.3), v′imag,TS15 (0.3)

R3 −33.7

−33.6

−3.7 −11.8

−3.7 −14.1

(vlow,r3 ′ = vlow,w2 ′ ) Ew6 (1 kcal/mol), v′low,w6 (0.3), α′w6 (0.4), βw6 (0.3) f VRC‑TST,TS16 ′ (0.4) E†TS18 (3 kcal/mol), v†low,TS18 ′ (0.3), vimag,TS18 ′ (0.3) E†TS19 (3 kcal/mol), v†low,TS19 ′ (0.3), vimag,TS19 ′ (0.3)

a Note that ′ indicates ln( ) of the quantity. Uncertainties for targets used for uncertainty quantification are listed in ( ). See text (section 3.2) for a description of the active parameters and their symbols listed above. See the Figure 1 caption for a description of the chemical nomenclature and section 3.3 for a description of the * symbol. bAll energies are relative to the R + O2 or QOOH + O2 reactants for each reaction. cTS-R3 refers to the irreversible sink corresponding to the bimolecular QOOH-1* + O2 reaction (section 3.3).

global measurements (e.g., time-dependent species concentrations, product yields, ignition delay times, and flame speeds). The modeling technique and its implementation presented here constitute a major extension and generalization of our earlier implementation of the multiscale informatics approach30 to treat systems of substantial complexity, where both “parametric” and “structural” uncertainties are present. In this context, parametric uncertainties refer to quantitative uncertainties in the model description due to uncertainties in the model parameters; structural uncertainties refer to qualitative uncertainties in the model description due to assumptions/limitations of the model structure itself (also termed “model inadequacy”68 or “model discrepancy error”, 69 e.g., missing reactions or assumption of reactants in a Boltzmann distribution, cf. sections 3.1 and 3.3 below). The impact of structural uncertainties due to these approximately treated aspects of the model is quantified here through exploration of limiting cases. In that way, the quantified impact of structural uncertainties can then be assessed in screening procedures to eliminate targets affected by these uncertainties, so that projection of information onto the model parameters in the optimization is minimally contaminated by these approximately treated aspects of the model.

The optimization procedure identifies a set of active model parameters, Xj, that minimizes the least-squared error to the set of weighted equations

Fi(Xj) = Yt, i ± Zi

(1)

where Fi(Xj) is the model prediction for a given Xj; Yt,i is the target value; and the weighting factor, Zi = σi/Wi, is equal to the uncertainty, σi, divided by an additional data set weighting factor, Wi. The data set weighting factor, Wi, is used to treat data sets with multiple measurements, as discussed in section 3 below. Employing a “surrogate model” of the system (as in solution mapping methods9) to represent the response of model predictions to variations in the active parameters in the neighborhood of X̃j Fi(Xj) ≈ Fi(X̃j ) +

∑ Sij(Xj − X̃j) j

(2)

with Sij = {∂Fi /∂Xj} Xj = X̃j , eq 1 can be expressed as the matrix equation

∑ Sij(Xj − X̃j) = Yi ± Zi j

(3)

with Yi = Yt, i − Fi(X̃j ). D

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Similar to the Active Thermochemical Tables,70 a standard weighted least-squares procedure is used to calculate the optimized values for the parameters and covariance matrix representing their constrained uncertainties as imposed by the set of targets. Using the nominal values as an initial guess, X̃j°, the following iterative procedure is implementedXj is calculated using eq 3, Fi and Sij are recalculated using the updated Xj → X̃j values, and the process is repeated until convergence is achieved. The covariance matrix is used to estimate propagated prediction uncertainties. Such an approach approximates all target distributions as independent and normal; it also approximates the response surface surrounding the optimized values as linear. (Uncertainty estimates that include higher-order terms or allow different assumed forms of the probability distribution function could be calculated,15,16,18,28,71−80 though they were not considered for the present purposes.) Uncertainty values here are intended to represent two standard deviations, though the fact that uncertainties are often not reported and/or not straightforward to estimate in both kinetics experiments (see section 2.6 of ref 81) and theoretical calculations prevents a more precise definition. Similarly, uncertainties in many model parameters are likely correlated, which can impact optimized values and quantified uncertainties,21,22 though the minimal information about their probability distribution functions also prevents proper quantification of those effects. As discussed below, the active parameters, Xj, include theoretical kinetics parameters, rate constant parameters, and physical model parameters. The targets, Yt,i, include (I) molecular properties from ab initio calculations, (II) rate constant measurements, (III) global measurements, and (IV) reported values for the experimental conditions. Target classes I and IV effectively serve as regularization terms, which impose the prior distributions on all model parameters. Consequently, the system is sufficiently constrained by target classes I and IV alone, and all other data impose further constraints on the model parameters. The a priori model is comprised of all parameters at their nominal values with uncertainties imposed by target classes I and IV alone. Multiscale informed models are comprised of all parameters optimized to and with uncertainties constrained by all target classes (I−IV). 2.1. Target Class I: Ab Initio Electronic Structure Calculations. Using molecular properties from ab initio electronic structure calculations (or conceivably spectroscopic measurements) as targets is accomplished via a trivial mapping Sij = δij

ab initio calculation, and the optimized model prediction error for the ith target may be nonzero, but the total model prediction error for all targets will be minimized. The situation is entirely analogous for target class IV discussed below. For example, a model parameter describing the initial temperature is analogous to a model parameter describing the barrier height; the reported initial temperature as measured by a thermocouple is analogous to the barrier height as calculated ab initio, and the uncertainty in the thermocouple measurement for the temperature is analogous to the uncertainty in the ab initio calculation for the barrier height. 2.2. Target Class II: Rate Constant Measurements. Using rate constant measurements for reaction n at given T/P/M conditions, kt,n(T,P,M), as targets is accomplished via a mapping from theoretical kinetics parameters to rate constant predictions, kp,n(T,P,M), where the subscripts t and p are used to denote the target and predicted rate constants. The mapping is based on an appropriate elementary kinetics theory for calculating the rate constants (e.g., transition state theory, master equation simulations, quantum scattering theory): Sij =

∂ ln k p, n(Ti , Pi , Mi) ∂Xj

(5)

Rate constant measurements impose constraints on combinations of theoretical kinetics parameters describing a particular reaction surface. Experimental measurements used to derive rate constants conceivably can be included as either rate constant measurements, class II targets, or global measurements, class III targets. Given that even experiments carefully designed to isolate a single elementary reaction and derive its rate constant involve uncertain rate constants for “secondary” reactions and/or uncertain physical model parameters (e.g., initial temperature, pressure, or mixture composition in the experiment) that influence interpretations to some nonzero extent (albeit sometimes negligible), the information content of the experimental measurements is more accurately quantified through their inclusion as global measurements (provided that sufficient details and raw data are available). Considering the measured observables from these experiments in that manner allows consideration of these uncertainties in the data interpretation in the kinetic and physical models. 2.3. Target Classes III and IV: Global Measurements and Experimental Conditions. Using observables from global experiments as targets is accomplished via a twofold mapping: the mapping from theoretical kinetics parameters to rate constant predictions as discussed above and a mapping from rate constants to predictions of global observables:

(4)

where δij is the Kronecker delta. For example, suppose that Xj is the model parameter describing the barrier height for a particular transition state and that the ith constraint is imposed by an ab initio calculation of that barrier height: Xj = Yt,i ± Zi, where Yt,i is the ab initio calculated value of the barrier height and Zi is the estimated uncertainty in the ab initio calculation, σi. If predictions of all other targets do not depend on that particular barrier height, the optimized jth model parameter (describing that barrier height) will be identical to the ab initio calculated barrier height, its constrained uncertainty will be equal to the estimated uncertainty in the ab initio calculation, and the optimized model prediction error for the ith target will be zero. If predictions of other targets depend on that particular barrier height (and thus provide further information about that barrier height), the optimized jth model parameter (describing that barrier height) may have a different value than the ab initio calculated value, its constrained uncertainty will be lower than the uncertainty in the

Sij =

∑ n

∂ ln k p, n(T , P , M ) ∂Fi ∂ ln k p, n(T , P , M ) ∂Xj

(6)

Measurements of global observables (e.g., time-dependent species concentrations or ignition delay times) impose constraints on combinations of many parameters, including theoretical kinetics parameters for more than one reaction surface, secondary reaction rate parameters, and physical model parameters. Since interpretations of experimental measurements are often influenced by uncertainties in physical model parameters (e.g., initial conditions),30,71,73 uncertain physical model parameters are included here among the active model parameters using their E

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 2. Schematic for the set of entirely automated scripts that construct the overall model representation (including kinetic and physical models) and its prediction for each observable from a given set of active parameters.

3.1. RMG-Generated Kinetic Model. A detailed chemical mechanism was compiled using RMG (Reaction Mechanism Generator) software,82 which couples direct input for the rate constants in a core model and a supplemental rate constant library with automated rate-rule-based estimates for other initially “missing” reactions. The core model was constructed by layering various submodels in the following order: the H2O2 mechanism from Burke et al.;30 the H2 mechanism from Burke et al.;83 R + O2 and QOOH + O2 reactions from Goldsmith et al.;59 and the C0−C3 mechanism from Miller et al.84 The methyl formate mechanism of Dooley et al.,85 which contains low- and high-temperature treatment of C0−C3 species, was used as a supplemental reaction rate constant library. A representative set of conditions in temperature, pressure, and mixture composition for n-R/i-R/C3H8/O2/He mixtures is then explored to search for new species and reactions to add to the core model. When an identified reaction is in the supplemental reaction rate constant library, it is added with the rate constant extracted from the library; otherwise, it is added with a rate constant estimated using RMG rate rules.82 The pressureindependent version of the program was used here (to avoid convergence issues experienced in the pressure-dependent version). Reactions for which automated rate constant estimates appeared unphysical were removed from the model. (In particular, reactions of the form X + YO2 = YOOH + X−H for X and Y = H, CH3, C2H5, n-C3H7, and i-C3H7 were removed because their automatically estimated rate constants based on the same rate rule were 3 orders of magnitude larger than found in recent theoretical studies on C2H5 and n/i-C3H7 + HO2.88 In addition, O + OH = HO2, n/i-C3H7O2 + CH3 = n/ i-C3H7OOCH3, and reactions of the forms X + HO2 = XOOH and XO2 + H = XOOH for X = H, CH3, C2H5, n-C3H7, and i-C3H7 were removed because automatically estimated high-pressurelimit rate constants were expected to be significant overestimates at the present conditions of interest.) Because photolytically generated Cl atoms were used to initiate propane oxidation in

nominal values as targets through eq 4. Treating physical model parameters as active parameters in this way is intended to achieve two purposes. First, treating uncertain experimental conditions as active model parameters is a more representative means of treating the uncertainty they introduce than simply propagating it to, and commensurately inflating, the uncertainty in the experimental observable, as is commonly done. For example, if an initial temperature in the experiment is different from the reported nominal temperature, it introduces systematic errors in the observables from that experiment that are better represented through considering the uncertainty in the initial temperature than increasing the uncertainty in the observable by an amount estimated through uncertainty propagation. Second, treating uncertain experimental conditions as optimization parameters conceivably allows for the compilation of posterior statistics on the experimental conditions (provided there are sufficiently many experiments of sufficient diversity). For example, if many measurements are performed for a wide variety of kinetic systems with the same experimental apparatus and instrumentation, inspection of the optimized values for the initial temperatures may reveal that the measurements have a particular systematic error and random error.

3. IMPLEMENTATION The overall model representation (including kinetic and physical models) and model prediction for the ith observable, Fi(Xj), is constructed on the fly through a series of entirely automated scripts (see Figure 2). The kinetic model is comprised of the following reaction sets: (1) active reactions for which rate constants are described by theoretical kinetics parameters (R1−R3, Figure 1, Table 1); (2) active secondary reactions (ASRs) for which rate constants are described by Arrhenius parameters (ASR1−ASR32, Table 2); (3) static secondary reactions (SSRs), inclusion of which are described by a binary inclusion parameter, ISSR; (4) static missing secondary reactions (MSRs), inclusion of which are described by a binary inclusion parameter, IMSR. F

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 2. List of Active Parameters for Secondary Reactions, Structural Assumptions, and Physical Modelsa active secondary reactions ASR1 Cl + C3H8 ↔ HCl + n-R ASR2 Cl + C3H8 ↔ HCl + i-R ASR3 OH + C3H8 ↔ H2O + n-R ASR4 OH + C3H8 ↔ H2O + i-R ASR5 n-R + HO2 ↔ OH + n-RO ASR6 i-R + HO2 ↔ OH + i-RO ASR7 Cl + n-R ↔ HCl + propene ASR8 OH + Cl2 ↔ HOCl + Cl ASR9 Cl2 + n-R → Cl + n-RCl ASR10 Cl2 + i-R → Cl + i-RCl ASR11 n-RO2 + HO2 ↔ n-ROOH + O2 ASR12 i-RO2 + HO2 ↔ i-ROOH + O2 ASR13 CH3 + HO2 ↔ CH3O + OH ASR14 HO2 + HO2 ↔ HOOH + O2 ASR15 n-RO + OH ↔ n-ROOH ASR16 i-RO + OH ↔ i-ROOH ASR17 OH + n-RO2 ↔ HO2 + n-RO ASR18 OH + i-RO2 ↔ HO2 + i-RO ASR19 n-RO2 + n-RO2 ↔ n-RO + n-RO + O2 ASR20 n-RO2 + CH3 ↔ n-RO + CH3O ASR21 n-RO2 + n-R ↔ n-RO + n-RO ASR22 n-RO2 + i-R ↔ n-RO + i-RO ASR23 i-RO2 + i-RO2 ↔ i-RO + i-RO + O2 ASR24 n-RO2 + i-RO2 ↔ n-RO + i-RO + O2 ASR25 i-RO2 + CH3 ↔ i-RO + CH3O ASR26 i-RO2 + i-R ↔ i-RO + i-RO ASR27 i-RO2 + n-R ↔ i-RO + n-RO ASR28 n-RO2 + n-RO2 ↔ propanal + n-propanol + O2 ASR29 i-RO2 + i-RO2 ↔ acetone + i-propanol + O2 ASR30 n-RO2 + i-RO2 ↔ acetone + n-propanol + O2 ASR31 i-RO2 + n-RO2 ↔ propanal + i-propanol + O2 ASR32 CH3 + C2H4 (+M) ↔ n-R (+M) model assumptions non-Boltzmann branching ratio product-excitation branching ratio static secondary reactions missing secondary reactions global experiment E1−E4 (COCl)2/Cl/C3H8/O2/He

rate constant parameters A′ASR1, EASR1 (0.2 kcal/mol) A′ASR2, EASR2 (0.2 kcal/mol) A′ASR3, EASR3 A′ASR4, EASR4 A′ASR5 (0.7) A′ASR6 (0.7) A′ASR7 (1.6) A′ASR8 (0.4) A′ASR9 (1.6) A′ASR10 (1.6) A′ASR11 (1.6) A′ASR12 (1.6) A′ASR13 (0.4) A′ASR14 (0.7) A′ASR15 (1.6) A′ASR16 (1.6) A′ASR17 (1.6) A′ASR18 (1.6) A′ASR19 (1.6) A′ASR20 (1.6) A′ASR21 (1.6) A′ASR22 (1.6) A′ASR23 (1.6) A′ASR24 (1.6) A′ASR25 (1.6) A′ASR26 (1.6) A′ASR27 (1.6) A′ASR28 (1.6) A′ASR29 (1.6) A′ASR30 (1.6) A′ASR31 (1.6) A′ASR32 (0.7) structural parameters γn‑B γp‑e ISSR IMSR physical model parameters T′e (0.02), P′e (0.02), M′O2,o,e (0.05), M′C3H8,o,e (0.05)

E5−E11

Te′ (0.02), Pe′ (0.02), M′O2,o,e (0.05), M′C3H8,o,e (0.05)

M′(COCl)2,o,e (0.05), MCl,o,e ′ (0.3), σOH,e ′ (0.1) Cl2/Cl/C3H8/O2/He

e = 1...4

M′Cl2,o,e (0.05), M′Cl,o,e (0.3), σ′OH,e (0.1)

e = 5...11

E12−E13

Cl2/Cl/C3H8/O2/He

Te′ (0.02), Pe′ (0.02), M′O2,o,e (0.05), M′C3H8,o,e (0.05)

E14−E15

C6F5C4H9/n-R/O2/He

T′e (0.02), P′e (0.02), M′O2,o,e (1.6), M′n‑R,o,e (1.6)

e = 14...15

E16

precursor/i-R/O2/He

T′e (0.02), P′e (0.02), M′O2,o,e (1.6), M′i‑R,o,e (1.6)

e = 16

E17−E19

(COCl)2/Cl/C3H8/O2/He

Te′ (0.02), Pe′ (0.02), M′O2,o,e (0.05), M′C3H8,o,e (0.05)

M′Cl2,o,e (0.05), MCl,o,e ′ (0.7), σ′HO2,e (0.1)

e = 12...13

M′(COCl)2,o,e (0.05), MCl,o,e ′ (0.3), σpropene,e ′ (0.2), σoxetane,e ′ (0.2) σ′acetone,e (0.2), σ′propanal,e (0.2), σ′methyloxirane,e (0.2)

e = 17...19

Note that ′ indicates ln( ) of the quantity. Uncertainties for targets used for uncertainty quantification are listed in ( ). See text (section 3.2) for a description of the active parameters and their symbols listed above. a

many of the experiments considered, reactions involving Cl-containing species and other secondary reactions from DeSain et al.52 were appended to the RMG-generated model. After the above-mentioned construction of the model, a few reactions suspected to be of importance to modeling the observables considered as targets were then assigned rate constants

from more recent investigations and/or added to the kinetic model. These included ASR1 and ASR2 (Cl + C3H8 ↔ HCl + n-R/i-R),86 ASR3 and ASR4 (OH + C3H8 ↔ H2O + n-R/i-R),87 ASR5 and ASR6 (n-R/i-R + HO2 ↔ OH + n-RO/i-RO),88 n-R and i-R decomposition,27 and n-RO and i-RO decomposition.89 Additional RO 2 + RO2 ↔ R −HO + ROH + O 2 reactions G

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Active physical model parameters included initial temperatures, Te, initial pressures, Pe, initial mole fractions, Me,o, and absorption/photoionization cross sections, σ, for each species (Table 2). Structural parameters included the static secondary reaction inclusion factor, ISSR, and the missing secondary reaction inclusion factor, IMSR, discussed in the section above and nonBoltzmann reactant distribution factors, γn‑B, and product excitation factors, γp‑e, discussed in the section below. 3.3. Non-Boltzmann Reaction Sequences in Second O2 Addition. The n-R + O2 and QOOH-1 + O2 reaction surfaces are given a different treatment than that used in Goldsmith et al.59 in order to address the potential for non-Boltzmann reaction sequences31−40 (Figure 3). Namely, given that the O2

(ASR28−ASR31) were included with the same rate constants for the RO2 + RO2 ↔ RO + RO + O2 reactions (ASR19, ASR23, and ASR24) of the model of Dooley et al.,85 as motivated by the nearly equal branching ratios between RO + RO + O2 and R−HO + ROH + O2 channels recommended by Atkinson et al.86 at 298 K. The reaction rate expression for vinoxy + O2 ↔ CH2O + CO + OH was adopted from Dooley et al.85 Reactions of vinoxy + O2 yielding all other products were removed, as motivated by a recent theoretical study90 that suggests CH2O + CO + OH is the primary product channel at low pressures and/or high temperatures. Thermochemistry was adopted from a recent compilation of Goldsmith et al.91 based on bond-additivity-corrected ab initio calculations. Heats of formation were replaced with values from the Active Thermochemical Tables of Ruscic et al.92 for each of the species where available. Thermochemical values for all other species (mostly Cl-containing species) were taken from the database of Goos et al.93 with updated treatment of ClO.94 Finally, the reactions of the n-R + O2 (R1), i-R + O2 (R2), and QOOH-1 + O2 (R3) systems of Table 1 and active secondary reactions of Table 2 were replaced with the active model treatment based on the active parameters described below. Reactions for which rate constants are extracted from published kinetic models belong to the set of secondary reactions (ASRs, SSRs). Reactions for which rate constants are RMG estimates belong to the set of missing secondary reactions (MSRs). Inclusion of the static secondary reaction and missing secondary reaction sets are governed by binary inclusion parameters, ISSR and IMSR. The active secondary reaction set (ASR), which is included for all predictions, was chosen such that it encompassed all secondary reactions that contributed to predictions of the targets considered. For nearly all targets (and especially those used in the optimization), predictions using only the ASR set (ISSR = 0, IMSR = 0) and predictions using all secondary reactions (ISSR = 1, IMSR = 1) were essentially identical. 3.2. Active Parameters. For the multiwell reaction surfaces of n-R + O2 (R1), i-R + O2 (R2), and QOOH-1 + O2 (R3), important uncertain parameters corresponding to each reactant, well, and transition state of the theoretical kinetics calculations were treated as active model parameters (Table 1). For each set of bimolecular reactants, these included a scaling factor for the lowest five harmonic frequencies of R or QOOH, vlow,r, to represent uncertainties in the partition function calculation. For each well, these included a scaling factor for the lowest five harmonic frequencies of the well, vlow,w, to represent uncertainties in the partition function calculation; the well depth relative to R + O2 or QOOH + O2 reactants, Ew; and scaling factors for the magnitude and temperature exponent of the energy transferred per collision for each well, αw and βw, according to ΔEd,w = αwTβw. For each tight transition state, these included a scaling factor for the lowest five harmonic frequencies of the transition state, v†low, to represent uncertainties in the partition function calculation; the barrier height of the transition state relative to R + O2 or QOOH + O2 reactants, E†; and a scaling factor for the imaginary frequency in the Eckart tunneling correction, vimag. For each barrierless transition state, these included a variable reaction coordinate transition state theory (VRC-TST) correction factor to account for uncertainty in the number of states in the transition state, f VRC‑TST. For the active secondary reactions (ASR1−ASR32), Arrhenius parameters were treated as active parameters for their utility for target-screening purposes (Table 2).

Figure 3. Schematic of the principal non-Boltzmann reaction sequence for n-R + O2 (+O2). Symbols are used for 2-formyl-ethoxy (OQ′O) and as indicated in the caption for Figure 1. Not all pathways considered are shown above (see Figure 1).

mole fraction in some of the considered experimental targets reaches ∼15%, there is significant potential for reactive collisions with O2 to occur on similar time scales as the unreactive, thermalizing collisions necessary to achieving a Boltzmann distribution in the rovibrational modesviolating a critical assumption underlying all (thermal) kinetic models that bimolecular reactions only occur between Boltzmann-distributed reactants. Of particular concern in the present system, QOOH-1* is initially formed with significant rovibrational excitation. (Here, starred species names, e.g., QOOH-1*, are used to refer to molecular ensembles whose rovibrational energy may or may not follow a Boltzmann distribution, whereas unstarred species names, e.g., QOOH-1, are reserved for molecular ensembles whose rovibrational energy is well described by a Boltzmann distribution.) As discussed in greater detail elsewhere,39,40 when n-R reacts with O2, it forms a rovibrationally excited n-RO2*, which in turn can isomerize to an excited QOOH-1*. QOOH-1* formed in this manner can react with O2 prior to stabilization, thus forming an excited O2QOOH*. O2QOOH* at higher energies decomposes more rapidly to OH + OQ′OOH* with more energy contained in the fragments (than does O2QOOH* at lower energies). OQ′OOH* at sufficiently high energies can decompose directly to OH + OQ′O*. Even (Boltzmanndistributed) OQ′O proceeds rapidly thereafter through a series of reactions involving another O2 addition to ultimately form OH + 2CH2O + CO. This overall non-Boltzmann reaction sequence can occur with different rates and yield different products than typical thermal sequences. Here, reaction of QOOH-1* with O2 is included within the reaction master equation for n-R + O2 by adding an irreversible H

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 4. Procedure for the assignment of branching ratios for effective phenomenological reactions (describing non-Boltzmann reaction sequences) from various thermal reactants (indicated in bold on the left side) to thermal products (indicated in bold on the right side). Overall effective rate constants from each thermal reactant to each thermal product are calculated by taking the product of the rate constant for the irreversible sink corresponding to the second O2 addition process (e.g., kR+O2(+O2)) and the appropriate branching ratios defined above.

detailed calculation strategies.39,40) It is perhaps worth noting that the exact values that the non-Boltzmann and productexcitation factors, γn‑B and γp‑e, take on within their range of limiting cases strongly depend on temperature, pressure, O2 mole fraction, and other conditions.39,40 3.4. Automated Model Construction and Prediction. In automated construction of the overall model representation and model prediction for the ith observable, Fi(Xj), a series of scripts is used (see Figure 2). The kinetics input files for the three master equation systems (R1−R3) are generated from the set of active parameters; the kinetics code (VARIFLEX97) is called to calculate phenomenological rate constants; the kinetics output files are parsed and combined with other active parameters to construct chemical mechanism input files; the physical model input files are constructed from the set of active parameters; the physical model code (SENKIN) is called to calculate timedependent species concentrations; and the physical model output files are parsed and combined with other active parameters to calculate the model prediction for the ith observable, Fi(Xj). Brute force sensitivity coefficients of the model prediction for a given observable to each active parameter were calculated to construct the solution map. Theoretical kinetics calculations were performed using generalized transition state theory98 and multiwell master equation approaches99 in an updated version of the VARIFLEX code.97 Phenomenological rate constants were extracted using the long-time method.99 To handle situations where reactions (unimolecular99 or bimolecular39) are occurring on the same time scale as thermalization, two species within the master equation are considered to be “merged”100 when the eigenvalue corresponding to their conversion is within a factor of 10 of the lowest eigenvalue corresponding to internal energy relaxation. Physical model calculations were performed using the homogeneous, constant-pressure model in SENKIN.

product channel, TS-R3, with an energy-independent rate equal to k3,tot(T) [O2], similar to other studies.31,36−38 Branching ratios among products in effective phenomenological reactions involving this channel are treated differently for different reactants (e.g., n-R + O2 (+O2) vs n-RO2 (+O2)) to account for their different levels of rovibrational excitation relative to the QOOH-1 well (Figures 3 and 4). For example, QOOH-1* + O2 reaction that occurs via phenomenological reactions from n-RO2, i.e., n-RO2 (+O2), is assumed to yield the same products as from (thermal) QOOH-1, i.e., QOOH-1 + O2. (n-RO2 (+O2) can also contribute to non-Boltzmann reaction sequences but significantly less so compared to n-R + O2 (+O2).39) However, QOOH-1* + O2 reaction that occurs via phenomenological reactions from n-R + O2, i.e., n-R + O2 (+O2), is assumed to yield a combination of OH + OH + OQ′O and the products from (thermal) QOOH-1, i.e., QOOH-1 + O2, characterized by a nonBoltzmann reactant branching factor, γn‑B. (The decomposition rate of O2QOOH* → HO2 + Q−HOOH is sufficiently small at all energies relative to the total decomposition rate of O2QOOH* that its contribution to the non-Boltzmann reaction sequences is not considered here.) This branching factor, γn‑B, is treated as a structural parameter to assess the role of non-Boltzmann reactant distribution effects in the n-R + O2 (+O2) reaction (Figure 4). Similarly, branching ratios among products in effective phenomenological reactions involving O2QOOH* decomposition to OH + OQ′OOH* are treated differently for different reactants to account for the fact that their different levels of rovibrational excitation relative to the O2QOOH well produce different levels of rovibrational excitation in the nascent OQ′OOH* (Figures 3 and 4). O2QOOH* decomposition to OH + OQ′OOH* that occurs via phenomenological reactions from thermal O2QOOH is assumed to yield only OH + OQ′OOH. (O2QOOH thermal decomposition can also contribute to non-Boltzmann reaction sequences but significantly less so compared to chemically activated O2QOOH*.40) However, O2QOOH* decomposition to OH + OQ′OOH* that occurs via phenomenological reactions from QOOH-1 + O2 is assumed to yield a combination of OH + OQ′OOH and OH + OH + OQ′O with a branching factor, γp‑e. The branching factor, γp‑e, is treated as a structural parameter to assess the role of product excitation effects in O2QOOH* decomposition from the QOOH-1 + O2 reaction (Figure 4). Although these structural parameters in the above-mentioned treatment do not capture the full complexity of these effects, they are useful here to assess whether, and under what conditions, more detailed calculations39,40 (enabled by a new master equation code95,96) are likely necessary. (In fact, the results from the present study largely motivated the development of the

4. SPECIFICATION OF THE TARGETS 4.1. Target Class I: Ab Initio Electronic Structure Calculations. Molecular properties from the bond-additivitycorrected ab initio electronic structure calculations of Goldsmith et al.59 were used as targets with the uncertainties given in Table 1. Uncertainties in the various theoretical parameters were given estimated typical values based on the level of electronic structure theory and level of conformational treatment. While a detailed consideration of correlations among theoretically calculated molecular properties is outside the scope of the present study, correlation constraints on properties of analogous transition states on the n-R + O2 and i-R + O2 reaction surfaces for HO2 + propene channels (TS4, TS12), OH + methyloxirane I

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Table 3. List of Optimization Targets Considereda

channels (TS3, TS11), and OH + propanal/acetone (TS5, TS13) were also included (viz., ΔE†TS4 − ΔE†TS12 = 0 ± 1 kcal mol−1, Δv†low,TS4 ′ − Δv†low,TS12 ′ = 0 ± 0.1, and Δv′imag,TS4 − Δvimag,TS12 ′ = 0 ± 0.1) to account for some correlated errors in their calculation. For an analysis on the relationship of theoretical parameter variation on the rate constants of multiwell reactions, we refer the reader to the work of Goldsmith et al.101 It is worth noting here that a 20% difference in the lowest five harmonic frequencies yields ∼50−80% differences in the partition functions under the conditions investigated, such that the uncertainties assigned here for v†low are in accordance with the factor of 2 reported uncertainties in the torsional treatment.59 4.2. Target Class II: Rate Constant Measurements. Rate constant measurements for ASR3 and ASR4,87 k298K and Ea for ASR1 and ASR2,86 and nominal values with estimated uncertainty bounds for ASR5−ASR32 were included as targets (Tables 2 and 3). For data sets for rate constant measurements containing multiple data points, weighting factors were taken to be Wi = 1/N0.5, where N is the number of points in the data setequivalent to counting the entire data set as one datum. It is worth noting that individual constraints on the Arrhenius parameters A′ASR1, A′ASR2, A′ASR3, EASR3, A′ASR4, and EASR4 were not necessary, since sufficient constraints from rate constant measurements (Table 3) were included. 4.3. Target Classes III and IV: Global Measurements and Experimental Conditions. Experimental data treated as measured global observables (refs 26, 44, 45, 50, and 61 and present work) and the associated nominal values for physical model parameters in Table 3 (listed in Table S2 in the Supporting Information) were considered as targets. These data include measured OH time profiles in (COCl)2/C3H8/O2/He mixtures at 10 Torr following pulsed photolysis to produce Cl;26 measured OH time profiles in Cl2/C3H8/O2/He mixtures at 30 Torr following pulsed photolysis to produce Cl (present study, see below); measured HO2 time profiles in Cl2/C3H8/O2/ He mixtures at ∼50 Torr following pulsed photolysis to produce Cl;50 n-R logarithmic decay times (i.e., time constants) of n-R precursor/O2/He mixtures at 1−7 Torr following pulsed photolysis to produce n-R;44,45 i-R logarithmic decay times of i-R precursor/O2/He mixtures at 1−4 Torr following pulsed photolysis to produce i-R;45 and measured time profiles for various species using MPIMS in (COCl)2/C3H8/O2/He mixtures at 4 Torr following pulsed photolysis to produce Cl.61 While the n-R and i-R logarithmic decay times were originally reported as second-order rate constants for R + O2,44,45 they are treated here as global observables to allow model quantification of any possible deviations from second-order decay behavior (which indeed were predicted to be minimal using the present models) for these complicated multiwell reaction surfaces. The precursors for n-R (C6F5C4H9,44 1-C3H7I,45 1-C4H9ONO,45 and 1-C3H7NO245), the precursors for i-R (2-C3H7I,45 2-CH3-1C3H6ONO,45 and 2-C3H7NO245), and the associated coproducts of R were not included in the model predictions. Although nine different species were quantified (and even more detected) in the MPIMS measurements,61 we focus here on the stable primary products from the R + O2 reactions (propene, oxetane, acetone, propanal, and methyloxirane). (Note that the oxetane concentration was below the detection limit at 530 K such that oxetane could not be quantified. Also note that the temporal instrument response function applied to the data reported in our companion paper61 was not applied to the data used as targets and plotted below.)

I. Ab Initio Calculations see Table 1 II. Rate Constant Measurements T (K) P (Torr) 298 300−800 III. Global Experiments T (K) P (Torr) E1 426 10 E2 507 10 E3 594 10 E4 698 10 E5 300 30 E6 412 30 E7 472 30 E8 546 30 E9 615 30 E10 687 30 E11 733 30 E12 645 56 E13

683

60

E14 E15 E16 E17

297 549 297 530

4 7 4 4

E18

600

4

E19

670

4

IV. Experimental Conditions see Table 2

source see text measurement kASR1 ′ (0.1), kASR2 ′ (0.1) k′ASR3 (0.1), k′ASR4 (0.1)

source 86 87

measurement MOH,E1 ′ (t) (0.1)§ M′OH,E2(t) (0.1)§ M′OH,E3(t) (0.1)§ MOH,E4 ′ (t) (0.1) MOH,E5 ′ (t) (0.1) M′OH,E6(t) (0.1) MOH,E7 ′ (t) (0.1) MOH,E8 ′ (t) (0.1)§ MOH,E9 ′ (t) (0.1)§ M′OH,E10(t) (0.1) MOH,E11 ′ (t) (0.1) MHO ′ 2,E12(t) (0.2)§,b

source 26 26 26 26 present present present present present present present 50

M′HO2,E13(t) (0.2)

50

{d(ln[n-R])/dt}E14 (0.7)§ {d(ln[n-R])/dt}E15 (0.7)§ {d(ln[i-R])/dt}E16 (0.7)§ Mpropene,E17 ′ (t) (0.1)§ M′oxetane,E17(t) Macetone,E17 ′ (t) (0.3) Mpropanal,E17 ′ (t) (0.2) M′methyloxirane,E17(t) (0.1)§ M′propene,E18(t) (0.1)§ Moxetane,E18 ′ (t) (0.5)§ Macetone,E18 ′ (t) (0.2) M′propanal,E18(t) (0.2) M′methyloxirane,E18(t) (0.1)§ Mpropene,E19 ′ (t) (0.1)§ Moxetane,E19 ′ (t) (0.2)§ M′acetone,E19(t) (0.4) Mpropanal,E19 ′ (t) (0.2) Mmethyloxirane,E19 ′ (t) (0.1)§

44 44 45 61

61

61

source see text

a Note that ′ indicates ln( ) of the quantity. Uncertainties are listed in ( ). Global experiments used for optimization are denoted by §. See text for a description of the targets listed above. bHO2 measurements only before 2.5 ms are used for optimization.

Uncertainties were typically assigned estimated values or reported values, when available. The uncertainty for measured n-R logarithmic decay times44 of precursor/O2/He mixtures following pulsed photolysis to produce n-R was estimated on the basis of the scatter in the measurements. The uncertainty for measured n-R and i-R logarithmic decay times45 of precursor/ O2/He mixtures following pulsed photolysis to produce n-R and i-R was estimated on the basis of the variability in the measurements for different precursors. Uncertainties in the MPIMS measurements reflect uncertainties in the fits to the photoionization spectra and the scatter in the data. Again, weighting factors were taken to be Wi = 1/N0.5. Within the current framework, acetone and propanal are predicted to be produced largely by RO2 + RO2 reactions (ASR28−ASR31), particularly at lower temperatures, with J

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A potential contributions from i-R + O2 ↔ OH + acetone (via TS13) and n-R + O 2 ↔ OH + propanal (via TS5) only at higher temperatures. Trial optimizations that included acetone and propanal at 670 K as targets (in addition to those indicated in Table 3) required adjustments at the limits of or beyond estimated uncertainties in TS5 and T13 (3.8 kcal mol−1 lower E†TS5 and 3.0 kcal mol−1 lower E†TS13). We conducted limited electronic structure calculations to search for potentially missing low-lying transition states for concerted H atom transfer and OH loss from QOOH radicals (similar to those found in alcohol oxidation102). These explorations did not identify any such lowlying transition statesbarriers for identified transition states for QOOH-1 ↔ OH + propanal, QOOH-2 ↔ OH + propanal, and QOOH-3 ↔ OH + acetone were calculated to be 22.2, 21.1, and 23.7 kcal mol−1 above the R + O2 reactants, respectively, using M062X/6−311++G(d,p)but were not exhaustive. These calculations were performed using Gaussian 09.103 Although there are some measurements for total rate constants for n-RO2 + n-RO2 and i-RO2 + i-RO2 near 300 K, as discussed in Atkinson et al.,86 product branching ratios as well as the temperature and pressure dependence of these reactions are considerably less characterized. The RO2 + RO2 reactions proceed through a tetroxide intermediate, ROOOOR, which has two possible dissociation channels: (i) RO + RO + O2 and (ii) alcohol + aldehyde + O2.67 The difficulty in calculating temperature- and pressure-dependent rate constants for these reactions is that the second channel is spin forbidden, such that a more detailed study of the intersystem crossing is required before these reactions can be addressed quantitatively. Interpretation of acetone and propanal measurements will require better characterization of the product branching ratios as well as temperature and pressure dependence of these radical−radical reactions. Therefore, these measurements are not included as targets in the present optimization, though future modeling efforts that consider the propanal and acetone measurements would lead to a more comprehensive interpretation of available data. 4.4. New Measurements of OH Time Profiles. We carried out additional measurements of time-resolved OH concentration in Cl2/C3H8/O2/He mixtures at 30 Torr following pulsed photolysis to produce Cl, using an experimental protocol similar to that of Huang et al.26 Details of the experiment are provided in the Supporting Information.

5. RESULTS AND DISCUSSION 5.1. Target Screening and Selection for Optimization. The full set of global experiment targets listed in Table 3 was screened to determine whether predictions are affected by structural uncertainties (including static secondary reactions, RMG-identified “missing” secondary reactions, non-Boltzmann reactant distribution effects, and product excitation effects) and poorly constrained parametric uncertainties (uncertainties due to active secondary reactions for which direct rate constant measurements are not included, ASR5−ASR32). For example, the sources of variance in predictions of OH time profiles in Cl2/C3H8/O2/He mixtures at 30 Torr following pulsed photolysis to produce Cl are shown in Figure 5 for various temperatures. At low temperatures (300 K), the principal source of model variance is uncertainty in the treatment of nonBoltzmann reactant distribution effects. In particular, γn‑B, the branching ratio among products in QOOH-1* + O2 from n-R + O2, i.e., n-R + O2 (+O2), is responsible for orders of magnitude uncertainty in OH predictions at those conditions. At intermediate temperatures (615 K), the nearly exclusive source of

Figure 5. Sources of uncertainty in predictions of OH time profiles in Cl2/C3H8/O2/He mixtures at 30 Torr following pulsed photolysis to produce Cl at (a) 300 K, (b) 615 K, and (c) 733 K. The thin gray line represents the measured OH time profile from the present work.

model variance is from various theoretical kinetics parameters. At high temperatures (733 K), the primary source of model variance is from theoretical kinetics parameters, though there are substantial contributions from secondary reaction parameters. (In fact, previous interpretations of measured OH time profiles at high temperatures no longer reproduce measured OH when proper treatment of pressure dependence of n-R decomposition27 is included; see Figure 12d below.) Only targets for which predictions are unaffected by static and “missing” secondary reactions and for which uncertainties due to active secondary K

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A reactions and non-Boltzmann reaction sequences contribute less than 10% to the total model variance are included in the optimization. These targets, denoted by the § symbol in Table 3, provide the most direct information about the theoretical kinetics parameters of R1−R3 in a manner minimally affected by structural uncertainties and poorly constrained parametric uncertainties from secondary reactions. Uncertainty-weighted sensitivity analyses were performed for each of the optimization targets to characterize the information that each target provides on the model parameters, as shown in Figures 6 and 7, for example. Uncertainty-weighted sensitivity

Figure 6. Uncertainty-weighted sensitivity analysis for methyloxirane time profiles in (COCl)2/C3H8/O2/He mixtures at 4 Torr and 670 K following pulsed photolysis to produce Cl. The five largest contributors to uncertainty in methyloxirane predictions at 10 ms shown in the figure constitute 97% of the total a priori model variance. Parameters in the legend are listed in order of decreasing variance at 10 ms.

Figure 7. Uncertainty-weighted sensitivity analysis for OH time profiles in (COCl)2/C3H8/O2/He mixtures at 10 Torr and (a) 426 K and (b) 594 K following pulsed photolysis to produce Cl. The 10 largest contributors to uncertainty in peak OH shown in part a constitute 90% of the total a priori model variance; the 6 largest contributors to uncertainty in peak OH shown in part b constitute 90% of the total a priori model variance. Parameters in the legend are listed in order of decreasing variance at the time of peak OH.

analysis for methyloxirane time profiles in (COCl)2/C3H8/O2/ He mixtures at 4 Torr and 670 K following pulse photolysis to produce Cl is shown in Figure 6 using the multiscale informed model. Parameters related to the transition states for methyloxirane production (E†, v†low, and vimag for TS3) are the largest contributors to model prediction uncertainty. Because QOOH-2 quickly and exclusively decomposes to OH + methyloxirane, TS3 is rate-limiting for OH + methyloxirane production from n-RO2*. Additionally, parameters related to the transition states for HO2 + propene production (E† and v†low for TS4), which compete with TS3, are also substantial contributors to model prediction uncertainty. (Analogous channels in the i-R + O2 reaction surface are less important given the lower OH + methyloxirane yield compared to the n-R + O2 reaction surface.) TS4 parameters are among the largest contributors to model prediction uncertainty for a wide range of optimization targets, most notably propene measurements, whereas TS3 parameters are among the largest contributors to model prediction uncertainty only for methyloxirane and OH time profiles. In contrast to methyloxirane predictions, for which five parameters constitute 97% of the model variance, similar analyses for OH predictions reveal contributions from considerably more parameters, as shown in Figure 7. For example, uncertaintyweighted sensitivity analysis for OH time profiles in (COCl)2/ C3H8/O2/He mixtures at 10 Torr and 426 K following pulsed photolysis to produce Cl is shown in Figure 7a using the multiscale informed model. Parameters related to multiple transition states (e.g., E†, v†low, and vimag for TS3, TS2, and TS19) are major contributors to prediction uncertainty. As discussed above, TS3 is rate-limiting for OH + methyloxirane production

from n-RO2*. Because prompt formation of QOOH-1 exceeds its equilibrium constant with n-RO2, TS2 controls the rate at which prompt QOOH-1 is converted back to n-RO2. QOOH-1* can also react with O2 to form O2QOOH*, and TS19 controls the rate at which O2QOOH* decomposes to OH + OQ′OOH*. Similarly, uncertainties in the treatment of the non-Boltzmann reaction sequences, n-R + O2 + O2 → OH + OH + OQ′O (i.e., γn‑B), also affect predictions. Additionally, uncertainties related to partition function treatment for n-RO2 (vlow,w1), collisional energy transfer for n-RO2 (αw1), well depth for QOOH-1 (Ew3), and initial Cl concentration in the experiment (MCl,o) provide substantial contributions to the prediction uncertainty. Similar analysis for OH predictions at 594 K, shown in Figure 7b, reveals contributions from TS3 parameters (E†, v†low, and vimag), partition function and collisional energy transfer for n-RO2 (vlow,w1 and αw1), and experimental conditions (e.g., MCl,o). (Again, analogous channels in the i-R + O2 reaction surface are less important given the lower OH + methyloxirane yield compared to the n-R + O2 reaction surface.) Unlike previous modeling, where barrier heights for TS3 and TS11 of the R + O2 reaction surfaces were considered to be the predominant source of uncertainty in modeling the OH profiles,26 the present analysis shows that rigorous interpretation of the OH measurements requires consideration of many additional model aspectsincluding uncertainties in barrier L

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

consists of the active model with all parameters at their optimized values with prediction uncertainty bounds estimated through propagation using the covariance matrix. Figure 8 shows comparisons of various species time profiles measured using MPIMS in (COCl)2/C3H8/O2/He mixtures at 4 Torr at 530 and 670 K following pulsed photolysis to produce Cl,61 predictions using the rate constants derived in Huang et al.,26 and predictions using the a priori and multiscale informed models. (Figure S7 in the Supporting Information shows comparisons for 600 K, and our companion paper61 shows comparisons for all species at all temperatures.) Model predictions using the rate constants of Huang et al.26 reproduce propene reasonably well; however, they yield significantly less oxetane (by orders of magnitude) and significantly more methyloxirane (by up to a factor of 4) than observed. The a priori model predictions reproduce propene and oxetane reasonably well; however, they yield more methyloxirane (by up to a factor of 2) than observed. The multiscale informed model predictions show significant improvements over both models reproducing the measurements with reasonable consistency for propene, oxetane, and methyloxirane. A set of calculations was performed using the multiscale informed model with each parameter reset to the nominal value one at a time in order to discern the parameter adjustments responsible for the differences between the a priori and multiscale informed model predictions. The minor differences between a priori and multiscale model predictions of propene at 670 K are nearly exclusively attributable to a 15% decrease from a 15% lower [Cl]0 and an ∼5% increase from higher flux through TS4 due to a 1.3 kcal mol−1 lower E†TS4 and ∼10% lower v†low,TS4. The minor differences between a priori and multiscale model predictions of oxetane at 670 K are the result of large offsetting effects: a factor of 3 increase from a 1.5 kcal mol−1 lower E†TS6 roughly offsets a factor of nearly 3 decrease from higher flux through TS4. The factor of 2 lower methyloxirane predictions at 670 K for the multiscale informed model compared to the a priori model are the result of a slightly more than 2-fold decrease from the higher flux through TS4, an ∼40% increase from a 0.4 kcal mol−1 lower E†TS3, and an ∼15% decrease from an ∼15% lower [Cl]0. This combination of adjustments in TS3 and TS4 results in methyloxirane predictions that are a factor of 2 lower than a priori model predictions and reproduce the experimental measurements well, while also resulting in an ∼10% higher OH peak at 10 Torr and 594 K. Figure 9 shows comparisons of the measured HO2 time profiles in Cl2/C3H8/O2/He mixtures at 56 Torr and 645 K following pulsed photolysis to produce Cl,50 predictions using the rate constants derived in Huang et al.,26 and predictions using the a priori and multiscale informed models. Model predictions using the rate constants of Huang et al.,26 the a priori model, and the multiscale informed model yield somewhat more HO2 than observed, though the multiscale informed model reproduces the experimental measurements within estimated uncertainties. The differences between a priori and multiscale informed model predictions are attributable to multiple offsetting factors. The nearly 50% increase in predicted [HO2]/[Cl]0 at 5 ms from higher flux through TS4 and TS12 (E†TS4 lower by 1.3 kcal mol−1, E†TS12 lower by 0.7 kcal mol−1, and v†low,TS4 and v†low,TS12 lower by ∼5−10%) is largely offset by a corresponding decrease from differences in n- and i-RO2 properties (Ew1 and Ew4 lower by ∼0.5 kcal mol−1, vlow,w1 lower by ∼5%, and vlow,w4 lower by ∼15%), a 4 K lower Te12, and an ∼10% higher initial Cl concentration.

height, partition function treatment, and tunneling properties of multiple transition states in both R + O2 and QOOH + O2 reaction surfaces; uncertainties in secondary reactions; uncertainties in experimental conditions; and uncertainties in the treatment of non-Boltzmann reaction sequences. Therefore, the measured OH time profiles do not provide direct information regarding specific transition states responsible for OH production, but rather they provide coupled constraints on the properties of multiple transition states, wells, and other parameters. The sources of variance in predictions and uncertaintyweighted sensitivity coefficients for all global experiments are shown in Figures S1−S6 and Table S1 in the Supporting Information. The full set of optimization targets, which include measurements of OH, HO2, n-R, i-R, propene, oxetane, and methyloxirane time profiles, provide information on parameters describing numerous transition states (e.g., TS2, TS3, TS4, TS6, TS11, TS12, TS19), stable species (e.g., n/i-RO2), and experimental conditions (e.g., T, MCl,o). 5.2. Results with Multiscale Informed Model. The active parameters were optimized against and their uncertainties constrained by the targets indicated in Tables 1, 2, and 3. For the optimization, the non-Boltzmann and product-excitation factors, γn‑B and γp‑e, were set to zero with zero uncertainties. The set of optimized theoretical kinetics parameters, rate parameters, and physical model parameters from the present analysis using the targets from Tables 1−3 yield final values and predictions in reasonable consistency with ab initio calculations,59 rate constant measurements,86,87 and global measurements and experimental conditions (refs 26, 44, 45, 50, and 61 and present work) used as targets (with all |Yt,i − Fi(Xj,opt)| < Zi). In particular, as shown below, the multiscale informed model is consistent with earlier time-resolved measurements of OH at 10 Torr26 as well as the present time-resolved measurements of OH at 30 Torr and recent MPIMS measurements of methyloxirane at 4 Torr,61 for which modeling based on previous interpretations26 and recent theoretical kinetics calculations59 are inconsistent with experimental observations. A complete list of active model parameters, multiscale informed parameter values and their constrained uncertainties, and additional model predictions are provided in Table S3 and Figure S7 in the Supporting Information. Consistent with minimal effects of secondary reaction uncertainties on model predictions for the optimization targets, the optimized values of all rate constants are quite close to the nominal values (within ∼15% and generally much less) except for ASR5, ASR7, ASR11, ASR17, and ASR18 (within a factor of 2 and well within uncertainty bounds); final uncertainties in the Arrhenius parameters for ASR5−ASR32 are essentially identical to initial uncertainties. (The adjustments in secondary reaction rate constants yield minor differences in predictions of the targets. For example, a 50% higher rate constant for ASR5 yields ∼5% higher predicted OH peak concentrations, and rate constant differences for other reactions have less pronounced effects.) Below, we present comparisons of experimental measurements with model predictions from three different models. For predictions labeled “Huang et al. (2011)”, the model consists of the core model described in section 3.1 with the R1 and R2 treatment substituted by the rate constants for R1 and R2 derived in the study of Huang et al.26 using their OH measurements at 10 Torr. For predictions labeled “a priori model”, the model consists of the active model with all parameters at their nominal values. For predictions labeled “MS-informed model”, the model M

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 8. Species time profiles in (COCl)2/C3H8/O2/He mixtures at 4 Torr following pulsed photolysis to produce Cl for (a) propene, (b) oxetane, and (c) methyloxirane at 530 K and for (d) propene, (e) oxetane, and (f) methyloxirane at 670 K. The thin black line represents experimental data from Welz et al.;61 thicker lines represent model predictions. Dashed lines reflect uncertainty estimates of the multiscale informed model. The long thin dashed line in part b which represents the estimated detection limit for oxetane (rather than a measurement) is not visible because it directly overlaps the dashed line for the lower model uncertainty estimate.

Figures 10 and 11 shows comparisons of the measured logarithmic decay times for n-R and i-R in precursor/O2/He mixtures following pulsed photolysis to produce n-R and i-R44,45 and predictions using the a priori and multiscale informed models at various temperatures and pressures. The a priori and multiscale informed models yield comparable agreement with experimental data, though the multiscale informed model generally reproduces the experimental data more closelywell within a factor of 2 (the uncertainty assigned here). Differences between a priori and multiscale informed model predictions are primarily attributable to adjustments in the properties of the R + O2

reactants, RO2 radicals, collisional energy transfer, and transition states for barrierless addition. For n-R at 297 K and 4 Torr (550 K and 7 Torr), ∼30% higher f VRC‑TST,TS1 yields an ∼25% (∼15%) higher decay rate, an ∼15% higher vlow,r1 yields an ∼10% (∼30%) higher decay rate, and an ∼30% lower αw1 yields an ∼10% (∼20%) lower decay rate. Similarly, for i-R at 297 K and 4 Torr, an ∼15% lower vlow,w4 yields an ∼15% higher decay rate, an ∼15% higher f VRC‑TST,TS10 yields an ∼10% higher decay rate, and a 0.5 kcal mol−1 lower Ew4 yields an ∼5% higher decay rate. Figure 12 shows comparisons of the measured OH time profiles in (COCl)2/C3H8/O2/He mixtures at 10 Torr and N

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

same study (Huang et al.),26 and predictions using the a priori and multiscale informed models. Model predictions using the rate constants of Huang et al.26 reproduce the measurements reasonably well except at 698 K (where, as mentioned above, the originally derived rate constants no longer reproduce the measurements when pressure dependence of n-R decomposition is included). The a priori model yields OH predictions that are a factor of 2 lower than the measurements. The multiscale informed model shows significant improvementsreproducing the raw experimental data at least as well as or better than the rate constants derived originally.26 The factor of 2 higher peak OH predictions at 426 K for the multiscale informed model compared to the a priori model is the result of significant contributions from an ∼30% lower αw1, 0.5 kcal mol−1 lower E†TS3, 2.3 kcal mol−1 lower E†TS19, and 0.8 kcal mol−1 higher E†TS2. As discussed above regarding Figure 7a, a higher E†TS2 increases the fraction of prompt QOOH-1 that reacts with O2 instead of isomerizing back to n-RO2 and a lower E†TS19 increases the fraction of nascent O2QOOH* from QOOH-1 + O2 that decomposes directly to OH + OQ′OOH. A lower αw1 increases the fraction of prompt OH + methyloxirane, OH + oxetane, and QOOH-1, and a lower E†TS3 increases the production of OH + methyloxirane. Interestingly, the four adjusted parameters responsible for the vast majority of the differences between a priori and multiscale informed model predictions (αw1, E†TS3, E†TS19, E†TS2) are not simply the four greatest contributors to the model prediction uncertainty (Figure 7a) suggesting that the multiscale informed model interpretations heavily depend on the constraints imposed by other experimental and theoretical data. For example, as discussed above, the methyloxirane time profiles from the MPIMS experiments impose tight constraints on parameters describing TS3 and, therefore, restrict the potential interpretations of the OH measurements. The factor of 2 higher peak OH predictions at 594 K for the multiscale informed model compared to the a priori model are the result of comparable contributions from an ∼25% higher [Cl]0, an ∼30% lower αw1, and the combination of adjustments in TS3 and TS4 discussed above that result in methyloxirane predictions that are a factor of 2 lower than a priori model predictions and reproduce the experimental measurements well. Figure 13 shows comparisons of new measured OH time profiles in Cl2/C3H8/O2/He mixtures at 30 Torr and 300−733 K following pulsed photolysis to produce Cl, predictions using the rate constants of Huang et al.,26 and predictions using the a priori and multiscale informed models. The rate constants of Huang et al.26 reproduce the 30 Torr measurements reasonably well at intermediate temperatures; however, they predict over an order of magnitude lower OH at lower temperatures and 50% higher OH at higher temperatures than measured. Similar to the results for the 10 Torr measurements, the a priori model yields OH predictions that are significantly lower than the measurements by up to an order of magnitude at the lowest temperatures. The multiscale informed model predicts OH time profiles consistent with all experimental measurements at 10 and 30 Torr except perhaps at 30 Torr and 300 K and the highest temperatures of both data sets. At 30 Torr and 300 K, uncertainties in the treatment of non-Boltzmann reaction sequences (i.e., γn‑B and γp‑e) yield order of magnitude uncertainties in OH predictions. In fact, detailed calculations for this condition suggest non-Boltzmann reaction sequences may result in up to 3 times more OH than calculations based on traditional thermal assumptions.39,40 At the highest temperatures, secondary reaction uncertainties are

Figure 9. HO2 time profiles in Cl2/C3H8/O2/He mixtures at 56 Torr and 645 K following pulsed photolysis to produce Cl. The thin black line represents experimental data from DeSain et al.;50 thicker lines represent model predictions. Dashed lines reflect uncertainty estimates of the multiscale informed model.

Figure 10. Scaled logarithmic decay times for n-R in precursor/O2/He mixtures at 297, 400, and 550 K following pulsed photolysis to produce n-R as a function of pressure. Symbols represent experimental data from Slagle et al.44 and Ruiz and Bayes;45 thicker lines represent model predictions. Error bars reflect uncertainty estimates of the multiscale informed model.

Figure 11. Scaled logarithmic decay times for n-R and i-R in precursor/ O2/He mixtures at 298 K following pulsed photolysis to produce n-R and i-R as a function of pressure. Symbols represent experimental data from Ruiz and Bayes;45 thicker lines represent model predictions. Error bars reflect uncertainty estimates of the multiscale informed model.

426−698 K following pulsed photolysis to produce Cl,26 predictions using the rate constants originally derived in the O

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 12. OH time profiles in (COCl)2/C3H8/O2/He mixtures following pulsed photolysis to produce Cl at 10 Torr and (a) 426 K, (b) 507 K, (c) 594 K, and (d) 698 K. The thin black line represents experimental data from Huang et al.;26 thicker lines represent model predictions. Dashed lines reflect uncertainty estimates of the multiscale informed model. (Note that for predictions labeled “Huang et al. (2011)” the model consists of the rate constants for R1 and R2 derived in Huang et al.26 combined with the present core model for all other reactions including recent pressure-dependent expressions for n-R decomposition.27)

reactions (20%), pathways leading to OH + oxetane from the n-R + O2 reaction surface (7%), and pathways from the QOOH + O2 reaction surface (4%) are responsible for a substantial fraction of the total OH production. Nevertheless, the present, reasonably consistent, multiscale informed model interpretation of the available experimental and theoretical data may still not be unique. In fact, careful inspection of the results suggests that the present data set may impose tight constraints on combinations of parameters instead of individual parameterssuch that similar consistency with data may still be possible with different sets of parameters. For example, the constrained uncertainties for all barrier heights in the multiscale informed model remain above 1 kcal mol−1; constrained uncertainties for E†TS2 and E†TS19, in particular, remain above 2 kcal mol−1, and model prediction uncertainties for OH time profiles remain significantly larger (by a factor of 4) if crosscorrelations are excluded in the uncertainty propagation using the constrained model. Furthermore, interpretations of the available data appear to be quite strongly dependent on the assumptions in the treatment of non-Boltzmann reaction sequencesmodel predictions that assume γn‑B = 100% (instead of 0%) yield up to a factor of 2 more OH at conditions included as optimization targets (and up to an order of magnitude more OH at conditions excluded as optimization targets). Further studies are necessary to determine whether the parameter adjustments in

significant contributors to the model prediction variance responsible for up to 40% of the total model prediction variance at 30 Torr and 733 K. In fact, OH production from X + HO2 reactions (where X = n-R, i-R, CH3, Cl, and H) exceeds OH production from X + O2 reactions (where X = n-R, i-R, and QOOH-1) after 0.7 ms at 10 Torr and 698 K and after 0.4 ms at 30 Torr and 733 K. Consequently, proper interpretations of these experiments will require more sophisticated consideration of non-Boltzmann reaction sequences and secondary reactions. Interestingly, while the multiscale informed model reproduces the previous OH measurements as closely as the modeling of Huang et al.26 for the conditions of their experiments, it suggests a different interpretation for the exact sources of OH production (Figure 14). In the original interpretations (which were based on a more limited experimental data set),26 pathways leading to OH + methyloxirane from the R + O2 reaction surfaces were responsible for nearly all (over 90%) of the OH production at the peak at 426 and 594 K. By contrast, in the multiscale informed model interpretations, pathways leading to OH + methyloxirane from the R + O2 reaction surfaces are responsible for roughly half (54 and 65%) of OH production at 426 and 594 K. The multiscale informed model interpretation at 426 K suggests a substantial contribution from QOOH + O2 pathways responsible for 28% of OH production. Similarly, in the multiscale informed model interpretation at 594 K, R + HO2 ↔ OH + RO P

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 13. OH time profiles in Cl2/C3H8/O2/He mixtures following pulsed photolysis to produce Cl at 30 Torr and (a) 300 K, (b) 412 K, (c) 546 K, (d) 615 K, (e) 687 K, and (f) 733 K. The thin black line represents experimental data from the present work; thicker lines represent model predictions. Dashed lines reflect uncertainty estimates of the multiscale informed model.

measurements on a consistent basis to form a robust description of highly complex systems. Correspondingly, we expect the present multiscale model to provide the best current representation of low-temperature propyl oxidation kinetics. Future modeling efforts for low-temperature propane oxidation would benefit from a more complete description of active secondary chemistry, a more rigorous treatment of nonBoltzmann reaction sequences, more experimental data, and a better understanding of the inherent uncertainties in the active model parameters and cross-correlations among them. Among

the multiscale informed model are indeed physically realistic or whether the parameter adjustments suggest deficiencies in the current level of understanding for the elementary kinetics treatment of R1−R3 reaction surfaces (e.g., unidentified transition states), global kinetics model treatment (e.g., missing secondary reactions), or physical model treatment (e.g., wall reactions). While extracting direct, unambiguously interpretable information from the present data set appears very difficult, the results here demonstrate that the present approach is capable of incorporating both theoretical calculations and experimental Q

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 14. Sources of OH production at the time of the OH peak in (COCl)2/C3H8/O2/He mixtures following pulsed photolysis to produce Cl at 10 Torr and (a) 426 K and (b) 594 K using the rate constants for R1 and R2 from Huang et al.26 and the multiscale informed model.

the most valuable theoretical calculations, higher level ab initio calculations for key aspects of R1−R3 would allow determination of whether the parameter adjustments in the multiscale informed model are physically realistic; T/P/M-dependent rate constants for radical−radical reactions (e.g., RO2 + RO2, R + RO2) would enable further insight into acetone and propanal production as well as many other species (e.g., CH3, C2H4, CH2O, C2H5O2, acetaldehyde, and ROOH) from MPIMS measurements not considered here as targets; recently developed approaches to quantify and understand non-Boltzmann reaction sequences in propyl oxidation can aid in establishing a more rigorous consideration of these effects.39,40 Moreover, as a general experimental philosophy, in addition to the conventional, direct approach to studying elementary reactions through designing experiments to maximally isolate those reactions within experimental constraints, the present approach allows a different paradigm, where large amounts of experimental data are gathered across a wide range of conditions (regardless of whether they isolate particular reactions); then, optimization approaches are instead used to deconvolve complex constraints into physically meaningful parameter adjustments for multiple elementary reactions that comprise the model. Intuition would suggest that particularly valuable experiments for the present system would be measurements (e.g., both OH or HO2 laserspectroscopic measurements or multispecies MPIMS measurements) over a wider range of pressures, temperatures, initial O2 concentrations (and O2 mole fractions), initial radical concentrations, initial fuel concentrations, and initial fuel types (e.g., CH4, C2H6) that would provide a greater degree of kinetic variation (e.g., variation in the ratio of thermalizing collision time scales relative to isomerization time scales, ratio of thermalizing collision time scales relative to reactive collision time scales, extent of secondary reactions, the importance of specific secondary reactions, and relative isolation of specific subsystems) to allow discrimination among multiple interpretations and redundancy checks to ensure reliability of the model interpretations. More directed quantitative design of experiments104,105 and theoretical calculations106 to identify optimal experiments and calculations to restrict uncertainties for the low-temperature propane oxidation system are underway.

representative system for low-temperature fuel oxidation. In contrast to previous modeling, the resulting multiscale informed model provided a consistent quantitative explanation of ab initio calculations and time-resolved measurements for OH, HO2, n-propyl, i-propyl, propene, oxetane, and methyloxirane at pressures from 4 to 60 Torr and temperatures from 300 to 700 K. Previous modeling26 (based on adjustments in the R + O2 ↔ OH + methyloxirane reaction to reproduce OH measurements26) yields predictions of methyloxirane (previously thought to be the primary OH coproduct) that are up to a factor of 4 higher than MPIMS measurements. However, with the constraints imposed by the MPIMS measurements on the coproducts of OH and HO2 from the R + O2 reaction surfaces, the multiscale informed model identifies an alternative interpretation of the earlier OH measurements that is consistent with other available theoretical and experimental information, including MPIMS measurements of methyloxirane. Whereas previous interpretations26 suggested that pathways leading to OH + methyloxirane in the R + O2 reaction surfaces were responsible for nearly all (over 90%) of OH production, interpretations using the multiscale informed model suggest that QOOH + O2 reaction pathways, R + HO2 ↔ OH + RO reactions, and pathways leading to OH + oxetane from the n-R + O2 reaction surface may also constitute roughly half (up to 46%) of the OH production at experimental conditions. More generally, the present results suggest that interpretations of OH measurements even in dilute, isothermal, and isobaric alkyl oxidation experiments are significantly more complicated than previously thoughtinterpretations may depend on additional theoretical parameters for R + O2 reactions, secondary reactions, QOOH + O2 chemistry, and treatment of non-Boltzmann reaction sequences. Extraction of physically rigorous information from those measurements may require more sophisticated treatment of all of those model aspects as well as additional experimental data over a wider range of conditions to discriminate among possible interpretations and ensure model reliability.



ASSOCIATED CONTENT

S Supporting Information *

Figures S1−S6 and Table S1: Sources of variance in model predictions and uncertainty-weighted sensitivity coefficients for global experiments. Figure S7: Comparisons of experimental measurements of propene, oxetane, and methyloxirane at 600 K from Welz et al.61 and model predictions. Table S2: Nominal experimental conditions for global observable targets. Table S3: Optimized values and constrained uncertainties for multiscale informed model parameters. The Supporting Information is

6. CONCLUSIONS A major extension and generalization of the multiscale informatics approach30 was presented to treat systems of much greater complexityinvolving multiwell reactions, potentially missing reactions, nonstatistical product branching ratios, and non-Boltzmann reactant distributions. The methodology was implemented for low-temperature propane oxidation, as a R

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

(13) Qin, Z.; Lissianski, V. V.; Yang, H.; Gardiner, W. C.; Davis, S. G.; Wang, H. Combustion Chemistry of Propane: A Case Study of Detailed Reaction Mechanism Optimization. Proc. Combust. Inst. 2000, 28, 1663−1669. (14) Davis, S. G.; Joshi, A. V.; Wang, H.; Egolfopoulos, F. An Optimized Kinetic Model of H2/CO Combustion. Proc. Combust. Inst. 2005, 30, 1283−1292. (15) Frenklach, M.; Packard, A.; Seiler, P.; Feeley, R. Collaborative Data Processing in Developing Predictive Models of Complex Reaction Systems. Int. J. Chem. Kinet. 2004, 36, 57−66. (16) Sheen, D.; You, X.; Wang, H.; Løvås, T. Spectral Uncertainty Quantification, Propagation and Optimization of a Detailed Kinetic Model for Ethylene Combustion. Proc. Combust. Inst. 2009, 32, 535− 542. (17) Nagy, T.; Turányi, T. Uncertainty of Arrhenius Parameters. Int. J. Chem. Kinet. 2011, 43, 359−378. (18) Turányi, T.; Nagy, T.; Zsély, I. Gy.; Cserháti, M.; Varga, T.; Szabó, B.; Sedyó, I.; Kiss, P. T.; Zempléni, A.; Curran, H. J. Determination of Rate Parameters Based on Both Direct and Indirect Measurements. Int. J. Chem. Kinet. 2012, 44, 284−302. (19) Varga, T.; Nagy, T.; Olm, C.; Zsély, I. G.; Pálvölgyi, R.; Valkó, É; Vincze, G.; Cserháti, M.; Curran, H. J.; Turányi, T. Optimization of a Hydrogen Combustion Mechanism Using Both Direct and Indirect Measurements. Proc. Combust. Inst. 2015, 35, 589−596. (20) Sheen, D. A.; Sheen, D. A.; Rosado-Reyes, C. M.; Tsang, W. Kinetics of H Atom Attack on Unsaturated Hydrocarbons Using Spectral Uncertainty Propagation and Minimization Techniques. Proc. Combust. Inst. 2013, 34, 527−536. (21) Prager, J.; Najm, H. N.; Sargsyan, K.; Safta, C.; Pitz, W. J. Uncertainty Quantification of Reaction Mechanisms Accounting for Correlations Introduced by Rate Rules and Fitted Arrhenius Parameters. Combust. Flame 2013, 160, 1583−1593. (22) Cai, L.; Pitsch, H. Mechanism Optimization Based on Reaction Rate Rules. Combust. Flame 2014, 161, 405−415. (23) Harding, L. B.; Wagner, A. F. The Reaction of Atomic Hydrogen with the Formyl Radical. Proc. Combust. Inst. 1988, 22, 983−989. (24) Golden, D. M.; Smith, G. P.; McEwen, A. B.; Yu, C.-L.; Eiteneer, B.; Frenklach, M.; Vaghjiani, G. L.; Ravishankara, A. R.; Tully, F. P. OH(OD)+CO: Measurements and an Optimized RRKM fit. J. Phys. Chem. A 1998, 102, 8598−8606. (25) Miller, J. A.; Klippenstein, S. J.; Robertson, S. H. A Theoretical Analysis of the Reaction between Ethyl and Molecular Oxygen. Proc. Combust. Inst. 2000, 28, 1479−1486. (26) Huang, H.; Merthe, D. J.; Zádor, J.; Jusinski, L. E.; Taatjes, C. A. New Experiments and Validated Master-Equation Modeling for OH Production in Propyl + O2 Reactions. Proc. Combust. Inst. 2011, 33, 293−299. (27) Miller, J. A.; Klippenstein, S. J. Dissociation of Propyl Radicals and Other Reactions on a C3H7 Potential. J. Phys. Chem. A 2013, 117, 2718− 2727. (28) Prager, J.; Najm, H. N.; Zádor, J. Uncertainty Quantification in the ab initio Rate-Coefficient Calculation for the CH3CH(OH)CH3 + OH -> CH3C(OH)CH3 + H2O Reaction. Proc. Combust. Inst. 2013, 34, 583−590. (29) Eskola, A. J.; Carr, S. A.; Shannon, R. J.; Wang, B.; Blitz, M. A.; Pilling, M. J.; Seakins, P. W.; Robertson, S. H. Analysis of the Kinetics and Yields of OH Radical Production from the CH3OCH2 + O2 Reaction in the Temperature Range 195−650 K: An Experimental and Computational study. J. Phys. Chem. A 2014, 118, 6773−6788. (30) Burke, M. P.; Klippenstein, S. J.; Harding, L. B. A Quantitative Explanation for the Apparent Anomalous Temperature Dependence of OH + HO2 = H2O + O2 through Multi-Scale Modeling. Proc. Combust. Inst. 2013, 34, 547−555. (31) Olzmann, M.; Gebhardt, J.; Scherzer, K. An Extended Mechanism for Extended Activation Systems. Int. J. Chem. Kinet. 1991, 23, 825−835. (32) Bohn, B.; Siese, M.; Zetzsch, C. Kinetics of the OH + C2H2 Reaction in the Presence of O2. J. Chem. Soc., Faraday Trans. 1996, 92, 1459−1466.

available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.5b01003.



AUTHOR INFORMATION

Corresponding Author

*Address: 500 West 120th Street, 228 Mudd Building, MC 4703, New York, NY 10027. E-mail: [email protected]. Phone: +1(212) 851-0782. Present Address ⊥

O.W.: Institute for Combustion and Gas Dynamics, University of Duisburg-Essen, Duisburg, Germany Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by a Director’s Postdoctoral Fellowship from Argonne National Lab and by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences, under Contract No. DE-AC02-06CH11357 at Argonne as part of the ArgonneSandia Consortium on High-Pressure Combustion Chemistry, (ANL FWP # 59044; SNL FWP # 014544). Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the National Nuclear Security Administration under contract DE-AC04-94-AL85000.



REFERENCES

(1) Hwang, W.; Dec, J.; Sjöberg, M. Spectroscopic and ChemicalKinetic Analysis of the Phases of HCCI Autoignition and Combustion for Single- and Two-Stage Ignition Fuels. Combust. Flame 2008, 154, 387−409. (2) Wang, H.; Yao, M.; Reitz, R. D. Development of a Reduced Primary Reference Fuel Mechanism for Internal Combustion Engine Combustion Simulations. Energy Fuels 2013, 27, 7843−7853. (3) Reitz, R. Directions in Internal Combustion Research. Combust. Flame 2013, 160, 1−8. (4) Som, S.; Liu, W.; Zhou, D. Y.; Magnotti, G. M.; Sivaramakrishnan, R.; Longman, D. E.; Skodje, R. T.; Davis, M. J. Quantum Tunneling Affects Engine Performance. J. Phys. Chem. Lett. 2013, 4, 2021−2025. (5) Sivaramakrishnan, R.; Comandini, A.; Tranter, R. S.; Brezinsky, K.; Davis, S. G.; Wang, H. Combustion of CO/H2 Mixtures at Elevated Pressures. Proc. Combust. Inst. 2007, 31, 429−437. (6) Burke, M. P.; Chaos, M.; Dryer, F. L.; Ju, Y. Negative Pressure Dependence of Mass Burning Rates of H2/CO/O2/diluent Flames at Low Flame Temperatures. Combust. Flame 2010, 157, 618−631. (7) Burke, M. P.; Dryer, F. L.; Ju, Y. Assessment of Kinetic Modeling for Lean H2/CH4/O2/diluent Flames at High Pressures. Proc. Combust. Inst. 2011, 33, 905−912. (8) Westbrook, C. K.; Dryer, F. L. Chemical Kinetic Modeling of Hydrocarbon Combustion. Prog. Energy Combust. Sci. 1984, 10, 1−57. (9) Frenklach, M. Systematic Optimization of a Detailed KineticModel Using a Methane Ignition Example. Combust. Flame 1984, 58, 69−72. (10) Frenklach, M.; Wang, H.; Rabinowitz, M. F. Optimization and Analysis of Large Chemical Kinetic Mechanisms Using the Solution Mapping Method − Combustion of Methane. Prog. Energy Combust. Sci. 1992, 18, 47−73. (11) Eiteneer, B.; Yu, C.-L.; Goldenberg, M.; Frenklach, M. Determination of Rate Coefficients for Reactions of Formaldehyde Pyrolysis and Oxidation in the Gas Phase. J. Phys. Chem. A 1998, 102, 5196−5205. (12) Smith, G. P.; Golden, D. M.; Frenklach, M.; Moriarty, N. W.; Eiteneer, B.; Goldenberg, M.; Bowman, C. T.; Hanson, R. K.; Song, S.; Gardiner, W. C., Jr.; et al. GRI-MECH 3.0; 1999. Available at: http:// www.me.berkeley.edu/gri_mech/. S

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

and Propyl + O2 Reactions (vol 107, pg 4415, 2003). J. Phys. Chem. A 2004, 108, 7127−7128. (54) Estupinán, E. G.; Klippenstein, S. J.; Taatjes, C. A. Measurements and Modeling of HO2 Formation in the Reactions of n-C3H7 and i-C3H7 Radicals with O2. J. Phys. Chem. B 2005, 109, 8374−8387. (55) Hughes, K. J.; Griffiths, J. F.; Fairweather, M.; Tomlin, A. S. Evaluation of Models for the Low Temperature Combustion of Alkanes through Interpretation of Pressure−Temperature Ignition Diagrams. Phys. Chem. Chem. Phys. 2006, 8, 3197−3210. (56) Estupinán, E. G.; Smith, J. D.; Tezaki, A.; Jusinski, L. E.; Klippenstein, S. J.; Taatjes, C. A. Measurements and Modeling of DO2 Formation in the Reactions of C2D5 and C3D7 Radicals with O2. J. Phys. Chem. A 2007, 111, 4015−4030. (57) Gallagher, S. M.; Curran, H. J.; Metcalf, W. K.; Healy, D.; Simmie, J. M.; Bourque, G. A Rapid Compression Machine Study of the Oxidation of Propane in the Negative Temperature Coefficient Regime. Combust. Flame 2008, 153, 316−333. (58) Huynh, L. K.; Carstensen, H.-H.; Dean, A. M. Detailed Modeling of Low-Temperature Propane Oxidation: 1. The Role of the Propyl + O2 Reaction. J. Phys. Chem. A 2010, 114, 6594−6607. (59) Goldsmith, C. F.; Green, W. H.; Klippenstein, S. J. Role of O2 + QOOH in Low-Temperature Ignition of Propane. 1. Temperature and Pressure Dependent Rate Coefficients. J. Phys Chem. A 2012, 116, 3325−3346. (60) Cord, M.; Husson, B.; Lizardo Huerta, J. C.; Herbinet, O.; Glaude, P. A.; Fournet, R.; Sirjean, B.; Battin-Leclerc, F.; Ruiz-Lopez, M.; Wang, Z.; et al. Study of the Low Temperature Oxidation of Propane. J. Phys. Chem. A 2012, 116, 12214−28. (61) Welz, O.; Burke, M. P.; Antonov, I. O.; Goldsmith, C. F.; Savee, J. D.; Osborn, D. L.; Taatjes, C. A.; Klippenstein, S. J.; Sheps, L. New Insights into Low-Temperature Oxidation of Propane from Synchrotron Photoionization Mass Spectrometry and Multiscale Informatics Modeling. J. Phys. Chem. A 2015, DOI: 10.1021/acs.jpca.5b01008. (62) Curran, H. J.; Gaffuri, P.; Pitz, W. J.; Westbrook, C. K. A Comprehensive Modeling Study of n-Heptane Oxidation. Combust. Flame 1998, 114, 149−177. (63) Peters, N.; Paczko, G.; Seiser, R.; Seshadri, K. Temperature CrossOver and Non-Thermal Runaway at Two-Stage Ignition of n-Heptane. Combust. Flame 2002, 128, 38−59. (64) Kazakov, A.; Chaos, M.; Zhao, Z.; Dryer, F. L. Computational Singular Perturbation Analysis of Two-Stage Ignition of Large Hydrocarbons. J. Phys. Chem. A 2006, 110, 7003−7009. (65) Zhao, P.; Law, C. K. The Role of Global and Detailed Kinetics in the First-Stage Ignition Delay in NTC-Affected Phenomena. Combust. Flame 2013, 160, 2352−2358. (66) Beeckmann, J.; Cai, L.; Berens, A.; Peters, N.; Pitsch, H. An Analytical Approximation for Low- and High-Temperature Autoignition for Dimethyl Ether−Air Mixtures. Proc. Combust. Inst. 2015, 35, 275− 281. (67) Zádor, J.; Taatjes, C. A.; Fernandes, R. X. Kinetics of Elementary Reactions in Low-Temperature Autoignition Chemistry. Prog. Energy Combust. Sci. 2011, 37, 371−421. (68) Kennedy, M. C.; O’Hagan, A. Bayesian Calibration of Computer Models. J. R. Stat. Soc. B 2001, 63, 425−464. (69) Sargsyan, K.; Najm, H. N.; Ghanem, R. On the Statistical Calibration of Physical Models. Int. J. Chem. Kinet. 2015, 47, 246−276. (70) Ruscic, B.; Pinzon, R. E.; Morton, M. L.; von Laszevski, G.; Bittner, S. J.; Nijsure, S. G.; Amin, K. A.; Minkoff, M.; Wagner, A. F. Introduction to Active Thermochemical Tables: Several “Key” Enthalpies of Formation Revisited. J. Phys. Chem. A 2004, 108, 9979− 9997. (71) Sheen, D. A.; Wang, H. Combustion Kinetic Modeling Using Multispecies Time Histories in Shock-Tube Oxidation of Heptane. Combust. Flame 2011, 158, 645−656. (72) Davis, M. J.; Skodje, R. T.; Tomlin, A. S. Global Sensitivity Analysis of Chemical-Kinetic Reaction Mechanisms: Construction and Deconstruction of the Probability Density Function. J. Phys. Chem. A 2011, 115, 1556−1578.

(33) Knyazev, V. D.; Tsang, W. Incorporation of Non-Steady-State Unimolecular and Chemically Activated Kinetics into Complex Kinetic Schemes. 1. Isothermal Kinetics at Constant Pressure. J. Phys. Chem. A 1999, 103, 3944−3954. (34) Maranzana, A.; Barker, J. R.; Tonachini, G. Master Equation Simulations of Competing Unimolecular and Bimolecular Reactions: Application to OH Production in the Reaction of Acetyl Radical with O2. Phys. Chem. Chem. Phys. 2007, 9, 4129−4141. (35) Asatryan, R.; da Silva, G.; Bozzelli, J. W. Quantum Chemical Study of the Acrolein (CH2CHCHO) + OH + O2 Reactions. J. Phys. Chem. A 2010, 114, 8302−8311. (36) da Silva, G. Reaction of Methacrolein with the Hydroxyl Radical in Air: Incorporation of Secondary O2 Addition into the MACR + OH Master Equation. J. Phys. Chem. A 2012, 116, 5317−5324. (37) Glowacki, D. R.; Lockhart, J.; Blitz, M. A.; Klippenstein, S. J.; Pilling, M. J.; Robertson, S. H.; Seakins, P. W. Interception of Excited Vibrational Quantum States by O2 in Atmospheric Association Reactions. Science 2012, 337, 1066. (38) Pfeifle, M.; Olzmann, M. Consecutive Chemical Activation Steps in the OH-Initiated Atmospheric Degradation of Isoprene: An Analysis with Coupled Master Equations. Int. J. Chem. Kinet. 2014, 46, 231−244. (39) Burke, M. P.; Goldsmith, C. F.; Georgievskii, Y.; Klippenstein, S. J. Towards a Quantitative Understanding of the Role of Non-Boltzmann Reactant Distributions in Low-Temperature Oxidation. Proc. Combust. Inst. 2015, 35, 205−213. (40) Goldsmith, C. F.; Burke, M. P.; Georgievskii, Y.; Klippenstein, S. J. Effect of Non-Thermal Product Energy Distributions on Ketohydroperoxide Decomposition Kinetics. Proc. Combust. Inst. 2015, 35, 283−290. (41) Baker, R. R.; Baldwin, R. R.; Walker, R. W. Addition of C3H8 to Slowly Reacting Mixtures of Hydrogen and Oxygen at 480°C. Reactions of Propyl Radical. Trans. Faraday Soc. 1970, 66, 3016. (42) Baldwin, R. R.; Walker, R. W.; Yorke, D. A. Reaction of n-Propyl Radicals with Oxygen, Hydrogen, and Deuterium. J. Chem. Soc., Faraday Trans. 1 1973, 69, 826−832. (43) Baldwin, R. R.; Cleugh, C. J.; Walker, R. W. Reactions of i-Propyl Radicals with Oxygen, Hydrogen and Deuterium. J. Chem. Soc., Faraday Trans. 1 1976, 72, 1715. (44) Slagle, I. R.; Park, J.-Y.; Gutman, D. Experimental Investigation of the Kinetics and Mechanism of the Reaction of n-Propyl Radicals With Molecular Oxygen from 297 to 635 K. Proc. Combust. Inst. 1984, 20, 733−741. (45) Ruiz, R. P.; Bayes, K. D. Rates of Reaction of Propyl Radicals with Molecular Oxygen. J. Phys. Chem. 1984, 88, 2592−2595. (46) Slagle, I. R.; Ratajczak, E.; Heaven, M. C.; Gutman, D.; Wagner, A. F. Kinetics of Polyatomic Free-Radicals Produced by Laser Photolysis. 5. Study of the Equilibrium i-C3H7 + O2 ↔ i-C3H7O2 between 592 and 692 K. J. Am. Chem. Soc. 1985, 107, 1838−1845. (47) Gulati, S. K.; Walker, R. W. Arrhenius Parameters for the Reaction i-C3H7 + O2 = C3H6 + HO2. J. Chem. Soc., Faraday Trans. 2 1988, 84, 401. (48) Kaiser, E. W.; Wallington, T. J. Formation of C3H6 from the Reaction C3H7 + O2 and C2H3Cl from C2H4Cl + O2 at 297 K. J. Phys. Chem. A 1996, 100, 18770−18774. (49) Kaiser, E. W. Formation of C3H6 from the Reaction C3H7 + O2 between 450 and 550 K. J. Phys. Chem. A 1998, 102, 5903−5906. (50) DeSain, J. D.; Clifford, E. P.; Taatjes, C. A. Infrared FrequencyModulation Probing of Product Formation in Alkyl + O2 Reactions: II. The Reaction of C3H7 with O2 between 296 and 683 K. J. Phys. Chem. A 2001, 105, 3205−3213. (51) DeSain, J. D.; Taatjes, C. A.; Miller, J. A.; Klippenstein, S. J.; Hahn, D. K. Infrared Frequency-Modulation Probing of Product Formation in Alkyl + O2 Reactions - Part IV. Reactions of Propyl and Butyl Radicals with O2. Faraday Discuss. 2001, 119, 101−120. (52) DeSain, J. D.; Klippenstein, S. J.; Miller, J. A.; Taatjes, C. A. Measurements, Theory, and Modeling of OH Formation in Ethyl + O2 and Propyl + O2 Reactions. J. Phys. Chem. A 2003, 107, 4415−4427. (53) DeSain, J. D.; Klippenstein, S. J.; Miller, J. A.; Taatjes, C. A. Measurements, Theory, and Modeling of OH Formation in Ethyl + O2 T

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (73) Miki, K.; Prudencio, E. E.; Cheung, S. H.; Terejanu, G. Using Bayesian Analysis to Quantify Uncertainties in the H + O2 = OH + O Reaction. Combust. Flame 2013, 160, 861−869. (74) Tomlin, A. S. The Role of Sensitivity and Uncertainty Analysis in Combustion Modelling. Proc. Combust. Inst. 2013, 34, 159−176. (75) Zhou, D. D. Y.; Davis, M. J.; Skodje, R. T. Multitarget Global Sensitivity Analysis of n-Butanol Combustion. J. Phys. Chem. A 2013, 117, 3569−3584. (76) Tomlin A. S.; Turanyi, T. Investigation and Improvement of Reaction Mechanisms Using Sensitivity Analysis and Optimization. In Cleaner Combustion; Battin-Leclerc, F., Simmie, J. M., Blurock, E., Eds.; Springer: London, 2013; pp 411−445. (77) Braman, K.; Oliver, T. A.; Raman, V. Bayesian Analysis of Syngas Chemistry Models. Combust. Theor. Modell. 2013, 17, 858−887. (78) Mosbach, S.; Hong, J. H.; Brownbridge, G.; Kraft, M.; Gudiyella, S.; Brezinsky, K. Bayesian Error Propagation for a Kinetic Model of nPropylbenzene Oxidation in a Shock Tube. Int. J. Chem. Kinet. 2014, 46, 389−404. (79) Edwards, D. E.; Zubarev, D. Yu.; Packard, A.; Lester, W. A., Jr.; Frenklach, M. Interval Prediction of Molecular Properties in Parametrized Quantum Chemistry. Phys. Rev. Lett. 2014, 112, 253003. (80) Wang, H.; Sheen, D. A. Combustion Kinetic Model Uncertainty Quantification, Propagation and Minimization. Prog. Energy Combust. Sci. 2015, 47, 1−31. (81) Baulch, D. L.; Bowman, C. T.; Cobos, C. J.; Cox, R. A.; Just, T.; Kerr, J. A.; Pilling, M. J.; Stocker, D.; Troe, J.; Tsang, W.; et al. Evaluated Kinetic Data for Combustion Modeling: Supplement II. J. Phys. Chem. Ref. Data 2005, 34, 757−1397. (82) Green, W. H.; Allen, J. W.; Buesser, B. A.; Ashcraft, R. W.; Beran, G. J.; Class, C. A.; Gao, C.; Goldsmith, C. F.; Harper, M. R.; Jalan, A.; et al. RMG - Reaction Mechanism Generator v4.0.1; 2013; http://rmg. sourceforge.net/. (83) Burke, M. P.; Chaos, M.; Ju, Y.; Dryer, F. L.; Klippenstein, S. J. Comprehensive H2/O2 Kinetic Model for High-Pressure Combustion. Int. J. Chem. Kinet. 2012, 44, 444−474. (84) Miller, J. A.; Burke, M. P.; Glarborg, P.; Goldsmith, C. F.; Hansen, N.; Jasper, A. W.; Sivaramakrishnan, R.; Zádor, J.; Klippenstein, S. J. A Comprehensive Kinetic Model for C0-C3 Combustion. As updated Dec 2012. (85) Dooley, S.; Burke, M. P.; Chaos, M.; Stein, Y.; Dryer, F. L.; Zhukov, V. P.; Finch, O.; Simme, J. M.; Curran, H. J. Methyl Formate Oxidation: Speciation Data, Laminar Burning Velocities, Ignition Delay Times, and a Validated Chemical Kinetic Model. Int. J. Chem Kinet. 2010, 42, 527−549. (86) Atkinson, R.; Baulch, D. L.; Cox, R. A.; Crowley, J. N.; Hampson, R. F.; Hynes, R. G.; Jenkin, M. E.; Rossi, M.; Troe, J. Evaluated Kinetic and Photochemical Data for Atmospheric Chemistry: Volume II - Gas Phase Reactions of Organic Species. J. Atmos. Chem. Phys. 2006, 6, 3625−4055. (87) Droege, A. T.; Tully, F. P. Hydrogen-Atom Abstraction from Alkanes by OH: 3. Propane. J. Phys. Chem. 1986, 90, 1949−1954. (88) Jasper, A. W. Personal communication. (89) Zádor, J.; Miller, J. A. Unimolecular Dissociation of Hydroxypropyl and Propoxy Radicals. Proc. Combust. Inst. 2013, 34, 519−526. (90) Lee, J.; Bozzelli, J. W. Thermochemical and Kinetic Analysis of the Formyl Methyl Radical + O2 Reaction System. J. Phys. Chem. A 2003, 107, 3778−3791. (91) Goldsmith, C. F.; Magoon, G. R.; Green, W. H. Database of Small Molecule Thermochemistry for Combustion. J. Phys. Chem. A 2012, 116, 9033−9057. (92) Ruscic, B. Active Thermochemical Tables: Version Alpha 1.110 of the Core (Argonne) Thermochemical Network, release date April 2, 2011. (93) Goos, E.; Burcat, A.; Ruscic, B. Extended Third Millennium Ideal Gas and Condensed Phase Thermochemical Database for Combustion with Updates from Active Thermochemical Tables; ftp://ftp.technion. ac.il/pub/supported/aetdd/thermodynamics; mirrored at http:// garfield.chem.elte.hu/Burcat/burcat.html, accessed Feb 8, 2012. (94) Goos, E. Personal communication.

(95) Georgievskii, Y.; Miller, J. A.; Burke, M. P.; Klippenstein, S. J. Reformulation and Solution of the Master Equation for Multiple Well Reactions. J. Phys. Chem. A 2013, 117, 12146−12154. (96) Georgievskii, Y.; Jasper, A. W.; Zádor, J.; Miller, J. A.; Burke, M. P.; Goldsmith, C. F.; Klippenstein, S. J. PAPER, v1; 2015. (97) Klippenstein, S. J.; Wagner, A. F.; Dunbar, R. C.; Wardlaw, D. M.; Robertson, S. H.; Miller, J. A. VARIFLEX, version 1.14m; 2005. (98) Truhlar, D. G.; Garrett, B. C.; Klippenstein, S. J. Current Status of Transition-State Theory. J. Phys. Chem. A 1996, 100, 12771−12800. (99) Miller, J. A.; Klippenstein, S. J. Master Equation Methods in Gas Phase Chemical Kinetics. J. Phys. Chem. A 2006, 110, 10528−10544. (100) Miller, J. A.; Klippenstein, S. J. Determining Phenomenological Rate coefficients from a Time-Dependent, Multiple-Well Master Equation: ‘‘Species Reduction’’ at High Temperatures. Phys. Chem. Chem. Phys. 2013, 15, 4744−4753. (101) Goldsmith, C. F.; Tomlin, A. S.; Klippenstein, S. J. Uncertainty Propagation in the Derivation of Phenomenological Rate Coefficients from Theory: A Case Study of n-Propyl Radical Oxidation. Proc. Combust. Inst. 2013, 34, 177−185. (102) Welz, O.; Klippenstein, S. J.; Harding, L. B.; Taatjes, C. A.; Zádor, J. Unconventional Peroxy Chemistry in Alcohol Oxidation: The Water Elimination Pathway. J. Phys. Chem. Lett. 2013, 4, 350−354. (103) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; et al. Gaussian 09, revision D.01; Gaussian, Inc.: Wallingford, CT, 2009. (104) Scire, J. J., Jr.; Yetter, R. A.; Dryer, F. L. Flow Reactor Studies of Methyl Radical Oxidation Reactions in Methane-Perturbed Moist Carbon Monoxide Oxidation at High Pressure with Model Sensitivity Analysis. Int. J. Chem. Kinet. 2001, 33, 75−100. (105) Sheen, D. A.; Manion, J. Kinetics of the Reactions of H and CH3 Radicals with n-Butane: An Experimental Design Study Using Reaction Network Analysis. J. Phys. Chem. A 2014, 118, 4929−4941. (106) Skodje, R. T.; Tomlin, A. S.; Klippenstein, S. J.; Harding, L. B.; Davis, M. J. Theoretical Validation of Chemical Kinetic Mechanisms: Combustion of Methanol. J. Phys. Chem. A 2010, 114, 8286−8301.

U

DOI: 10.1021/acs.jpca.5b01003 J. Phys. Chem. A XXXX, XXX, XXX−XXX