Detailed reaction models for C 1 -C 4 hydrocarbons usually contain over 100 species in about 10 0 0 elementary reactions [1][2][3][4][5][6] . Direct application of such models is limited only to simple, canonical combustion simulations because of their tremendous computational cost. Various reduction techniques have been developed to accommodate realistic fuel chemistry simulations, and to capture intricacies of chemical kinetics in complex multidimensional combustion systems. As the first step in developing model reduction, it is important to extract a subset of the detailed reaction model, skeletal model , by eliminating unimportant species and reactions [7,8] . Local sensitivity analysis (SA), reaction flux analysis [9][10][11] , and directed relation graph (DRG) and its variants [12][13][14][15] have often been utilized for skeletal model reduction.
Local SA, which is the subject of the present work, explores the response of model output to a small change of a parameter from its nominal value [16] while global sensitivity analysis is useful for studying uncertainty of kinetic parameters ( i.e. collision frequencies and activation energies) which propagate through model and non-linear coupling effects [2,[17][18][19][20][21][22][23][24] .
Model reduction with local SA contains methods such as PCA [1,[25][26][27][28][29][30][31] , and construction of a species ranking [32] . In local SA, the sensitivities are commonly computed either by finite difference (FD) discretizations, directly solving a sensitivity equation (SE), or by an adjoint equation (AE) [33] . The computational cost of using FD or SE, which are forward in time methods, scales linearly with the number of parameters making them impracticable when sensitivities with respect to a large number of parameters are needed. On the other hand, computing sensitivities with AE requires a forward-backward workflow, but the computational cost is independent of the number of parameters as it requires solving a single ordinary/partial differential equation (ODE/PDE) [34][35][36] . The AE solution is tied to the objective function, and for cases where multiple objective functions are of interest, the same number of AEs must be solved. Regardless of the method of computing sensitivities, the output of FD, SE, and AE at each time instance is the full sensitivity coefficient matrix, which can be extremely large for systems with large number of parameters.
Recently, the forced optimally time dependent (f-OTD) decomposition method was introduced for computing sensitivities in evolutionary systems using a model driven low-rank approximation [33] . This methodology is the extension of OTD decomposition in which a mathematical framework is laid out for the extraction of the low-rank subspace associated with transient instability of the dynamical system [37] . The OTD approximates sensitivities with respect to initial conditions, while f-OTD approximates sensitivities with respect to external parameters, e.g. , forcing. As a consequence, in the formulation of f-OTD there is a two-way coupling between the evolution of the f-OTD modes and the f-OTD coefficients, whereas in OTD formulation, the evolution of the modes is independent of the coefficients. In forward workflow of f-OTD, the sensitivity matrix i.e. S(t ) ∈ R n eq ×n r is modeled on-the-fly as the multiplication of two skinny matrices
×r , and
×r which contain the f-OTD modes and f-OTD coefficients, respectively, where n eq is the number of equations (or outputs), n r is the number of independent parameters, r min{ n s , n r } is the reduction size, and
The key characteristic of f-OTD is that both U (t) and Y (t) are time-dependent and they evolve based on closed form evolution equations extracted from the model, and are able to capture sudden transitions associated with the largest finite time Lyapunov exponents [38] . The time-dependent bases have also been used for stochastic reduced order modeling [39][40][41][42][43] , and recently for on-the-fly reduced order modeling of reactive species transport equation [44] . In a nutshell, f-OTD workflow i) is forward in time unlike AE, ii) bypasses the computational cost of solving FD and SE, or other data-driven reduction techniques, and iii) stores the modeled sensitivities in a compressed format, i.e. we only store and solve for two skinny matrices U and Y instead of storing and solving the full sensitivity matrix S, as in FD, SE and AE.
The major advantage of PCA in skeletal reduction is to combine the sensitivity coefficients for a wide range of operating conditions ( e.g. equivalence ratio and pressure) [1] . The PCA finds the low-dimensional subspace of data gathered from different (temporal or spatial) locations by applying a minimization algorithm over the whole data at once. Therefore, PCA is a low-rank approximation in a time-averaged sense and may fail to capture highly transient finite-time events ( e.g. ignition). In order to resolve this issue, one needs to pre-recognize the locations of such events and use the data mainly from these locations. This requires extensive knowledge and/or expertise. Moreover, PCA modes are time invariant, and the process of selecting sufficient eigenvalues/eigenmodes to capture the essence of all observed phenomena ( e.g. ignition, flame propagation), is crucial but is usually done by trial and error. References [1,45] show that for certain problems, a skeletal model built solely upon the information conveyed by that first reaction group (first eigenmode) from PCA can fail to accurately reproduce the detailed model over the entire domain of interest. Therefore, one needs to deal with several eigenmodes with close eigenvalues and choose essential reaction groups among them [1] .
In order to resolve the drawbacks of current SA methods and PCA for skeletal model reduction, we use f-OTD methodology for both SA and skeletal model reduction. The applicability of our approach is demonstrated for ethylene-air burning with the University of Southern California (USC) chemistry model [2] as the detailed model. Adiabatic, constant pressure, spatially homogeneous ignition is the canonical problem; and the generated skeletal models with f-OTD are compared against detailed and several skeletal models.
The remainder of this paper is organized as follows. The theoretical description of PCA and f-OTD and their mathematical derivations for SA are presented in Section 2 . Model reduction with f-OTD is first described in Section 3 , with a simple reaction model for hydrogen-oxygen combustion, followed by skeletal model reduction with f-OTD for the more complex ethylene-air system in Section 4 . The paper ends with conclusions in Section 5 . All the generated models are supplied in supplementary materials section.
Consider a chemical system of n s species reacting through n r irreversible reactions,
where M k is a symbol for species k , and ν k j and ν k j are the molar stoichiometric coefficients of species k in reaction j. Changes of mass fractions ψ = [ ψ 1 , ψ 2 , . . . , ψ n s ] T and temperature T in an adiabatic, constant pressure p, and spatially homogeneous reaction system of ideal gases can be described by the following initial value problems (IVPs) [46] dψ k dt
where t ∈ [0 , t f ] is time, t f is the final time and W k and h k are the molecular weight and enthalpy of species k , respectively, and
Here, α = [1 , 1 , . . . , 1] ∈ R n r is the vector of perturbation parameters and K j is the rate constant of reaction j which is usually modeled using the modified Arrhenius parameters [47] for elementary reactions (Note: all reversible reactions are cast as irreversible reactions). In Eq. ( 2) ρ(T , ψ ) and c p (T , ψ ) = n s k =1 ψ k c p k (T ) are the density and specific heat at constant pressure of the mixture, respectively, where c p k (T ) is the specific heat at constant pressure of k th species given by the NASA coefficient polynomial parameterization [48] . Let ξ = [ ψ, T ] ∈ R n eq denote the vector of compositions and accordingly f = [ f ψ , f T ] where n eq = n s + 1 . Then the compositions IVP is governed by:
Since α = 1 , the perturbation with respect to α j amounts to an infinitesimal perturbation of progress rates Q j for j = 1 , 2 , . . . , n r .
The sensitivity matrix, S(t ) = [ s 1 (t) , s 2 (t ) , . . . s n r (t )] ∈ R n eq ×n r , contains local sensitivity coefficients, s j = ∂ ξ/∂ α j , and it can be calculated by solving the SE,
where
are the Jacobian and the forcing matrices, respectively.
where p 0 and p are unperturbed and perturbed normalized kinetic parameters, respectively; and p j = ln α j for j = 1 . . . n r . The integrated squared deviation is investigated on the interval [ z 1 , z 2 ] of the independent variable (time and/or space) [45] . It has been
shown [25] that G(p) can be approximated around the nominal parameter set ( p 0 ) as,
where p = pp 0 and
Here,
) on a series of m quadrature points on [ z 1 , z 2 ] to approximate the integral in Eq. ( 6) . Eigen decomposition of S T S = A A T results in:
where
Like PCA, our kinetic model reduction strategy is based on selecting reactions, whose perturbations grow most intensely in the composition evolution given by Eq. ( 4) . However, the selection of important reactions is made here based on instantaneous observation of modeled sensitivities, unlike PCA.
Imagine we perturb the composition evolution equation ( Eq. ( 4) ) by infinitesimal variations of α j = 1 to α j = 1 + δα j , where δα j 1 for j = 1 , 2 , . . . , n r . In f-OTD, we factorize the sensitivity matrix S(t ) into a time-dependent subspace in the n eqdimensional phase space of compositions represented by a set of f-OTD modes:
where δ i j is the Kronecker delta. The rank of S(t ) ∈ R n eq ×n r is d = min { n eq , n r } while the f-OTD modes represent a rank-r subspace, where r d. To this end, we approximate the sensitivity matrix via the f-OTD decomposition ( Fig. 1 ):
where
×r is the f-OTD coefficient matrix. The above decomposition is not exact as Eq. ( 10) is a low-rank approximation of the sensitivity matrix S(t ) . Note that in the above decomposition both U (t) and Y (t) are time dependent.
We drop the explicit time dependence on t for brevity. Fig. 1 shows the schematic of decomposition of S into f-OTD components U and Y . The evolution equation for U and Y are obtained by substituting Eq. (10) into Eq. ( 5) :
Projecting Eq. ( 11) to U results in
The f-OTD modes are orthonormal, U T U = I. Taking a time derivative of the orthonormality condition yields in:
As shown in Refs. [33,37] , any skew-symmetric choice of matrix ϕ, will lead to equivalent f-OTD subspaces. Here we choose ϕ = 0 . Using U T U = I and U T dU dt = 0 , Eq. ( 12) simplifies to
The evolution equation for U can be obtained by substituting dY T dt from Eq. ( 13) in Eq. ( 11) and projecting the resulting equation onto Y by multiplying Y from right
where Q = I -U U T is the orthogonal projection onto the space spanned by the complement of U and C = Y T Y ∈ R r×r is a correlation matrix matrix. Matrix C(t ) is, in general, a full matrix implying that the f-OTD coefficients are correlated. Eq. ( 13) can be written as
where L r = U T LU ∈ R r×r is a reduced linearized operator. Eqs. ( 14) and ( 15) are a coupled system of ODEs and they constitute the f-OTD evolution equations. The f-OTD modes align themselves with the most instantaneously sensitive directions of the composition evolution equation when perturbed by α. It is shown in Ref. [38] that when α is the perturbation to the initial condition, the OTD modes converge exponentially to the eigen-directions of the CauchyGreen tensor associated with the most intense finite-time instabilities. We refer to Ref. [37] for further details about OTD.
In our approach, instantaneous sensitivities are analyzed as opposed to PCA in which sensitivities are sampled at a few selected time instants. Computing instantaneous sensitivities would significantly increase the computational/bookkeeping costs especially if large detailed mechanisms are to be analyzed. The f-OTD leverages the fact that we are often interested in the leading sensitivity vectors, and it provides a computationally feasible approach for solving only those dominant sensitivity vectors, without requiring to explicitly compute the full sensitivities at any point. The reduction procedure is as follows:
1. Modeled sensitivities are computed in factorized format by solving Eqs. ( 14) and (15) . These two equations are evolved in addition to Eq. ( 4) , and the values of ξ , U, and Y are stored at resolved time steps t i ∈ [0 , t f ] . Eq. ( 4) is initialized with a combination of the initial temperatures, equivalence ratios ( φ 0 ), etc .
Each simulation with a different initial condition is denoted by a case here. Eqs. ( 14) and ( 15) are initialized by first solving the SE ( Eq. ( 5) ) for a few time steps, and then performing singular value decomposition of the sensitivity matrix S = B V T and assigning U = B and Y = V . 2. At each resolved time step and for each case, we compute the eigen decomposition of
The w vectors are basically the average of eigenvectors of ≈ S T ≈ S matrix weighted based on their associated eigenvalues. This prevents us from dealing with each eigenvector ( a i ) separately. We use the normalized sensitivities (
) of species with mass fractions greater than a threshold, e.g. 1.0e-6. This results in elimination of some of the rows of ≈ S associated with small mass fraction species at each time step. As shown in Sections 3 and 4 , the first sorted eigenvalue ( λ 1 ) is usually orders of magnitude larger than the The larger the w i value, the more important reaction i is. We define χ i as the highest value associated with w i through all resolved time-steps and cases.
3. The elements of χ vector are sorted in descending order to find the indices of the most important reactions in the detailed model. Species are also sorted based on their first presence in the sorted reactions, i.e. species who first shows up in a higher ranked reaction would be more important than a species who first participate in a lower ranked reaction. This results in a reaction and species ranking based on χ vector e.g. Fig. 2 (b).
4. In the last step, we choose a set of species by defining a threshold χ on the element of χ vector and eliminate species whose associated χ i are less than χ . We also get rid of the reactions which include the species. As our skeletal model reduction is reaction based, any non-reactive species with non-zero mass fraction in the initial condition should be manually added to the skeletal model.
In summary, we sort the reactions to find the most important species. These species and the reactions which connect them together, create our reduced models. In combustion systems, perturbations with respect to "fast'' reactions generate very large sensitivities for short time periods which vanish as t → ∞ . On the other hand, perturbations with respect to "slow'' reactions generate smaller and more sustained sensitivities. As our approach is based on the instantaneous observation of sensitivities, both "slow'' and "fast'' reactions can leave an imprint on the instantaneous normalized reaction vector ( w ) if their imprints are larger than the threshold value (χ ) . However, if the sensitivities asso- ciated with "fast'' and "slow'' reactions from all times and locations would be combined with each other before dimension reduction, as commonly done in PCA-type approaches, the smaller sensitivities associated with "slow'' reactions could have been outweighted by the large sensitivities associated with "fast'' reactions. In fact, this concern is the primary motivation for using f-OTD rather than PCA-type skeletal reduction approaches.
In this section, the process of eliminating unimportant reactions and species from a detailed kinetic model with f-OTD is described, and its differences with PCA are highlighted. The Burke model [49] for hydrogen-oxygen system which contains n s = 10 species, 1 and n r = 54 irreversible (27 reversible) reactions is considered as the detailed model. The reduction process is performed by analyzing only one case for the ignition of an adiabatic, stoichiometric hydrogen-oxygen mixture at atmospheric pressure and T 0 = 1200K, with integration of both the SE ( Eq. ( 5) ) and f-OTD equations ( Eqs. ( 14) and ( 15) ). Exact sensitivities from SE are computed for two purposes, i) finding PCA eigenmodes and eigenvalues, and ii) analyzing the performance of f-OTD by comparing the instant eigenvalues of ≈ S T ≈ S at each t, from f-OTD against those obtained by solving the SE. The latter is equivalent to performing in-stantaneous PCA (I-PCA) on the full sensitivity matrix. The I-PCA shows the optimal reduction of the time-dependent sensitivity matrix, and we show that the eigenvalues of f-OTD closely approximate the r most dominant eigenvalues of I-PCA.
Figure 2 (a) compares top eigenvalues of f-OTD with PCA (static) and I-PCA (instantaneous). It is shown that the top PCA eigenvalues are time invariant, and close to each other. In contrast the first f-OTD eigenvalue is orders of magnitude larger than the others during the course of ignition, i.e. from t= 0 to t= 30 µs (until most of the heat is released). Moreover, f-OTD eigenvalues match with I-PCA with increasing number of modes ( r). This means that the modeled sensitivities converge to the exact values by adding more modes, in this case addition of top 6 modes. The results also show that with f-OTD ( r= 5), the time variation of top eigenvalues is captured well, while the second dominant eigenvalues deviates from I-PCA solution in the main non-equilibrium reaction layer and the post heat release region. ample the first reaction (H+O 2 → O+OH) is the most important one at t= 5.0e-6, but by marching in time and passing the peak of the heat release region, other reactions ( e.g. reactions 2-10 which are specifically radical recombination reactions like OH+OH → O+H 2 O) also become important. Our model reduction approach ranks reactions based on their all time maximum value on w. The f-OTD approach detects reactions even if their importance become visible instantly. Such reactions cannot be detected in static model reduction techniques such as PCA.
Figure 5 shows the first 8 eigenmodes of S T S matrix ( Eq. ( 8) )
and
Several detailed kinetic models for ethylene-air burning are available in literature, and are developed at the University of California, San Diego (UCSD) [4] , the University of Southern California (USC -a subset of JetSurf) [2] , the KAUST (AramcoMech2) [50] , and the Politecnico of Milan (CRECK) [5] . Figure 6 indicates that the ignition delays as predicted by all these models are in a reasonable agreement with each other. Moreover, it is shown in Ref. [51] that USC ignition delays are closer to experimental data in comparison with the other three models. Therefore, here we consider USC as our detailed kinetic model, and extract a series of skeletal models via comparison with this baseline.
Simulations are conducted of an adiabatic, atmospheric pressure ( p= 1atm) reactor for 6 cases with different initial temperatures T 0 ∈ [140 0,20 0 0] and equivalence ratios φ ∈ [0.5,1.0,1.5] for ethylene-air mixture. The USC model [2] with 111 species and 1566 irreversible (784 reversible) reactions is the detailed model based on which all the skeletal models (f-OTDs) are generated. Only three f-OTD modes ( r= 3) are employed to model the sensitivity matrix. Simulations with SK31 [1] , SK32 [52] and SK38 [1] , which are also skeletal models generated from two versions of USC (optimized and unoptimized), are considered. The comparisons are made based on three criteria: i) ignition delay, ii) premixed laminar flame speed, iii) non-premixed extinction strain. The flame speeds and the extinction curves are generated by Cantera [48] . S for one of the cases with T 0 = 1400K and φ 0 = 1.0. The value of λ 1 (t) is two orders of magnitude larger than λ 2 (t) during the course of ignition, except around temperature inflection pointfoot_2 ( T in f ) in which λ 1 (t) is six orders of magnitude larger than λ 2 (t) . This means that S is dominant during the ignition phenomenon and contains more than 95% of the energy of the dynamical system ( λ 1 (t) / i λ i (t) > 0 . 95 ). Figure 7 (b) shows the species ranking based on the process described in Section 2.2.2 . Similar to the species ranking of hydrogen-oxygen system presented in Section 3 , H, OH, O and O 2 appear as the most important species based on their associated values on χ vector for the most sensitive reaction, which is H+O 2 → OH+O. Different skeletal models can be generated by putting different threshold χ on χ and eliminating species with χ < χ and their associated reactions. Our goal is to find a model which can reproduce the results of USC model based on the criteria mentioned in Section 4.1 with a pre-determined accuracy, e.g. less than 5% error. Table 1 provides the details of models generated with varying threshold values of χ .
Figures 8 (a) and 9 (a) demonstrate that f-OTD models with n s ≥ 32 perfectly predict the ignition delays for the stoichiometric mixture. The SK32 model under-predicts the ignition delays while SK38 and SK31 (with slightly different rate constants from USC) over-predict the ignition delays, and the relative error associated 1 with T 0 = 1400K, φ= 1.0, and at 1atm pressure. f-OTD models with n s ≥ 32 show strong ability in reproducing USC results.
with these three models are usually higher than f-OTD models with n s ≥ 32 . Figure 8 (b) and 9 (b) compare the laminar premixed flame speeds as predicted by different models initialized with T 0 = 300K. The f-OTD-2, SK31 and SK32 show largest deviations from the USC model. The f-OTD models with n s ≥ 38 and SK38 have the best flame speed predictions, while f-OTD models show better comparisons at lower and upper bounds of φ 0 . Figure 8 (c) compares the extinction strain rates. All f-OTD models show good agreements in estimating these rates. This is the toughest canonical flame feature to predict. The SK38 model under-predicts the maximum temperatures indicating the influence of optimized rate constants in Ref. [2] .
Figure 10 portrays the species mass fraction evolution for some key species in a mixture initialized with φ 0 = 1.0 and T 0 = 1400K.
This figure highlights the ability of f-OTD-2 (with 32 species) in predicting ignition. Moreover, all f-OTD models (with n s ≥ 32 ) provide a better estimate for the maximum mass fraction of species shown in Fig. 10 as compared with SK32 model. Observing the results in Fig. 8 , it is clear that f-OTD-1 is not a good skeletal model for USC. This is attributed to the elimination of C 2 O, CH 2 OCH 2 , CH 3 O, and H 2 O 2 in this 28 species model. Although f-OTD-1 cannot predict the ignition delay accurately, it performs reasonably well in estimating laminar flame speeds and maximum temperatures for extinction. As mentioned above, f-OTD-2 estimates the ignition delays and extinction strain rate reasonably. Moreover, its flame speed predictions also match those via SK32 model. We recommend this model when 10 % of relative error is tolerable in predicting ignition delays and laminar flame speeds.
Predictions with f-OTD-3, f-OTD-4, and f-OTD-5 for all the test cases are so close to USC, and become more precise by increasing n s . Comparing f-OTD-3 with f-OTD-5 based on their applications in reproducing USC model and their computational costs, we recommend the former. As Fig. 9 shows, f-OTD-3 predicts the results of USC with less than 5 % relative error. Here are some differences/similarities between participating species in f-OTD models and SK32 and SK38: • f-OTD-3 has 30 species in common with SK32.
Figure 11 demonstrates the ability of f-OTD-3 in predicting ignition delay, laminar flame speed and extinction strains for three different pressures: 0.5atm, 1.0atm, and 3.0atm. f-OTD-3 shows strong ability in reproducing USC results. In this regard, it should be mentioned that the USC model contains certain reactions which are tuned for atmospheric pressure conditions. Thus, any large excursions from 1atm, e.g. 0.1 or 10atm, require manual selection of alternate set of rate parameters; the process we have tried to avoid here. We believe the range of pressure, i.e. 0.5 to 3atm (pressure ratio of 6) provides a decent test of fall-off pressure effects in our results.
As mentioned earlier, f-OTD analyzes the eigen decompostion of ≈ S T ≈ S instantly instead of S T S in PCA ( Eq. 8 ). Using normalized sensitivities (
) can produce numerical issues where mass fractions ( ξ i ) approach zero. In this study we used the normalized sensitivities of species with mass fractions greater than a threshold, e.g. 1.0e-6. This results in elimination of some of the rows of ≈ S at each time step.
Another approach is to analyze the eigen decompostion of S T S instead of ≈ S T ≈ S in f-OTD. This approach does not need a threshold for species mass fractions, and results in almost exactly the same skeletal models (f-OTD * s models). f-OTD-1 * , f-OTD-3 * are exactly the same as f-OTD-1 and f-OTD-3, while f-OTD-2 * , f-OTD-4 * , and f-OTD-5 * have only one species difference with f-OTD-2, f-OTD-4, and f-OTD-5, respectively. f-OTD-2 * uses C 4 H 2 but f-OTD-2 uses C 2 O instead, f-OTD-4 * uses H 2 C 4 O but f-OTD-4 uses aC 3 H 4 instead, and f-OTD-5 * uses CH 2 OH but f-OTD-5 uses C 2 H 3 CHO instead. Moreover, the test results do not show significant differences for this one species difference.
Instantaneous sensitivity analysis with f-OTD is described and implemented for a systematic skeletal model reduction. A key feature of the f-OTD approach is that it factorizes the sensitivity matrix into a multiplication of two low-ranked time-dependent matrices which evolve based on evolution equations derived from the governing equations of the system. Modeled sensitivities are then normalized and the most important reactions and species of a detailed model are ranked in a systematic manner based on the correlations between normalized sensitivities. It is also shown that analyzing the correlations between the non-normalized sensitivities also result in almost the same skeletal models. The significance of the f-OTD approach in model reduction is described for hydrogen-oxygen combustion, and its application for skeletal reduction is demonstrated for ethylene-air burning. The generated skeletal models are compared based on their ability to predict ignition delays, flame speeds and diffusion flame extinction strain rates. f-OTD demonstrates strong ability in eliminating unimportant species and reactions from the detailed model in an efficient manner. We recommend using f-OTD-2 and f-OTD-3 as skeletal models for USC with 10 % and 5 % relative errors, respectively, in estimating USC ignition delays and laminar flame speeds.
The extension of this study would include sensitivity analysis based on the most effective thermochemistry parameters e.g. activation energies, formation enthalpies, and transport properties e.g. heat and mass diffusivities. Most importantly, as shown recently [33] , f-OTD can be used in solving PDEs for multi-dimensional combustion problems in a cost-effective manner -by exploiting the correlations between the spatiotemporal sensitivities of different species with respect to different parameters. This analysis can be especially insightful for problems containing rare events e.g. deflagration-detonation-transition by providing more insight about the global effective phenomena. Moreover, the f-OTD provides an excellent setting for development of reduced schemes in other mechanisms, for example thermonuclear reactions [53] .
This kinetic model has 13 species (H, H
, O, OH, H 2 O, O 2 , HO 2 , H 2 O 2 , N 2 , AR, HE, CO, CO 2 ) in which N 2 , CO, and CO 2 do not participate in the reactions.
Temperature at d T /d t| max