The upward propagation of a lean methane-air ame in a vertical tube is investigated. The shape and the velocity of the ame front are extracted from video records, and the velocity of the ow induced by the ame is measured with PIV using alumina particles to seed the gas. The thermophoretic drift of these particles relative to the gas is shown to cause an error in the PIV measurements in the transport region of the ame, where the temperature gradient is large. An iterative correction of the measured velocity is pro- posed that uses simpli ed, quasiunidimensional energy and species conservation equations to compute an approximation to the pro le of gas temperature across the ame in terms of the pro le of gas velocity. The correction is tested using synthetic velocity elds and applied along the axis of the tube. The results of the quasiunidimensional model also show that a conduction-radiation balance is approached in a region of low velocity that develops behind the ame front when the ammability limit is approached. Radiation losses from the combustion products in this region become then important and may cause the extinction of the ame. Keywords: Lean premixed ames; upward propagation; thermophoresis INTRODUCTION In an attempt to normalize the determination of ammability limits, which in practical premixed systems depend on the conditions of the ow in addition to the kinetics of the combustion reactions, Coward and Jones (1952) carried out a comprehensive analysis of previous work and proposed a standard test to characterize reactive mixtures. A vertical tube of 51 mm ID and 1.8 m long, open at the bottom and closed at the top, is lled with the mixture to be tested, which is ignited near the lower end of the tube. The mixture is said to be ammable if a ame is generated and propagates to the upper end of the tube. Levy (1965) noted that the shape and velocity of near limit ames propagating in this standard ammability tube are dominated by buoyancy forces and are similar to the shape and velocity of a bubble rising in the tube, except for extreme values of the Lewis number of the reactant that is depleted by the ame. Extinction at the ammability limit begins at the tip of the ame front. The limit concentration of the de citary reactant changes when the radius of the tube is changed, in a manner that depends on the Lewis number of this reactant and possibly on other factors (Coward and Jones 1952, Babkin et al. 1982). Lewis and von Elbe (1961) and Jarosinski et al. (1982) pointed to ame stretch as a cause of extinction at the ammability limit, and von Lavante and Strehlow (1983) approximately computed the ow of the fresh gas above the ame and found that the stretch is mainly due to the strain rate of the ow, which scales with the inverse of the square root of the tube radius, and that it is maximum at the tip of the ame front and of the right order to extinguish the ame (see Higuera 2011 for a more recent computation that includes the ow on both sides of the ame front). The issue of preferential di usion in stretched ames was taken up by Shoshin et al. (2008). Since a positive stretch strengthens the ame when the Lewis number of the reactant that is depleted is below unity (Williams 1985), the limit concentration of such reactant should Seeding system High speed camera FlameLaser sheet CH4 N2O2Impactor Figure 1: Sketch of the experiment decrease with the radius of the tube. However, from their analysis, Shoshin et al. (2008) concluded that these e ects cannot always explain the observed trends. Later, Shoshin and Jarosinski (2009) and Shoshin et al. (2010) proposed that heat losses from the ame due to the e ect of radiation in a low velocity region that appears below very lean ames may explain the observed extinction. Thermophoretic forces acting on small particles suspended in a gas whose temperature is not uniform cause these particles to drift away from hot regions of the gas with a velocity proportional to the local temperature gradient. When this velocity is not small compared to the velocity of the gas, the e ect is detrimental to the accuracy of LDV and PIV measurements, which rely on the assumption that seeding particles follow the motion of the gas.The importance of thermophoresis in ame-related studies was demonstrated by Rosner et al. (1992) and Gomez and Rosner (1993). Discrepancies between measured and computed velocity pro les in strained ames were investigated, among others, by Chelliah et al. (1990) for di usion ames and Sung et al. (1996) for premixed and di usion ames. In this paper, we report an experimental investigation of very lean methane-air ame fronts propagating upward in a tube. The shape and velocity of the ame front and the velocity of the gas around the ame are measured using photographic and PIV techniques. An approximate correction of the PIV results to account for the thermophoretic drift of the seeding particles is proposed that is based on simpli ed, quasiunidimensional model conservation equations for the energy and the reactants across the ame. The model determines an approximation to the temperature pro le and the thermophoretic drift across the ame in terms of the gas velocity pro le, which can thus be corrected in an iterative manner. The model is also used downstream of the ame, where it shows that radiation losses from the combustion products in a region of low velocity (relative to the ame) that appears behind very lean ames may account for the extinction of the ame at the ammability limit. EXPERIMENTAL SETUP AND TECHNIQUES The experimental apparatus has been described elsewhere (Muntean and Higuera 2015); only a brief summary is given here. A quartz tube of 54 mm internal diameter and 1.9 m long is lled from the upper end with a mixture of 99.5%-pure methane and dried air or synthetic air made from ultrapure oxygen and nitrogen. The mixture, with a selected equivalence ratio, ows through the tube for at least 3 minutes in order to evacuate residual gases from previous runs. Mass ow controllers (EL-Bronkhorst calibrated with Ritter TG05 or TG10 drum-type gas meters) are used to separately measure the ow rate of each gas, achieving an accuracy of 1% in the equivalence ratio. The upper end of the tube is closed after it is lled with the fresh mixture, and the gas in the lower 15 cm of the tube is replaced with a slightly enriched mixture that enters through the lower valve in the sketch of Fig. 1. The purpose of this mixture is to facilitate launching the ame near the ammability limit; see Levy (1965), Shoshin et al. (2008) and Shoshin and Goey (2004). When the process is complete, the mixture is allowed to settle for about two minutes and is ignited with an electrical resistance located near the lower end of the tube, which remains open. The shape and velocity of the ame front that rises in the tube are extracted from the video record of a camera that images a region of the tube about 15 cm long at 100 frames per second with a resolution of 15{20 pixel/mm. The Canny algorithm with a Laplacian of Gaussian lter is used to extract the contour of the ame from the recorded images. A PIV system consisting of a continuous 808 nm laser (Amtron LS453) with a cylin- drical sheet-forming lens and a high speed camera (NAC Memrecam HX-3) with a Nikkor Micro 60 mm objective is used to measure the velocity of the gas around the ame. The laser sheet is shone through the vertical axis of the tube, and the camera images a region of about 68 mm with a resolution of 32 pixel/mm. Alumina particles (Metalograph) are used to seed the gas. These particles are produced in the uidized bed sketched in the upper part of Fig. 1 and enter the tube with the fresh mixture after going through an im- pactor that suppresses the largest particles. Particles with nominal diameters of 0.3 m, 1m and 5m have been tried. The actual size distributions of the particles downstream of the impactor have been measured with DMA (TSI SMPS3938) and laser di raction (Sympatec) size spectrometers, and from SEM micrographs. These measurements reveal fairly wide distributions of diameter, with a single maximum at a diameter somewhat smaller than the nominal value given by the supplier. The distributions extend into the tens of nanometers range, but these small particles are not expected to play a role in our PIV measurements, as they are well into the range of Rayleigh scattering where the in- tensity of scattered radiation decreases as the sixth power of the diameter scaled with the wavelength of the laser light. Globular agglomerates of tens of particles, with diameters of the order of the nominal value, also appear in the SEM images, as well as some larger aggregates which for the most part should be suppressed by the impactor. Particles of 0.3m nominal diameter were discarded as they are poor scatterers of the laser light used. Thermophoretic e ects are important for all the particles used, which therefore fail to follow the gas in the transport region of the ame where the temperature gradient is large. An analysis of this e ect and a possible correction of the PIV measurement are presented in the following section. THERMOPHORESIS Thermophoretic force and drift velocity The thermophoretic force acting on a spherical particle of radius aimmersed in a gas is of the form FTP=afTPrlnT; (1) where==;,Tandare the density, temperature and viscosity of the gas; and fTP(Kn;) is a function of the Knudsen number, Kn = `=a, where`= 2=c with c=p 8RgT= is the mean free path of the gas molecules, and of the gas-to-particle con- ductivity ratio . Simple expressions for this function exist for small and large values of the Knudsen number, in the hydrodynamic and free molecular regimes, whereas correla- tions are available for Kn = O(1). A simple correlation proposed by Talbot et al. (1980) that spans the whole range of Knudsen numbers is fTP=12Cs( +CtKn) (1 + 3CmKn) (1 + 2 + 2 CtKn); (2) whereCs= 1:17,Ct= 2:18 andCm= 1:14. Alternative, more accurate results have been proposed by Beresnev and Chernyak (1985), but the simple expression (2) will be assumed here. The hydrodynamic force acting on the particle due to its relative motion to the gas is, for small values of the gas-to-particle density ratio and of the Reynolds number of the slip ow, Fv=6a fC(Kn)(vpv) withfC= 1 + Kn [A+Bexp (C=Kn)]; (3) where vpandvare the velocities of the particle and the gas, and fC(Kn), with A= 1:20, B= 0:41 andC= 0:88, is the Cunningham correction factor to Stokes law for non-small values of the Knudsen number. Using (1) and (3), the equations of motion of the particle become mdvp dt=Fv+FTP+mganddxp dt=vpwithm=4 3a3p; (4) where xpis the position of the particle, gis the acceleration of gravity, and pis the density of the solid. When the e ects of the inertia of the particle and the gravity can be neglected, the balance Fv+FTP= 0 gives the thermophoretic drift velocity of the particle relative to the gas as vTPvpv= TPrlnTwith TP=fTPfC 6: (5) The coecient TPtends to 0.55 in the free molecular limit Kn !1 and takes very small values in the hydrodynamic limit Kn !0 when  1. This coecient is shown in Fig. 2(a) as a function of the gas temperature for alumina particles of various radii (density 4000 kg/m3, thermal conductivity 30 W/m K). Here the approximation =u= (T=Tu), with= 0:75,Tu= 300 K and u= 1:78105kg/m s, has been used for the gas viscosity. The Knudsen number of a particle of 1 m radius increases from 6 :28102 to 0.47 when the gas temperature increases from 300 K to 1500 K, and TPincreases accordingly. In the transport region of a ame, rlnT1=, whereis the ame thickness, which satis es =ULin terms of the burning velocity UL(Williams 1985). Thus vTP TP= TPUL; the velocity of the particles signi cantly departs from the velocity of the gas when TP=O(1), which happens for small particles at high gas temperatures. Synthetic velocities To assess the e ect of thermophoresis in the conditions of our experiment, the trajectories of seeding particles have been computed by solving (4) with synthetic velocity and tem- perature elds generated by a model of the ow and the ame. The model assumes that (b) 0 0.1 0.2 0.3 0.4 0.5 1 2 3 4 5 0.8 0 1 0.5 1.2 1 1.4 1.6 -8 -4 0 4(a) 0.8 1.2 1.6 2 -8 -4 0 4(c) v/UL1 2 3 4 0 1 2 3 -8 -4 0 4 x/δLx/δLT/Tu x/δL θv/ULαT P ρv/ρ uUL(d) v/UL θ 05Figure 2: (a) Coecient TPas a function of the gas temperature scaled with Tu= 300 K for a= 1, 2, 5 an 10 microns, increasing from top to bottom. ( b) Pro les of gas velocity (black) and temperature ( , red, right-hand side scale) along the axis of the tube, computed from (6){(11) for = 0:546, and pro les of particle velocity for a= 1m (green) and 5 m (blue), computed from (4). The blue dashed curve is the velocity obtained from the balance Fv+FTP= 0 for 5m particles. Velocities are scaled with the reference burning rate UL= 4:63 cm/s, distance along the axis is scaled with the reference ame thickness L=DTu=UL= 0:43 mm, and = (TTu)=(TeTu), whereTeis the equilibrium temperature of the mixture computed with the Cantera software. ( c) Corrected gas velocity for = 0:546 from successive iterations for a= 1m seeding particles. The red curve, labeled 0, is the uncorrected (particle) velocity and the dashed curve is the true gas velocity. ( d) Measured (blue, solid) and corrected (blue, dashed) velocity pro les for = 0:54. Red and green curves are the temperature and mass ux pro les computed with the quasiunidimensional model for the measured (solid) and corrected (dashed) velocity pro les. combustion occurs through a single irreversible Arrhenius reaction, uses Fick's di usion law for all the species, and accounts for radiation losses from the CO 2and H 2O generated by the ame by assuming that the gas is optically thin and using the statistical narrow band model to evaluate the Planck mean absorption coecient (Fiveland 1984, Sou ani and Taine 1997). Assuming, in addition, that the ow is axisymmetric and stationary in a reference frame rising with the ame front, the governing equations take the form (see Higuera and Muntean 2014 for details) r(v) = 0; vrv=rp+r 0+g; T =uTu (6) cpvrT=r(cpDTrT) +qwL; vrYi=rDT LeirYi siw; (7) fori= CH 4, O 2, CO 2and H 2O, with the boundary conditions x!1 :v=U0i; T=Tu; YCH4= 17:39 + ; YO2=4 17:39 + ; YCO2=YH2O= 0 (8) x!1 :@v @x= 0;@T @x=@Yi @x= 0 (9) r=R:v=U0i; T=Tu;@Yi @r= 0; (10) and w=Bexp (Ta=T)YCH4YO2; u=DT uDTu=T Tu L= 4K T4T4 u withK=R0T KCO2YCO2 WCO2+KH2OYH2O WH2O! :9 >>>= >>>;(11) Here,T,pandvare the density, temperature, pressure and velocity of the gas; 0is the viscous stress tensor given by the Navier-Poisson law; Yi,Wiand Leiare the mass fractions, molecular masses and Lewis numbers of the di erent species; andDTare the gas viscosity and thermal di usivity, which are taken to increase as powers of the temperature, and cpis the speci c heat, which is taken to be constant; xandrare downward distance along the axis of the tube and distance to the axis, respectively, and iis a unit vector pointing downwards; the subscript udenotes conditions of the fresh gas above the ame, and = 4YCH4u=0:23 1YCH4u is the equivalence ratio; U0is the upward velocity of the ame front, which is an eigenvalue of the problem; Ris the radius of the tube; qis the heat released per unit mass of methane burnt, and sCH4= 1,sO2= 4, sCO2=11=4 andsH2==9=4; the frequency factor Band the activation temperature Ta are chosen for the single reaction model to reproduce the planar ame velocity computed with the GRI-Mech 3.1 scheme in the range of equivalence ratios of interest (Higuera and Muntean 2014); KCO2andKH2Oare functions of the temperature computed and tabulated by Ju et al. (1999); and andR0are the Stefan-Boltzmann and universal gas constants. Values of the parameters involved are DTu= 2105m2/s,cp= 103J/kg K, LeCH4= 1, LeO2= 1:05, LeCO2= 1:39, LeH2O= 0:83,q= 5:01104kJ/kg,Ta= 18750 K, Bsuch thatUL= 4:63 cm/s for YCH4u= 0:03 (Davis et al. 2002), and R= 27 mm. Figure 2(b) shows pro les of gas velocity (black) and temperature (red) along the axis of the tube in the region of the ame for a value of the equivalence ratio = 0:546. The velocity of the gas rst decreases on approaching the ame front from above, then increases on crossing the ame front due to the gas thermal expansion, decreases again downstream of the ame due to the buoyancy force acting on the burnt gas, and nally increases at some distance downstream of the ame front (not shown in the gure) due to the convergence of the burnt gas that crosses the ame at di erent points of the front. The velocity of a seeding particle along the axis is computed from (4) assuming that the initial particle velocity (relative to the ame) is equal to velocity U0of the gas. The result is also shown in Fig. 2(b) for particles of radii a= 1m (green curve) and 5 m (blue curve). As can be seen, the smaller particles are signi cantly delayed by thermophoresis inside the ame, where the temperature rapidly increases in the streamwise direction and FTP<0, and they move slightly faster than the gas downstream of the ame, where the temperature slowly decreases in the streamwise direction due to radiation losses. The e ects of the inertia and the gravity are negligible for these particles; the green curve is indistinguishable from the velocity obtained from the balance Fv+FTP= 0. The larger particles are less a ected by thermophoresis, but the e ects of the inertia and the gravity are important for them. These particles always fall faster than the gas, and the velocity di erence is largest upstream of the ame front, where the gas undergoes a continuous deceleration that the large particles fail to follow. The blue dashed curve shows the velocity obtained from the balance Fv+FTP= 0. Correction of thermophoretic drift PIV results obtained with small seeding particles for which the e ects of the inertia and the gravity are negligible could be corrected for thermophoretic drift if the temperature distribution of the gas were known. In the absence of this information, a simple procedure is proposed here that allows us to compute an approximation to the gas temperature using simpli ed forms of the energy and species conservation equations. These equations require the velocity of the gas as an input which should come from PIV measurements. Since the correct gas velocity is not known beforehand, an iteration is used to simultaneously re ne the correction of the measured velocity and the approximate temperature distribution on which this correction is based. The simpli ed conservation equations can in principle be applied to any section of the ame front, using information on the local normal. However, they are used here only at the uppermost point of the front, where the normal and the velocity of the gas, v(x) in what follows, are vertical. The simpli ed equations are derived from equations (7) taking advantage of the small thickness of the ame compared to the curvature radius of the ame front (of the order of the radius of the tube) to neglect radial di usion of species everywhere in the ow eld and radial conduction of heat in the transport region of the ame, while the radial conduction term of the energy equation is approximated downstream of the ame by cpDT (TbT)=R2, whereTbis the nal combustion temperature, Ris the radius of the tube, and = 15 is a model parameter whose value is chosen for the temperature obtained from the simpli ed equations to match the solution of (6){(11) in a number of cases (see Muntean and Higuera 2015). The simpli ed equations are cpvdT dx=d dx cpDTdT dx +cpDT H(T)TbT R2+qwL; vdYi dx=d dxDT LeidYi dx siw;(12) whereH(T) is an indicator function equal to zero upstream of the temperature maximum an equal to one downstream. These equations are to be solved with the boundary condi- tions (8) and (9) for the temperature and the species mass fractions, and with the reaction and radiation terms given by (11). For a given v(x), the equivalence ratio is iteratively computed for the mass ux vto vary smoothly across the ame. The di erence between the value of the equivalence ratio thus determined and its actual experimental value is a measure of the accuracy of the model. The iterative correction of the PIV measured velocity proceeds as follows. The mea- sured velocity, vm(x) say, is rst used for v(x) in (12). The temperature pro le computed from these equations is used to evaluate the thermophoretic drift velocity vTP(x) from (5) and compute a corrected gas velocity v(x) =vm(x)vTP(x). This is carried again to the simpli ed equations (12) and the computation is iterated until the corrected gas velocity converges. As a test, the iterative procedure has been applied to the synthetic particle velocity determined from (4) with the gas velocity and temperature elds obtained from (6){(11). Assuminga= 1m, the velocity of the particles is given by the green curve of Fig. 2(b), which plays the role of vm(x), while the true gas velocity is given by the black curve. Figure 2(c) shows the corrected gas velocities obtained in successive iterations. As can be seen, the corrected gas velocity is indistinguishable from the true velocity (dashed curve in Fig. 2(c)) after six iterations. RESULTS AND DISCUSSION Steadily rising ame fronts consisting of a rounded cap and a nearly cylindrical skirt are obtained after an initial transient when wall temperature perturbations and residual motions of the gas are minimized. As it has been often noted, the shape of the upper cap and the velocity of very lean ame fronts are similar to those of a bubble rising in the tube. The measured ame front velocity is U0= 24 cm/s and depends little on the equivalence ratio in a range of values of this variable around the ammability limit, which is about = 0:525 in our experiment. For comparison, the velocity of a bubble rising in this tube is 22 cm/s according to Davies and Taylor (1950) formula for a heavy-to-light uid density of 5, while the burning velocity ULof a planar ame with the same equivalence ratio is much smaller, of only a few centimeters per second. This disparity of velocity scales allows to explain the observed shape and velocity of the ame front. Buoyancy force pushes the burnt gas upward so that, in a reference frame moving with the ame front, the incoming fresh gas is decelerated and moves radially outward on approaching the front, much as the ow above a bubble. Also as for a bubble, pressure variations in the burnt gas behind the ame are small compared to the dynamic pressure of the fresh gas, due to the small density of the burnt gas. Finally, due to the small burning velocity of the ame, its tip lies close to the stagnation point of the equivalent bubble ow, and the rest of the ame is close to a streamline through this point. The length of the skirt increases when the equivalence ratio, and thus the burning velocity of the ame, is decreased. Figure 3 shows the velocity of the gas relative to the ame for = 0:53 (near the ammability limit) and 0.54. These velocity elds are obtained by substracting the upward velocity of the ame front from the PIV measured velocity. The results are not corrected for thermophoresis, whose e ect, however, is small outside of the ame. Comparison of the two elds shows that the ow above the ame front displays the features discussed above and is nearly independent of the equivalence ratio, while the ow below the ame front develops a region of low velocity (relative to the ame) as the ammability limit is approached. The size of this region scales with the radius of the tube, and in some cases, like the one shown in Fig. 3(a), includes recirculation of the burnt gas. The residence time of the gas in the region of low velocity can be very large. As will be seen below, this increases the importance of radiation losses (mainly from CO 2and H 2O) to such an extent as to a ect the ame and eventually cause its extinction, despite the small concentrations of the radiating species. The measured velocity of the gas along the axis of the tube has been corrected for r (cm) x (cm) −2 −1 0 1 243210−1 −2 −1 0 1 243210−1 r (cm) x (cm)(a) (b)Figure 3: PIV velocity eld of the gas relative to the ame for = 0:53 (a) and 0.54 ( b). The black curve is the contour of the luminous region of the ame extracted from the video records. The red contour encloses the region of low velocity where jvj<6 cm/s. The velocity is negative (upward) in the hatched region of gure ( a). -0.004-0.002 0 0.002 0.004 0 10 20 30 40 x/δLConductionRadiation Convection Figure 4: Convection (red), conduction (green) and radiation (blue) terms of the energy equation (12), all three scaled with uULcp(TeTu)=L, as functions of dimensionless distance from the point of maximum temperature, for = 0:53 (solid curves) and 0.55 (dashed curves). thermophoretic drift using the procedure described in the previous section. For this correction, the radius of the particles is taken to be a= 1m, as particles of this size are both abundant and more ecient laser light scatterers than smaller particles also present in the seeding. Figure 2(d) shows the e ect of the correction for = 0:54. The solid blue curve is the measured (uncorrected) velocity, and the dashed blue curve is the corrected velocity. The iteration scheme converged after four iterations. The solid and dashed red curves are the temperature pro les given by the quasiunidimensional model for the uncorrected and corrected gas velocity, respectively. The solid and dashed green curves are the corresponding mass uxes, v. As can be seen, the magnitude of the velocity correction is signi cant within the ame, and the corrected mass ux is smoother than the uncorrected one. The values of the equivalence ratio yielded by the quasiunidimensional model are model = 0:521 and 0.539 for the uncorrected and corrected velocity pro les. The second one is nearer the real experimental value, suggesting that the correction has improved the velocity pro le. Computations with the quasiunidimensional model can be extended downstream of 0 0 0.001 1 0.002 2 0.003 3 0.004 4 0.005 5 -6 -4 -2 0 2 4 6 T/T u v/UL x/δLωiMCH4δL//parenleftBig ρuULYCH4u/parenrightBigFigure 5: Pro les of temperature (red, right-hand scale) and rates of the radical production (green) and the radicals recombination (magenta) steps of the four step reduced scheme (left-hand scale) computed as functions of downward distance using the corrected velocity pro les (dashed blue, right-hand scale) for = 0:53 and 0.55, increasing from bottom to top or as indicated by the arrow. The values of the equivalence ratio determined from the quasiunidirectional model are model = 0:52 and 0.53, respectively. Temperatures are scaled with Tu= 300 K, velocities with the reference burning rate UL= 4:63 cm/s, distances with the reference ame thickness L=DTu=UL= 0:43 mm, with DTu= 0:2 cm2/s, and reaction rates with u YCH4u=WCH4 UL=L. the ame. Figure 4 shows the convection, conduction and radiation terms of the en- ergy equation (12) evaluated in this region for = 0:53 (solid curves) and 0.55 (dashed curves). At the lower value of the equivalence ratio, the small velocity of the gas down- stream of the ame front renders convection negligible and reduces the energy equation to the conduction-radiation balance rq+L0, with q=cpDTrT, in a large region of the ow eld. This balance is quite di erent from the convection-radiation bal- ance characteristic of a planar ame with heat losses (Williams 1985), a balance that is already approached by the results for the higher equivalence ratio in Fig. 4. The order of the heat ux in the conduction-radiation region is qbLbR, whereLbis the radia- tion loss term in (11) evaluated in the burnt gas. The decrease of the nal combustion temperature obtained by using this estimation in an enthalpy balance across the ame is T=O(LbR=ucpUL), whereULis the ame burning velocity (see Muntean and Higuera 2015). This temperature decrease may easily be of the order of the Frank-Kamenetskii temperature T2 b=Taand may account for the extinction of the ame at the ammability limit. The quasiunidimensional model of the previous section has been extended to replace the single step Arrhenius reaction by a more realistic four step reduced kinetic mecha- nism; see S anchez et al. (2000) and