Dianailys Nuñez-Reyes

Kevin M Hickson

A Low Temperature Investigation of the Gas-Phase N( 2 D) + NO Reaction. Towards a Viable Source of N( 2 D) Atoms for Kinetic Studies in Astrochemistry

The gas-phase reaction of metastable atomic nitrogen N( 2 D) with nitric oxide has been investigated over the 296 -50 K temperature range using a supersonic flow reactor. As N( 2 D) could not be produced photolytically in the present work, these excited state atoms were generated instead through the C( 3 P) + NO ® N( 2 D) + CO reaction while C( 3 P) atoms were created in-situ by the 266 nm pulsed laser photolysis of CBr4 precursor molecules. The kinetics of N( 2 D) atoms were followed on-resonance by vacuum ultraviolet laser induced fluorescence at 116.7 nm. The measured rate constants for the N( 2 D) + NO reaction are in excellent agreement with most of the earlier work at room temperature and represent the only available kinetic data for this process below 296 K. The rate constants are seen to increase slightly as the temperature falls to 100 K with a more substantial increase at even lower temperature; a finding which is not reproduced by theoretical work. The prospects for using this chemical source of N( 2 D) atoms in future studies of a wide range of N( 2 D) atom reactions are discussed.

Introduction

Excited state metastable nitrogen atoms N( 2 D) are thought to play important roles in the chemistry of a range of planetary atmospheres. In the terrestrial upper atmosphere these atoms are produced by various processes including the dissociative excitation of N2 by energetic electrons and by short wavelength N2 photolysis 1 . The reaction between N( 2 D) with atmospheric O2 is considered to be a major source of both NO and metastable O( 1 D) atoms 2 in the Earth's lower thermosphere around 110 km, while the reaction between N( 2 D) and NO limits the total concentration of N( 2 D) at these altitudes. As N2 is an important component of the atmospheres of Mars and Titan, the reactions of N( 2 D) are also important processes in the photochemistry of both these planets. Indeed, N( 2 D) reacts rapidly with CO2 in the Martian atmosphere to form NO and CO 3 whereas the reactions between N( 2 D) and the various saturated and unsaturated hydrocarbons present in the atmosphere of Titan are considered to be major sources of N-bearing organic molecules [4][5][6] . From an experimental point of view, there are several major difficulties associated with the investigation of the kinetics of N( 2 D) reactions.

Firstly, convenient photolytic sources of N( 2 D) atoms are very rare, usually requiring radiation in the vacuum ultraviolet (VUV) wavelength range to dissociate nitrogen containing precursor molecules. Husain et al. 7,8 employed the flash photolysis of N2O at wavelengths greater than 105 nm in their experiments whereas Umemoto et al. 9 used the multiphoton dissociation of NO molecules with photons around 275.3 nm as the source of N( 2 D) in their experiments. On the detection side, the radiative transition from the 2 D state to the ground 4 S state is too weak to produce measureable emission in the visible region [START_REF] Ralchenko | NIST Atomic Spectra Database (ver. 4.0.1)[END_REF] so that other methods are generally employed to follow the progress of N( 2 D) reactions such as electron spin resonance [START_REF] Evenson | [END_REF]12 , mass spectroscopy 13 and resonance absorption 7,14 and emission spectroscopy 15,16 . More recently, Umemoto et al. 17 detected N( 2 D) atoms by the 269 nm two-photon excitation of N atoms through the 2p 3 2 D5/2 -2p 2 3p 2 S transition. Rapid relaxation to the 2p 2 3s 2 P state followed by VUV fluorescence at 149 nm allowed these authors to determine rate constants for various N( 2 D) reactions at room temperature. One system that has received considerable experimental and theoretical attention in the past is the reaction between N( 2 D) and NO

N( 2 D) + NO ® N2 + O( 3 P, 1 D, 1 S) (1)
On the experimental side, earlier measurements at room temperature 7,8,15,[17][18][19] have determined rate constants for this process in the range (0.35 -1.8) × 10 -10 cm 3 s -1 with a recommended value of 𝑘 #( % &)(#) (298 K) = 6 × 10 -11 cm 3 s -1 20 . The equivalent reaction of ground state atomic reaction N( 4 S) is thought to be slower at room temperature and below 21 . To the best of our knowledge, no previous experimental studies of reaction (1) have been performed at other temperatures. On the theoretical side, Gonzalez et al. 22 surface. They obtained a rate constant value for this channel of (2.63 ± 0.04) × 10 -12 cm 3 s -1 , more than an order of magnitude lower than the lowest experimental value for the total rate constant at room temperature. More recently, Li et al. 23 performed QCT calculations over the 5 -3000 K range using a double many-body expansion potential energy surface for ground state N2O( 1 A¢) with an accurate description of the long-range interactions 24 . The derived rate constants were seen to be larger than those obtained by Gonzalez et al., with

𝑘 #( % &)(#)® )( / &)(# % (300 K) = (2.03 ± 0.01) × 10 -11 cm 3 s -1 .
The rate constant was found to decrease with temperature, in contrast with the work by Gonzalez et al. 22 who determined the opposite temperature dependence.

Here we report the first measurements of the total rate constant for the N( 2 D) + NO reaction over the 50 -296 K range using a continuous supersonic flow reactor. In the absence of efficient photolytic schemes for N( 2 D) production in accessible wavelength ranges, these atoms were generated instead by chemical reaction. N( 2 D) atoms were detected directly by pulsed laser induced fluorescence in the vacuum ultraviolet wavelength range. The methods used in this investigation are described in Section 2. Section 3 presents the data analysis methodology and the results obtained. The conclusions are outlined in Section 4.

Experimental

A continuous supersonic flow (Laval nozzle) reactor 25 was used in the present investigation.

Detailed descriptions of the main features of the apparatus have been provided in earlier work 26,27 with other modifications specifically related to the efficient generation of short wavelength tunable VUV radiation for the detection of atomic species reported in later studies [28][29][30][31] a The errors on the Mach number, density, temperature and mean flow velocity (1σ) are calculated from separate measurements of the impact pressure using a Pitot tube as a function of distance from the Laval nozzle and the stagnation pressure within the reservoir.

In addition, measurements were also performed at room temperature (without a Laval nozzle) through a procedure which used lower flow velocities. In these types of experiments, the minor reagent radical species are normally produced in-situ within the supersonic flow through the photolysis of suitable precursor molecules. The photolysis laser is aligned along the axis of the flow to create a column of radicals with uniform density along the entire length of the flow (assuming only weak absorption along the length of the flow). In the case of N( 2 D) atoms however, there are no known precursor molecules in the UV wavelength range, except through multiphoton processes. In this way, Umemoto et al. 17 increasing to slightly to 2.0 × 10 -10 cm 3 s -1 at 50 K. 34 According to previous theoretical work, channel (2a) is expected to decrease as the temperature falls. 35 C( 3 P) atoms were created by the 10 Hz pulsed sequential multiphoton dissociation of carbon tetrabromide, 36 CBr4, at 266 nm with pulse energies around 24 mJ. In this case, the photolysis laser was unfocused. CBr4 molecules were introduced into the reactor by flowing a small fraction of the carrier gas (Ar or N2) over solid CBr4 held in a separate flask at a fixed temperature of 298 K. Gas-phase concentrations of CBr4 were lower than 3.2 × 10 13 cm -3 in the supersonic flow, as estimated using the saturated CBr4 vapour pressure at room temperature. Previous work performed under similar conditions 28 showed that the major photolysis product is C( 3 P), with a small non- The PMT output was fed into a fast preamplifier which was connected to a boxcar module for signal processing and integration. The delay time between photolysis and probe lasers was scanned by a digital delay generator. 30 laser shots were recorded at each delay time with at least 90 time points for each decay profile. Laser shots were also recorded at negative time delays to set the fluorescence signal baseline level. All gases used in the experiments were flowed from cylinders with no further purification prior to use, while the flows of the carrier gases and reactive gases were regulated by calibrated mass-flow controllers. The suppliers and purities of the gases used are listed below (Linde Ar 99.999%, Xe 99.999%, Air Liquide N2 99.999%, NO 99.9%).

Results and Discussion

In order to apply the pseudo-first-order approximation, the concentration of NO was maintained

I #( % &) = A(exp(-𝑘 < = 𝑡) -exp(-𝑘 ? = 𝑡)) (3)
where 𝐼 #( % &) is the time dependent N( 2 D) VUV LIF signal, A is a constant, 𝑘 < = is the pseudofirst-order rate constant for N( 2 D) loss, 𝑘 ? = is the pseudo-first-order rate constant for N( 2 D) formation and t is the time. As the N( 2 D) VUV LIF signal levels were small compared with those recorded during our earlier studies of C( 3 P), O( 1 D) and H( 2 S) reaction kinetics, it was necessary to amplify the PMT output signal using a fast preamplifier. Although the PMT was supposed to be insensitive to radiation in the UV range, with an expected cutoff around 200 nm, scattered light (and/or window fluorescence) generated by the photolysis laser lead to saturation of the PMT for the first 15 microseconds after the photolysis laser pulse.

Consequently, the expected rise of the N( 2 D) VUV LIF signal could not be recorded close to zero. Instead, a single exponential decay function of the form

𝐼 #( % &) = Aexp(-𝑘 < = 𝑡) (4)
was used to fit the N( 2 D) VUV LIF signals as shown in Figure 1. Here,

𝑘 < = = 𝑘 #( % &)(#) [NO] + 𝑘 #( % &)(F [X] + 𝑘 HIJJ,#( % &) where 𝑘 #( % &)(#)
is the second-order rate constant for reaction (1), 2) is thought to be somewhat larger than the rate constant for reaction (1) (at least at room temperature 20,34 ), simulations indicated that N( 2 D) formation was complete (assuming a delay corresponding to four half-lives had elapsed) for all but a few kinetic profiles within the first 100 microseconds after C( 3 P) formation.

𝑘 #( % &)(F [X]
Decays similar to those shown in Figure 1 were recorded for a range of [NO], allowing a wide range of 𝑘 L = values to be derived. Plots of 𝑘 L = as a function of the corresponding [NO] at any given temperature allowed the second-order rate constant to be determined through weighted linear least-squares fits to the datapoints. Figure 2 shows examples of these second-order plots obtained at 296 K, 75 K and 50 K. The solid red, green and blue lines represent linear least-squares fits to the datasets, weighted by the statistical uncertainties of the individual datapoints (obtained through single exponential fits of expression (4) to curves similar to those shown in Figure 1).

The second-order rate constants are displayed as a function of temperature in Figure 3 alongside earlier experimental and theoretical work. c Uncertainties on the measured rate constants represent the combined statistical (1s) and estimated systematic (10%) errors.

Both statistical and systematic uncertainties contribute to the error limits on these values.

Statistical errors were obtained through the weighted linear-least squares fitting procedure used to derive second-order rate constants (see Figure 2). An estimated systematic error of 10% was also included, a value which was chosen to account for the possible contributions of errors in the calibration of the mass-flow controllers and the subsequent determinations of the concentrations of the reactants and the supersonic flow density as well as possible errors in pressure gauges and other instruments.

On the experimental side, the measured rate constant at 296 K, 𝑘 #( % &)(#) (296 K) = (7.3 ± 0.7)

× 10 -11 cm 3 s -1 is in good agreement with all but two of the previous determinations of the rate constant at room temperature, 7,8,17,19 with values in the range (5.9-8.3) × 10 -11 cm 3 s -1 using a variety of methods to follow the progress of the reaction. Indeed, this good agreement with the majority of earlier studies goes a long way to validating the novel strategy for N( 2 D) production developed in the present investigation. The measurements by Black et al. 18 and Sugawara et al. 15 of 𝑘 #( % &)(#) (298 K) = (1.8 ± 0.5) × 10 -10 cm 3 s -1 and (3.5 ± 0.3) × 10 -11 cm 3 s -1 derived rate constants that are considerably higher and lower respectively than all other values and the recommended room temperature value of 6 × 10 -11 cm 3 s - 1 20 . The low rate constant value reported by Sugawara et al. 15 has been attributed to possible errors in the determination of the NO concentration.

Below room temperature, the present study indicates that the rate constant displays a slight negative temperature dependence over the 296 K -75 K range followed by a more marked increase at even lower temperature, reaching a value of (1.3 ± 0.1) × 10 -10 cm 3 s -1 at 50 K.

On the theoretical side, it can be seen from Figure 3 that the QCT study of the N( 2 D) + NO reaction by Gonzalez et al. 22 predicts a rate constant that is more than an order of magnitude smaller than all the experimental values at room temperature. This work was based on the ground singlet 1 1 A¢ potential energy surface of N2O alone, whereas there are in fact five potential energy surfaces which adiabatically correlate with O( 1 D) + N2 products. In order to rationalize their results, Gonzalez et al. 22 attempted to consider the possible involvement of the other four excited potential energy surfaces. They concluded that the 2 1 A¢ and 1 1 A¢¢ surfaces might also contribute to the reactivity which could increase the rate constant to 7.9 × 10 -12 cm 3 s -1 assuming similar values of the rate constant for these two surfaces. As further contributions to the overall rate constant are also expected from the channels leading to O( 3 P) + N2 and O( 1 S) + N2 as products, Gonzalez et al. 22 stated that inclusion of these other channels should lead to even better agreement between the experimental and theoretical rate constants. The more recent theoretical investigation by Li et al. 23 quadrupole and the NO quadrupole, Li et al. 23 attributed the improvement over the calculations of Gonzalez et al. 22 as being due to an improved description of the long-range forces.

Nevertheless, the calculated rate constants of Li et al. 23 allow the second-order rate constants to be determined from the slope in an identical manner to the usual analysis (see Figure 2). In this instance, the y-axis intercept value will correspond to the term 𝑘 #( % &)(#) [NO] (assuming that diffusional losses are small in comparison).

Preliminary measurements in our laboratory indicate that should be possible to use this procedure to study the kinetics of a wide variety of N( 2 D) atom reactions down to low temperature if the excess reagent species does not react too rapidly with ground state atomic carbon.

Conclusions

Here we present an experimental investigation of the temperature dependent kinetics of the reactions of metastable nitrogen atoms in the 2 D state with nitric oxide. A Laval nozzle reactor was used to study this gas-phase reaction at temperatures as low as 50 K. As N(

  investigated the channel leading to N2 and O( 1 D) as products, through ab initio CASSCF/CASPT2 electronic structure calculations of the ground 1 A¢ potential energy surface coupled with quasiclassical trajectory (QCT) calculations to obtain the rate constant over the temperature range 300 -1500 K. They identified several different reactive mechanisms including attack by the N( 2 D) atom on both the O-atom and N-atom ends of the NO molecule to form products N2 and O( 1 D) over the barrierless 1 A¢

  produced N( 2 D) atoms in their kinetic experiments through the multiphoton dissociation of NO, by focusing a tuneable laser at 275.3 nm into their reactor. Unfortunately, this technique cannot be applied in the present experiments given the requirement to produce a uniform density of radicals along the length of the flow. Instead, N( 2 D) atoms were generated here through the reaction of ground state carbon atoms C( 3 P) with NO molecules C( 3 P) + NO ® N( 2 D, 4 S) + CO (2a) ® O( 3 P) + CN (2b) where the branching ratio [N( 2 D) + N( 4 S)]/[O( 3 P)] has previously been estimated to be 1.5 ± 0.3 at room temperature 33 with a total rate constant 𝑘 2( 3 4)(#) (298 K) = 1.5 × 10 -10 cm 3 s -1

  negligible production of C( 1 D) atoms. Although C( 1 D) atoms also react with NO molecules with rate constants which are smaller than those of the equivalent ground state carbon reaction,37 the products are likely to be similar to those of the C( 3 P) + NO reaction, with the added possibility of O( 1 D) formation. Among the potential reaction products, only secondary reactions of CN radicals and O( 1 D) atoms could potentially interfere with our kinetic study of the N( 2 D) + NO reaction (the O( 3 P) + NO association reaction is slow38 ). The CN + NO reaction is thought to proceed mostly by association at low temperature,39 so this secondary process is also unlikely to interfere with the present measurements. Moreover, the only exothermically available products of the O( 1 D) + NO reaction are N( 4 S) + O2, which are only accessible through nonadiabatic transitions.40 Neither of these products should influence the N( 2 D) kinetics measurements.N( 2 D) atoms were detected directly in the present experiments by on-resonance pulsed laser induced fluorescence in the vacuum ultraviolet wavelength range (VUV LIF) through the 2s 2 2p 3 2 D° -2s 2 2p 2 ( 3 P)3d 2 F transition at 116.745 nm. To obtain tunable light at this wavelength, the 700.5 nm output of a pulsed tunable dye laser was first injected into a BBO crystal, allowing radiation around 350.25 nm to be generated by second harmonic generation with a pulse energy of approximately 8 mJ. After separation of the fundamental radiation, the UV beam was steered and focused into a cell containing 40 Torr of Xe and 560 Torr of Ar for phase matching purposes, generating tunable VUV radiation by frequency tripling. The cell itself was attached to a 75 cm long side arm positioned at right angles to the reactor at the level of the detection region. The side arm was constantly flushed by either Ar or N2 to avoid attenuation of the VUV beam by reactive gases in the chamber. A MgF2 lens at the exit of the cell allowed the VUV beam to be collimated while most of the divergent residual UV radiation was trapped by circular baffles placed along the length of the sidearm. VUV emission from excited N( 2 D) atoms in the cold supersonic flow was collected by a LiF lens and focused onto the photocathode of a solar blind photomultiplier tube (PMT). The PMT was protected from reactive gases in the chamber by a LiF window while the region between the PMT and the LiF window was evacuated to prevent supplementary VUV absorption losses by atmospheric O2.

  in excess with respect to both C( 3 P) and N( 2 D) atoms in these experiments. The variation of the N( 2 D) VUV LIF signal as a function of time is displayed in Figure 1 for two different values of [NO] recorded at 50 K.

Figure 1 N

 1 Figure 1 N( 2 D) VUV LIF signal as a function of time recorded at 50 K. (Red solid circles)

  is a term corresponding to other first-order reactive losses of N( 2 D) atoms (such as through the N( 2 D) + CBr4 reaction for example) and 𝑘 HIJJ,#( % &) represents the first-order diffusional loss of N( 2 D) atoms. In this instance, the starting point for the fits had to be chosen carefully to avoid fitting to the rising part of the temporal profiles. In practice, this was achieved relatively easily by displaying the N( 2 D) VUV LIF signals on a logarithmic scale and by carefully selecting the 'linear' part of the decay curves. Fortunately, as the rate constant for reaction (

Figure 2

 2 Figure 2 Second-order plots for the N( 2 D) + NO reaction. (Red solid circles) data recorded at

Figure 3

 3 Figure 3 Rate constants for the N( 2 D) + NO reaction as a function of temperature. Experimental

  based on a new ground state 1 A¢ potential energy surface is seen to be in better agreement with the measured values. As N2O intermediate formation follows a capture-type mechanism, with the dominant contribution being the electrostatic interactions between the N( 2 D) quadrupole and the NO dipole in addition to the N( 2 D)

  are still four times smaller than those obtained during the present investigation, a discrepancy that almost certainly indicates the role of other channels leading to O( 1 D) + N2, O( 3 P) + N2 and O( 1 S) + N2 in the overall reactivity. Interestingly, the larger measured increase in reactivity below 75 K determined in the present work could be indicative of the increased contribution of another mechanism at low temperature. Indeed, Gonzalez et al. 22 hint at the possible existence of crossing seams between the lowest singlet and triplet adiabatic potential energy surfaces which might become more important as the intermediate N2O lifetime increases. Although the N( 2 D) + NO reaction itself has little or no known application in planetary atmospheres at the present time (other than the Earth's atmosphere), this study could represent an important breakthrough in the investigation of other more important reactions for planetary atmospheres due to its potential as a source of N( 2 D) atoms. Although the production of N( 2 D) atoms through reaction (2a) and their subsequent removal through reaction (1) are coupled processes, we can nonetheless envisage using these reactions to investigate other N( 2 D) reactions such as N( 2 D) + hydrocarbons that are considered to be major sources of N-bearing organic molecules in Titan's atmosphere. Crucially, none of the most abundant saturated hydrocarbon molecules in Titan's atmosphere (CH4 and C2H6) react with ground state atomic carbon, our source for N( 2 D) atoms, so that the yield of N( 2 D) atoms will remain relatively unaffected by the addition of these molecules to the supersonic flow. If the NO concentration is maintained at a constant (small) value it should be possible to extract pseudo-first-order rate constants for the specific hydrocarbon reaction by varying its excess concentration. Here, 𝑘 L = = 𝑘 #( % &)(#) [NO] + 𝑘 #( % &)(N2 [HC] + 𝑘 HIJJ,#( % &) where the term 𝑘 #( % &)(N2 [HC] is now the pseudo-first-order loss rate of N( 2 D) through reaction with the hydrocarbon. Consequently, a plot of the pseudo-first-order loss rate of N( 2 D) as a function of the corresponding [HC] should

  2 D) is very difficult to produce photolytically, a chemical source of these atoms was used in the present work, involving the generation of C( 3 P) atoms through the multiphoton dissociation of CBr4 molecules at 266 nm followed by their reaction with NO. N( 2 D) atoms were detected directly through vacuum ultraviolet laser induced fluorescence. The derived rate constants are seen to be in excellent agreement with most of the previous studies at room temperature. Below room temperature the rate constant displays a slight negative temperature dependence down to 75 K, with a more pronounced increase at lower temperature. This work demonstrates the potential importance of this chemical source of N( 2 D) atoms for future investigations of astrochemically relevant N( 2 D) reactions such as those with saturated hydrocarbons.

Table 1

 1 Characteristics of the supersonic flows

	. In

Table 2

 2 Temperature dependent rate constants for the N( 2 D) + NO reaction Uncertainties on the calculated temperatures represent the statistical (1s) errors obtained from Pitot tube measurements of the impact pressure. b Number of individual measurements.

	T / K	N b	[NO]/ 10 14 cm -3	𝑘 #( % &)(#) / cm 3 s -1

Acknowledgements K.M.H. and D.N.R acknowledge support from the French program ''Physique et Chimie du Milieu Interstellaire" (PCMI) of the CNRS/INSU with the INC/INP co-funded by the CEA and CNES as well as funding from the "Program National de Planétologie" (PNP) of the CNRS/INSU.


Dorian Pena -Éric

Laurendeau -Yannick Hoarau

Development of an icing simulation code in the NSMB flow solver

Keywords:

An icing model is developed in the NSMB (Navier-Stokes Multi-Block) compressible solver.The code implemented in this study is based on an Eulerian formulation for droplets tracking solved implicitly by means of CGSTAB or SIP methods, a modified iterative Messinger model using an improved water runback scheme for ice thickness calculation and three-dimensional mesh deformation to track the ice/air interface through time. The whole process is parallelized with MPI for complex multi-block configurations. The validated on several test cases including NACA23012 and NACA00012 airfoils and ONERA-M6 swept wing. Pressure coefficient around iced 23012 airfoil is compared with experimental data. Rime icing is computed on a DLR-F6 wing/body configuration. In all studied cases, the results are in good agreement with the literature.

I Introduction

Statistics establish that icing and its effects remain one of the main cause of flight accidents. The modification of geometry induced by accreted ice results in large aerodynamic performance degradation. Numerous wind tunnel experiments on iced wings establish a loss of lift, a decrease of stall angle and an increase of drag. The accreted ice eventually breaks and impacts components such engines or fuselage. Anti-icing equipment as well as antennas and Pitot tubes must be particularly well designed to avoid any icing effect.

The combination of time and cost constraints with a more and more demanding certification process has conducted the development of icing simulation tools.

Icing codes are commonly based on inviscid panel or Euler methods to simulate the fluid, a Lagrangian formulation to predict the trajectories of water droplets and impingement areas, and a classical Messinger Model to compute the ice thickness. It is evident that a full fluid solution provided by a Navier-Stokes solver in combination with turbulence model would results in more accurate icing simulations. Icing codes such LEWICE3D (NASA) [START_REF] Bidwell | Users manual for the nasa lewis three-dimensional ice accretion code (lewice 3d)[END_REF],CANICE2D-NS [START_REF] Hasanzadeh | Wing aerodynamic performance analysis and stall prediction using canice2d-ns icing code[END_REF] ( École Polytechnique Montréal),ONERA3D (ONERA) [START_REF] Montreuil | ECLIPPS: 1. three-dimensional CFD prediction of the ice accretion[END_REF] and FENSAP-ICE (McGill University) [START_REF] Beaugendre | Development of a second generation in-flight icing simulation code[END_REF] were developed for that matter.

Lagrangian approach for droplets tracking solves Newton equation of motion for each individual droplet by means of Runge-Kutta integration. This approach requires non-automatic and well defined initial positions and amount of incoming droplets to obtain accurate impingement areas.

Bourgault et al. [START_REF] Bourgault | A finite element method study of eulerian droplets impingement models[END_REF] developed in FENSAP-ICE an Eulerian model for droplets tracking . In this PDE-based formulation, droplet volume fraction and velocity field are solved in the whole computational domain. Solving droplets equation on the same computational grid permits straightforward integration with flow solver.

The key element of ice accretion is the determination of impingement collection efficiency. While a Lagrangian formulation requires a statistical averaging process to evaluate it, the determination of this quantity in Eulerian formulation is direct. In addition, the resolution of droplet equations in the Eulerian frame open up access to the large panel of computational fluid dynamics accelerating techniques. For reasons stated above, the advantages of choosing Eulerian framework are obvious.

In FENSAP-ICE, droplet equations are solved on computational grid by means of weak Galerkin finite element method and implicit time marching using GMRES [START_REF] Bourgault | A finite element method study of eulerian droplets impingement models[END_REF]. Eulerian model for droplets tracking in finite volume formulation have been investigated in many studies [START_REF] Kim | Ice accretion modeling using an eulerian approach for droplet impingement[END_REF][START_REF] Jung | Numerical modeling for eulerian droplet impingement in supercooled large droplet conditions[END_REF][START_REF] Cao | Numerical simulation of ice accretions on an aircraft wing[END_REF][START_REF] Jung | Development of eulerian droplets impingement model using hllc riemann solver and pod-based reduced order method[END_REF], most of the approaches presented in these papers solve equations explicitly by means of Multi-Step Runge-Kutta integration. Depending on the airflow temperature and velocity, impinging water droplets on the airfoil surface fully freeze on impact (rime ice conditions), or freeze partially (glaze ice conditions). In glaze type conditions, the remaining liquid water flows on the body surface as runback and eventually freezes downstream where surface temperature is below freezing point.

Ice shape computation relies on mass and energy conservation. The so-called Messinger Model [START_REF] Messinger | Equilibrium temperature of an unheated icing surface as a function of air speed[END_REF] predicts in each control volume located on the body surface the ice accumulation, the liquid water runback mass rate and the surface temperature. The conservation equations are solved in a serial way starting from the stagnation point to the trailing edge on the upper and the lower part of the airfoil. Inflow liquid water runback in the adjacent control volume is equalized with outflow liquid water runback of current control volume. This model is currently used in many ice accretion codes such NASA LEWICE [START_REF] Wright | Users manual for the improved nasa lewis ice accretion code lewice 1[END_REF],LEWICE3D, CANICE2D-NS, ONERA, TRAJICE(DERA),CIRA.

According to our knowledge, only two PDE-based approaches for ice thickness computation were developed. The PDE-based formulation developed by Myers [START_REF] Myers | Extension to the messinger model for aircraft icing[END_REF] takes into account the water film height as an additional property to the classical Messinger runback model. In the PDE model developed in FENSAP-ICE [START_REF] Beaugendre | FENSAP-ICE's three-dimensional in-flight ice accretion module: ICE3d[END_REF], the velocity of the runback water flowing on the surface is taken into account as well, allowing three dimensional icing computation in an unsteady fashion.

A non-PDE extension of the classic Messinger runback scheme to three dimension was proposed by [START_REF] Zhu | 3d ICE ACCRETION SIMULATION FOR COM-PLEX CONFIGURATION BASING ON IMPROVED MESSINGER MODEL[END_REF]. In that formulation liquid water film runback direction is mainly determined by the friction direction.

Comparison of PDE-based formulation with classic quasi-steady Messinger approach on experimental wind tunnel test in glaze ice conditions do not showcase dramatic improvements for either one. Physics of glaze ice accretion remains extremely complex and icing codes nowadays are only able to obtain qualitative results in such icing conditions.

In this paper we present the icing code developed in the structured compressible and multi-block Navier-Stokes solver NSMB [START_REF] Vos | Recent advances in aerodynamics inside the nsmb (navier-stokes multi-block) consortium[END_REF]. Supercooled water droplets field and impingement areas are obtained by means of an Eulerian approach solved implicitly using solvers BiGSTAB or SIP. The ice thickness is calculated based on an iterative Messinger model upgraded by an improved water runback scheme for threedimensional and multi-block ice accretion simulations. Mesh deformation is used to track the ice/air interface through time for multi-steps calculations.

The whole process is solved on the same mesh used for the fluid computation and parallelized with Message Passing Interface for complex multi-block configurations. Global overview of module interactions in NSMB is shown on Figure 1.

II Numerical method and governing equations A Modeling of Droplet Field

An Eulerian approach is developed for modeling droplet trajectories and impingement on a solid structure. The Eulerian model solves a droplet velocity field and a volume fraction distribution in the whole computational domain. Droplets are considered sphericals, with constant diameter and not to be subject to deformation or breaking. The droplet phase is sufficiently diluted to neglect interactions between droplets. Temperature of the droplet phase is set to be constant and equal to free-stream temperature, heat transfer with the surrounding air is neglected. The only forces acting on droplets are drag, gravity and buoyancy.

The governing equations for the conservation of mass and momentum of the droplets are written as follows [START_REF] Bourgault | A finite element method study of eulerian droplets impingement models[END_REF] :

{ ∂α ∂t + ∇.(αu) = 0 ∂αu ∂t + ∇. (αu ⊗ u) = α CDRe d 24K (u a -u) + α ( 1 -ρa ρ ) 1 F r 2 g (1)
α is defined as the non-dimensionalized volume fraction of water and u the non-dimensionalized velocity field of droplets. u a is the non-dimensionalized velocity of air, ρ the density of water, ρ a the density of air, g the gravity vector. F r = U ∞ /sqrt(Lg) is the Froude number, U ∞ the freestream speed of air, L the characteristic length (typically the airfoil chord length), K = ρdU ∞ /18Lµ is an inertia parameter and the dynamic viscosity of air, d is the median diameter of the droplets. The first term on the right-hand-side of the momentum equation accounts for the drag acting on the droplet or particle based on low-Reynolds number behaviour for spheres. The droplets Reynolds number ( Re d ) is defined based on the slip velocity between the air and droplet and the droplet diameter.

Re d = ρd|u a -u|U ∞ µ (2)
The drag coefficient of the droplets C d is given by :

     C d = 0.44 Re d , if Re d > 1000 C d = 24 Re d , if Re d < 0.1 C d = 24×(1+0.15×Re 0.657 d ) Re d (3)
For modeling the capture of incoming droplets on body surface, boundary conditions evolves during the iterative resolution such that Neumann boundary conditions are applied when the incoming droplet flux is positive or Dirichlet boundary conditions otherwise.

{ u.n < 0, Impingement, Neumann BC :

∂u ∂n = 0, ∂α ∂n = 0 u.n > 0, No Impingement, Dirichlet BC : u = 0, α = 0 (4)
The freestream values of droplet velocity and volume fraction are imposed as far field boundary conditions. Freestream volume fraction is set equal to α ∞ = 1 and droplet velocity vector equals to air freestream velocities.

The ability of the configuration to capture incoming droplets on solid walls is defined by the collection efficiency β, employing the non-dimensionalized formulation for droplet velocities and volume fraction defined above, the impingement efficiency factor writes :

β = αu.n (5)
A multi-block strategy paralelized with MPI is implemented for exchange of droplet state-vector values at block connectivity boundaries. The equations are treated as a generic transport problem for structured grids and recast in a finite volume form. A first order Upwind differencing Scheme (UDS) is used to evaluate the convective fluxes at interfaces. A blending Central Scheme/UDS scheme was implemented as well but causing convergence stability issues. A First order forward Euler scheme is employ to discretize the time derivative term allowing implicit resolution of the equations.

After residual linearisation, the following matrix system is obtained :

( V ∆t + ∂R ∂u ) ( u n+1 -u n ) ∆t = -R n (6)
Explicit terms are moved to the residual in the right-hand-side of the equation. Implicit boundary conditions are applied depending the droplet flux direction and system is solved. For memory requirements in case of large configurations, the sparse matrix obtained is inverted by iterative methods. Two algorithms have been implemented with success : the Stone's Method (Strongly Implicit Procedure) based on an incomplete LU decomposition and the biconjugate gradients stabilized method (CGSTAB). Convergence behaviour for both algorithms is almost identical but shows a slight improvement on CPU-time in favor of the CGSTAB method.

B Modeling of Ice Thickness

Ice thickness computation is based on mass conservation and energy conservation as stated by laws of thermodynamics. The so-called classic Messinger model solves the mass and energy conservation in control volumes located on surface body. The classic runback model consider inflow mass rate of liquid water in adjacent control volume equal to the outflow mass rate of runback liquid water film in the current control volume. This model need to be modified for three dimensional icing cases and multiblock grids. A modified iterative Messinger model is presented in this section.

The energy and mass flow rates conservation in control volume writes :

{ ṁimpinging + ∑ ṁinflow -ṁice -ṁevap - ∑ ṁoutflow = 0 qimpinging + ∑ qinflow -qevap -qice - ∑ qoutflow -qconv = 0 (7)
ṁinflow and ṁoutflow are respectively the inflow and the outflow mass rates of runback liquid water in the control volume. Collection efficiency β permits the calculation of the impinging mass rate of liquid water ṁimpinging . The evaporating mass rate ṁevap is determined by a parametric model. ṁice is the resultant mass rate of ice in the control volume.

qinflow and qoutflow are the energy rates of the runback liquid water film entering and leaving the control volume. The energy rate of impacting droplets is defined by qimpinging . qevap is the evaporating energy rate term. qice is the energy per second lost by solidification of the liquid water. The main contribution in the energy balance comes from the convection term qconv . The heat transfer coefficient necessary for the calculation of the convection term is obtained using Kays and Crawford [START_REF] Kays | Convective heat and mass transfer[END_REF] correlation. Conductive and radiative terms are neglected as they account for low contribution in the energy balance equation [START_REF] Wright | Users manual for the improved nasa lewis ice accretion code lewice 1[END_REF].

Freezing fraction f = ṁice ṁimpinging+ ṁinflow and surface temperature are obtained by means of Newton's iterations. A complete derivation and resolution method of equations is described in [START_REF] Wright | Users manual for the improved nasa lewis ice accretion code lewice 1[END_REF].

The extension to three-dimensional icing configurations requires a modification of the runback model to distinguish amount of liquid water flowing out in north-south and west-east directions. The model implemented in the present study is based on [START_REF] Zhu | 3d ICE ACCRETION SIMULATION FOR COM-PLEX CONFIGURATION BASING ON IMPROVED MESSINGER MODEL[END_REF] and consider that the liquid water film flows in the friction direction. The amount of liquid water flowing out towards control volumes located in west-east and north-south directions is determined by the following equations, subscript P denotes the current control volume :

Water runback in west/east direction :

ṁoutflow,we = C f , we |C f , we| + |C f , sn| . ∑ ṁinflow,P (8)
Water runback in south/north direction :

ṁoutflow,sn = C f , sn |C f , we| + |C f , sn| . ∑ ṁinflow,P (9)
The resolution procedure treats initially all control volumes as stagnation point by imposing a guess inflow liquid water mass rate of ṁinflow = 0. The liquid water film properties are then solved iteratively by the Messinger model and the runback scheme. In such formulation the treatment of stagnation points is straightforward and should permit to handle multiple stagnation points problems that are a common feature of glaze type icing. The iterative procedure permits multiblock and parallelization implementation and is an overall improvement to classic and serial ice accretion codes. Summary of the resolution procedure is the following : (1) initialize ṁinflow = 0 in all control volumes , (2) compute Messinger model, (3) Compute Runback distribution and update ṁinflow in control volumes, (4) repeat steps (2) and (3) until convergence.

By knowing the mass rate of ice in each control volume, the computation of the volume of ice V ice is straightforward : with ρ ice the density of ice, dt the icing timestep and ∆s the surface of the control volume.

V ice = ṁice dt ρ ice ∆s (10)
The ice is considered growing in the normal direction of the surface of the control volume.

Computed ice thickness δ are provided to the ALE mesh deformation algorithm available in NSMB that generates a new mesh.

III Results

Simulations are performed on a two-dimensional NACA23012 & NACA00012 aifoils, an Onera M6 swept wing and a Wing/Body DLR-F6 configuration. The flow computation is accomplished implicitly using LU-SGS method. Spatial discretization is achieved using second order or fourth order central scheme with matrix artificial dissipation. Spalart-Allmaras is used as turbulence model. Table 1 summarizes the icing cases studied.

A Two-Dimensional Configurations

A.1 Rime Ice Accretion Simulation

Rime icing is simulated on a two-dimensional NACA23012 airfoil and compared with experiments of Broeren et al. [START_REF] Broeren | Effect of high-fidelity ice-accretion simulations on full-scale airfoil performance[END_REF] The choice of these particular validation cases was motivated by comparing our numerical results with recent high fidelity experimental ice shapes.

The O-grid generated by NS-GRID software (citation) is shown in Fig. 2. The NACA23012 mesh has a chord of c = 1.828m and contains 512 × 256 cells divided into two blocks. The parameters of the simulations ( Cases 1 & 2) are summarized on Table 1. Visualization of Liquid Water Content distribution for case 2 is shown in Fig. 3. Shadow zones observed evince dried areas where liquid water content is nearly zero. The Collection efficiency curves are shown in Fig. 4. Ice shapes obtained after a single icing step are compared with experimental results in Figs. 5 and6. The simulated ice shapes are in good accordance with experimental data. Taking into account a corrected angle of attack and a droplet size distribution to simulate the spray is expected to provide more precise results on the upper impingement limit where a slight underestimation of the ice thickness is obtained.

In addtion, flow is simulated on clean and automatically deformed mesh for case 2 at an angle of attack of AOA = 11.9 • , a Mach Number of Ma = 0.20 and a Reynolds number of Re = 16.1×10 6 . Pressure Coefficients are compared for both clean and iced configurations with experimental data in Fig. 7. Pressure coefficient is highly increased at location where ice is accreted resulting in leading edge stall. A separation bubble occurs in the vicinity of the accreted ice as well as a flow reattachment downstream. Comparison of results from Fig. 7 showcases the validity of present code to predict aerodynamic behavior resulting from rime icing.

A.2 Glaze Ice Accretion Simulation

Glaze icing is simulated on a two-dimensional NACA00012 airfoil and compared with LEWICE [START_REF]LEWICE Validation CD-ROM[END_REF]. Numerical Results of LEWICE were obtained by means of a panel method to simulate the flow, a Lagrangian formulation for droplet tracking, a boundary layer method to evaluate the heat transfer coefficient and a classic Messinger model to compute the ice thickness.

The O-grid mesh generated by NS-GRID has a chord of c = 0.5334m and contains 256×128 cells.

The sand grain roughness used in the computation is set to k s = 5mm. Results of impingement efficiency distribution along the airfoil surface is shown on Fig. 8. Present results showcase a higher value of the lower impingement limits and a slight increase of peak value of collection efficiency. These observations are in accordance with results of several studies of the Eulerian droplet model [START_REF] Jung | Numerical modeling for eulerian droplet impingement in supercooled large droplet conditions[END_REF][START_REF] Bourgault | A finite element method study of eulerian droplets impingement models[END_REF]. The calculated glaze-type ice shape is compared with LEWICE in Fig. 9 after a single icing step. A Horn forms on the upper part of the airfoil in accordance with results of LEWICE. However a discrepancy is observed on the lower part of the ice shape where a secondary horn start to form in LEWICE results while the ice thickness remains smooth in present results.This difference likely account for the different Messinger-based models used in both simulations. Test case was repeated on a NACA00012 extruded mesh. Identical results were obtained along wing span validating the three-dimensional implementation.

B Onera M6 Swept Wing

Droplet transport equations are solved on an Onera M6 swept Wing. Impingement distribution is compared with results from Sang et al. [START_REF] Sang | Numerical simulation of icing effect and ice accretion on threedimensional configurations[END_REF]. The wing has a span of 1m and a mean aerodynamic chord of MAC=0.53m. The structured mesh used contains approximately 884000 cells and is shown in Fig. 10. The grid is available on CFL3D validation website http://cfl3d. larc.nasa.gov/Cfl3dv6/cfl3dv6_testcases.html. Collection efficiency distribution on airfoil is shown in Fig. 11. Impingement efficiency gets higher as the chord gets smaller which results in spare accreted ice close to the tip. This phenomena is a common result for swept wings [START_REF] Vargas | Ice accretions on a swept glc-305 airfoil[END_REF]. Impingement curves are plotted on sections A minor deviation of results is obtained near to the stagnation point. The origin of this discrepancy is unclear but similar collection efficiency behavior with low droplet diameters using Eulerian droplet model was obtained by Kim et al. on a two-dimensional MS-317 airfoil [START_REF] Kim | Ice accretion modeling using an eulerian approach for droplet impingement[END_REF].

C DLR-F6 Wing Body Configuration

The last test case concerns the simulation of rime icing on a DLR-F6 wing/body configuration. The grid is available from the DPW3 Drag Prediction Workshop website http://aaac. larc.nasa.gov/tsab/cfdlarc/aiaa-dpw/Workshop3/workshop3.html. The structured mesh contains over 8.8 × 10 6 divided in 152 blocks. The mean aerodynamic chord is MAC= 0.14m. Simulated rime ice accretions are shown in Fig. 14. A global overview of icing on the geometry is displayed as well as the ice shapes on Sections A,B,C & D. Results shows that ice only accumulates on the nose, the cockpit and the leading edge of the wing. In accordance with collection efficiency distribution, ice thickness gets higher on the leading edge near the tip indicating that anti-icing equipment must be particularly effective in that location. At the opposite, no ice accretion occurs close to the wing root leading edge. The deformed mesh at the tip vincinity is shown in Fig. 15. NSMB is validated with experimental and numerical data on NACA23012 and NACA0012 twodimensional configurations. Results showcase global good agreement to predict droplet trajectories and ice shape. Pressure coefficient on iced airfoil is compared with experimental data, similar aerodynamics behavior is obtained validating the methodology.

For three-dimensional validation, collection efficiency distribution is compared with numerical results on an ONERA M6 swept wing with good agreement.

Rime icing is computed of a wing/body DLR-F6, results are qualitatively in accordance with study on a wing/body High-Lift configuration.

Limitations of structured automatic mesh movement arises in complex glaze-type ice shapes. In future work, efforts will be made on this particular issue and alternative methods will be investigated.

  Droplet and ice thickness modules are presented in Sect. II and validated on several test cases including NACA23012 and NACA00012 airfoils in Sect. III.A and an ONERA M6 wing in Sect. III.B. Collection efficiency and rime ice accretion are computed on a Wing / Body DLR-F6 configuration in Sect. III.C.

Figure 1 :

 1 Figure 1: Module Interactions in NSMB

Figure 2 :Figure 3 :

 23 Figure 2: NACA23012 O-grid, 512 × 258 cells

Figure 4 :

 4 Figure 4: Case 1 & Case 2 -NACA23012 -Collection Efficiency distribution

Figure 5 :Figure 6 :

 56 Figure 5: Case 1 -NACA23012 -Comparison of ice shape with experimental data of Broeren et al.

Figure 7 :

 7 Figure 7: Case 2 -NACA23012 -Comparison of pressure coefficient with experimental data for clean and iced airfoils -Ma= 0.20, Re = 16.1 × 10 6 , AOA = 11.9 • .

Figure 8 :Figure 9 :

 89 Figure 8: Case 3 -NACA00012 -Comparison of collection efficiency with LEWICE.

  Flow is computed at an angle of attack of AOA = 6 • , a freestream velocity of U ∞ = 50m/s and a temperature of T = 263.15K. Icing parameters are droplet diameters of d = 20µm, a liquid water content of LW C = 1g.m 3 .

Figure 10 : 6 Figure 11 : 2 -

 106112 Figure 10: Case 4 -ONERA M6 -Visualization of the grid.

Figure 12 :

 12 Figure 12: Case 4 -ONERA M6 -Comparison of Collection Efficiency with results from Sang et al. on sections 1 & 2.

Figure 13 :

 13 Figure 13: Case 5 -DLR-F6 -Visualization of collection efficiency on the fuselage and plot of impingement curves on sections A,B,C & D.

Figure 14 :

 14 Figure 14: Case 5 -DLR-F6 -Global View of ice accretion, ice shapes on sections A,B,C & D.

Figure 15 :

 15 Figure 15: Case 5 -DLR-F6 -View of the deformed mesh in the vicinity of wing tip leading edge.

Table 1 :

 1 Icing Cases Studied

	Case	Configuration	AOA	Mach	Reynolds	T	d	LW C∞	Icing time
	no.		[ • ]		[×10 6 ]	[K]	[µm]	[g/m 3 ]	[s]
	1	NACA23012	2.0	0.24	12.60	247.85	30	0.55	600
	2	NACA23012	2.0	0.32	16.23	252.45	15	0.30	1200
	3	NACA00012	4.0	0.32	4.44	262.04	20	0.55	65
	4	ONERA M6	6.0	0.15	2.167	263.15	20	1.0	//
	5	DLR-F6	0.0	0.19	3.172	249.15 14.5	0.55	120

  1 & 2 respectively at 60% and 90% of span length and compared with Sang et al. in Fig. 12. Results from Sang et al. were obtained by means of an Euler Flow Solver and a Lagrangian method for droplets tracking and do not include impingement limits of collection efficiency. Comparison with our results showcases a global good agreement.

V Acknowledgments

The authors would like to thank Compute Canada for computational ressources, Kazem Hasanzadeh for grids generation and Thibaut Deloze for assistance and discussions.


Dianailys Nuñez-Reyes

Jean-Christophe Loison

Kevin M Hickson

Michel Dobrijevic

A low temperature investigation of the N( 2 D) + CH4, C2H6 and C3H8 reactions

The gas-phase reactions between metastable nitrogen atoms, N( 2 D) and saturated hydrocarbons CH4, C2H6 and C3H8 have been investigated using a supersonic flow reactor over the 296 -75 K temperature range. N( 2 D) was generated as a product of the C( 3 P) + NO ® N( 2 D) + CO reaction, with C( 3 P) atoms created in-situ by pulsed laser photolysis of CBr4.

The kinetics of N( 2 D) loss were followed by vacuum ultraviolet laser induced fluorescence.

The measured rate constants for the N( 2 D) + CH4 reaction are in good agreement with earlier work and extend the available kinetic data for this process down to 127 K. The measured rate constants for the N( 2 D) + C2H6 and N( 2 D) + C3H6 reactions are in reasonable agreement with previous work at room temperature and extend the available kinetic data for these processes down to 75 K. The rate constants for all three reactions decrease as the temperature falls, indicating the presence of activation barriers for all three processes. While the recommended values for the low temperature rate constants of the N( 2 D) + CH4 reaction are close to those reported here, the previous recommendations for the other saturated hydrocarbon reactions significantly overestimate the rate constants for these processes. The effects of the new rate constants on a coupled ion-neutral photochemical model of Titan's atmosphere are discussed.

Introduction

Nitrogen is the fifth most abundant element in the Universe behind hydrogen, helium, carbon and oxygen. 1 The reactions of atomic nitrogen are important processes in a wide range of astrochemical environments. In the interstellar medium, N-atoms in their ground electronic 4 S state participate in a series of reactions with other nitrogen containing radicals, effectively controlling the transformation from N to N2 in dense interstellar clouds. 2,3 In these regions, the specific form of the major elemental nitrogen reservoir (N or N2) is critically important for the overall complexity of the molecules present, considering the highly unreactive nature of molecular nitrogen. Excited electronic states of atomic nitrogen, N( 2 D) and N( 2 P), play little or no role in the chemistry of the gas-phase interstellar medium as the timescale for radiative relaxation is shorter than the collision timescale. Nevertheless, in contrast to cold interstellar environments, excited state nitrogen atoms participate in the gas-phase chemistry of planetary atmospheres, where large fluxes of energetic photons and high densities lead to numerous collisions with the various atmospheric components. In the Earth's atmosphere, atomic nitrogen is produced through the short wavelength photodissociation of N2 and through N2 collisions with photoelectrons in the upper atmosphere. 4 Indeed, the reactions of N( 4 S) and N( 2 D) with O2 are both thought to be major sources of NO in the thermosphere, while subsequent reactions of N( 4 S) and N( 2 D) with NO 5,6 limit the overall odd nitrogen concentration at these altitudes. In Titan's upper atmosphere, both N( 4 S) and N( 2 D) atoms are considered to be important reactive species for the formation of nitrogen-bearing compounds.

While N( 4 S) atoms react with radicals such as CH3, 7 to form nitrogen containing hydrocarbon molecules such as H2CN (and ultimately HCN), 8 these ground state atoms are largely unreactive with the most abundant closed-shell hydrocarbon molecules present in Titan's atmosphere, CH4, C2H2, C2H4 and C2H6. Instead, the more reactive metastable N( 2 D) atoms react with these molecules to form a range of amines, imines and nitriles 9 (N( 2 P) is much less reactive than N( 2 D), 10 so that the fate of these atoms is quenching to the 2 D state through collisions with N2, rather than reaction). Although the major reaction channels for the N( 2 D) + CH4, [11][12][13][START_REF] Rosi | Computational Science and Its Applications -ICCSA 2012[END_REF][START_REF] Rosi | Computational Science and Its Applications -ICCSA 2013[END_REF] N( 2 D) + C2H6, 16 N( 2 D) + C2H4 17,18 and N( 2 D) + C2H2 [19][20][21] reactions have been identified theoretically and through crossed molecular beam experiments, there are few kinetics measurements of these processes and none at temperatures pertinent to Titan's atmosphere. For the N( 2 D) + CH4 reaction, the derived room temperature rate constants [22][23][24][25] fall in the range (3.0 -5.4) × 10 -12 cm 3 s -1 . Black et al. 22 generated N( 2 D) atoms in their reactor through N2O photolysis at 147 nm and 123.6 nm. N( 2 D) removal in the presence of coreagent CH4 was followed using the chemiluminescent emission from NO(B 2 P®X 2 P) produced by the N( 2 D) + N2O reaction. In their study, Fell et al. 23 used a fast flow system where N( 2 D) atoms were produced through the microwave discharge of N2 diluted in He. The N( 2 D) concentration was monitored through electron spin resonance spectroscopy. They measured rate constants for N( 2 D) removal for numerous coreagent species including CH4, C2H6 and C3H8. Umemoto et al. 24 created N( 2 D) atoms in their experiments through the multiphoton dissociation of NO at 275.3 nm. These authors detected N( 2 D) through twophoton excitation around 269 nm via the 2p 3 2 D5/2 -2p 2 3p 2 S transition. Rapid relaxation to the 2p 2 3s 2 P state followed by VUV fluorescence at 149 nm allowed these authors to determine rate constants for the reactions of N( 2 D) with various species at room temperature including CH4, C2H6 and C3H8. Takayanagi et al. 25 used the pulsed radiolysis -resonance absorption method to perform kinetic measurements over the 292 K-223 K range, reporting a positive temperature dependence for the N( 2 D) + CH4 reaction, with the rate constant decreasing to a value of (2.3 ± 0.2) × 10 -12 cm 3 s -1 at 223 K. The Arrhenius parameters of A = 4.8 × 10 -12 cm 3 s -1 and E/R = 750 recommended by Herron 10 were based on the temperature dependence of the Takayanagi et al. 25 study and the earlier room temperature work (with a correction factor applied to the Fell et al. rate constants -see Herron 10 for further justification). For the N( 2 D) + C2H6 and N( 2 D) + C3H8 reactions, kinetic studies have only been performed at room temperature. Fell et al. 23 and Umemoto et al. 24 obtained rate constants of 2.7 × 10 -11 cm 3 s -1 and 2.1 × 10 -11 cm 3 s -1 for the N( 2 D) + C2H6 reaction and 4.6 × 10 -11 cm 3 s -1 and 3.1 × 10 -11 cm 3 s -1 for the N( 2 D) + C3H8 reaction respectively. The recommended rate constants 10 for these two reactions are 1.9 × 10 -11 cm 3 s -1 (N( 2 D) + C2H6) and 2.9 × 10 -11 cm 3 s -1 (N( 2 D) + C3H8) at 298 K, considering the correction factor of 0.6 that was systematically applied to the Fell et al. results. Given the lack of temperature dependent data, current models of Titan's atmosphere 26,27 consider temperature independent rate constants for these two processes, equal to the recommended room temperature values.

This study presents kinetic measurements for the N( The experiments performed in this study employed an existing continuous supersonic flow reactor whose main features have been described previously. [28][29][30] Further modifications to the apparatus allowing the generation of narrowband tunable radiation in the vacuum ultraviolet wavelength range for the detection of atomic radicals in ground and excited electronic states are reported in later work (C( 3 P), [31][32][33] H( 2 S) 31,[33][34][35][36][37][38][39][40][41] and D( 2 S), 33,42 O( 1 D), 38,41,[43][44][45][46][47] and N( 2 D)). 6 During this work, both Ar and N2 were used as carrier gases as they result in only slow quenching of N( 2 D) atoms. 48,49 The properties of the Laval nozzles employed in this work are summarized in Table 1. The ratio [N( 2 D) + N( 4 S)]/[O( 3 P)] has been estimated as 1.5 ± 0.3 at room temperature 50 while previous theoretical work predicts that the branching ratio for channel (1a) should decrease by approximately 50 % between 300 K and 50 K. 51 Ground state C( 3 P) atoms were produced by the pulsed multiphoton dissociation of carbon tetrabromide, CBr4, at 266 nm with pulse energies around 26 mJ. C( 1 D) atoms were also produced during the photolysis process. 31 Possible interferences from secondary reactions are discussed below. To introduce these molecules into the system, a small flow of the carrier gas was diverted into a vessel containing solid CBr4 at a known fixed (room) temperature and pressure before rejoining the main flow upstream of the Laval nozzle reservoir. The CBr4 concentration in the flow was estimated to be less than 3.2 ´ 10 13 cm -3 as determined by the product of the ratio of the CBr4 The PMT was protected from reactive gases within the chamber by a LiF window while the zone separating the PMT from the reactor (containing a LiF lens to focus the VUV light onto the photocathode) was continuously evacuated by a dry pump. As signal levels were small compared with our earlier studies of atomic carbon and oxygen reactions, spurious light sources (residual UV light, scattered VUV light and window fluorescence) were reduced to a minimum, while it was also necessary to connect a preamplifier module to the output of the PMT before signal processing and integration by a boxcar system. Although this provided a substantial improvement in signal-to-noise levels, scattered light and/or window fluorescence from the photolysis laser which could not be entirely removed was also detected. This resulted in saturation of amplified signal for the first 15 microseconds following the photolysis pulse, preventing us from recording VUV LIF signals in this time period throughout the experiments. A delay generator was used to scan the time between photolysis and probe lasers to obtain the temporal profile of N( 2 D) atoms. 30 points were recorded at each time point with at least 70 time points recorded for each kinetic profile. Points recorded at negative time delays (that is with the probe laser pulse occurring before the photolysis laser one) allowed the baseline level to be established. Gases (Linde Ar 99.999%, Xe 99.999%, CH4 99.9995%, C2H6 99.5%, C3H8 99.5%, Air Liquide N2 99.999%, NO 99.9%) were used without further purification. Gas cylinders were attached to calibrated mass-flow controllers, allowing the various carrier gas and reagent flows to be regulated precisely.

Potential secondary chemistry

Although several secondary reactions could be occurring in the reactor during these experiments, these are not expected to interfere with the present measurements to any significant extent considering that N( 2 D) atoms are detected directly. A discussion of the influence of potential secondary reactions has already been provided by Nuñez-Reyes & Hickson 6 for the simpler case of the N( 2 D) + NO reaction. In the present case, additional secondary processes arising from the reactions of photolysis products C( 3 P) and C( 1 D), quenching product N( 4 S), as well as reaction products CN and O( 1 D) with the saturated hydrocarbon reagents also need to be considered. Both C( 3 P) and N( 4 S) atoms are unreactive with CH4, C2H6 and C3H8 at room temperature and below, so that these reactions do not need to be considered further. Although C( 1 D) will react with all three of these reagents, leading to atomic products such as H( 2 S) and hydrocarbon radicals such as CH and CH3, 39, 52 these species will be present in low concentrations and none will yield N( 2 D) as a product through reaction with NO, the only nitrogen containing excess reagent present in the chamber.

Similarly, although the CN and O( 1 D) products of the C + NO reaction react readily with saturated hydrocarbon molecules, 38,41,53 none of these processes are thought to lead to the formation of N( 2 D) atoms.

Results and Discussion

The pseudo-first-order approximation was applied in these experiments, by maintaining both NO and the hydrocarbon reagent CxH2x+2 in excess with respect to both C( 3 P) and N( 2 D).

Under these conditions, the temporal profile of N( 2 D) atoms, 𝐼 1( 2 3) , can be described by the following expression

𝐼 1( 2 3) = C(exp(-𝑘 ; < 𝑡) -exp(-𝑘 > < 𝑡)) (2)
where C is a constant, 𝑘 ; < is the pseudo-first-order rate constant for N(

to describe the N( 2 D) signal as a function of time. In this respect, the starting points of individual fits were carefully chosen to avoid the rising part of the temporal profile. Figure 1 shows two profiles of the N( 2 D) fluorescence intensity as a function of time recorded at 296 K in the presence and absence of C3H8. (solid black circle) Black et al. 22 ; (open red circle) Fell et al. 23 ; (solid purple circle) Umemoto et al. 24 ; Takayanagi et al. 25 The solid blue lines represent best fits to the present experimental data. The N( 2 D) + C2H6 reaction: (Green solid squares) this work; (black solid square) Fell et al. 23 ; (blue open square) Umemoto et al. 24 The N( 2 D) + C3H8 reaction: (Cyan solid triangles) this work; (red solid triangle) Fell et al. 23 ; (black solid triangle) Umemoto et al. 24 The dashed lines represent the previously recommended values for these reactions. Error bars on the present data represent the combined statistical (1s) and systematic errors (estimated to be 20 % of the nominal value of the rate constant except for measurements of the N( 2 D) + CH4 reaction at 177 K and 127 K where it was estimated to be 50 %).

The N( 2 D) + CH4 reaction At room temperature, experiments were performed with two different NO concentrations (2.2 × 10 14 cm -3 and 4.3 × 10 14 cm -3 ) to test for the influence of potential secondary reactions, yielding similar values for 𝑘 [( 2 \)?$] ' (296 K) = (2.1 ± 0.4) × 10 -12 cm 3 s -1 and (2.5 ± 0.5) × 10 -12 cm 3 s -1 . These values are slightly lower than those obtained by previous studies [22][23][24][25] with values in the range (3.0 -5.4) × 10 -12 cm 3 s -1 . The measured rate constants for the N( 2 D) + CH4 reaction below room temperature are slower than 10 -12 cm -3 s -1 ; a value which is normally considered as the lower limit for rate constant measurements through the Laval nozzle method (assuming that the minor reagent decay can be followed over at least two half-lives). If we consider the measurements performed at 177 K as an example, the nozzle used has a characteristic hydrodynamic time of ~250 microseconds (that is, the time over which the nominal characteristics of the supersonic flow are valid). In this case, the maximum [CH4] used in these experiments was 5.6 × 10 15 cm -3 which constitutes approximately 6 % of the total flow, and therefore approaching the upper limit for the added coreagent concentration before significant deviations from the nominal flow conditions occur (N2 carrier gas density = 9.4 × 10 16 cm -3 ). The measured rate constant for the N( 2 D) + CH4 reaction at 177 K is 4.7 × 10 -13 cm 3 s -1 allowing us to calculate a half-life time, 2), it is important to realize that the derived second-order rate constants at 177 K and at 127 K are characterized by substantially larger uncertainties than usually associated with this type of measurement. In this respect, we have increased our estimate of the systematic errors on these two measurements to 50 % of the nominal value, although an upper limit value of 3.2 × 10 -13 cm 3 s -1 could also be considered as appropriate for the measured rate constant at 127 K. Nevertheless, the measured rate constants are in reasonably good agreement with the recommended ones. 10 Indeed, the recommended Arrhenius parameters A = 4.8 × 10 -11 cm 3 s -1 and (Ea/R) = 750 K yield 𝑘 [( 2 \)?$] ' (177 K) = 6.9 × 10 -13 cm 3 s -1 and 𝑘 [( 2 \)?$] ' (127 K) = 1.3 × 10 -13 cm 3 s -1 which are close to the values obtained in this work (𝑘 [( 2 \)?$] ' (177 K) = (4.7 ± 2.4) × 10 -13 cm 3 s -1 and 𝑘 [( 2 \)?$] ' (127 K) = (1.7 ± 1.5) × 10 -13 cm 3 s -1 ). An Arrhenius fit to the present data yields the following Arrhenius parameters; A = 1.93 × 10 -11 cm 3 s -1 and (Ea/R) = 638 K. Given the slow rate constants at 177 K and 127 K, no attempt was made to measure rate constants at lower temperature.

The N( 2 D) + C2H6 reaction In a similar manner to the N( 2 D) + CH4 reaction, experiments were performed with two different NO concentrations (4.4 × 10 14 cm -3 and 6.6 × 10 14 cm -3 )

to test for the influence of potential secondary reactions, yielding similar values for

𝑘 [( 2 \)?$ 2 ] ^(296 K) = (1.2 ± 0.
3) × 10 -11 cm 3 s -1 and (1.3 ± 0.3) × 10 -11 cm 3 s -1 . These values are slightly lower than previous work by Fell et al. 23 and Umemoto et al. 24 which yielded values of (1.6 ± 0.6) × 10 -11 cm 3 s -1 (corrected value -see Herron 10 for more detail) and (2.1 ± 0.2) × 10 -11 cm 3 s -1 respectively. The rate constant is seen to decrease as the temperature falls, reaching a value of (1.1 ± 0.3) × 10 -12 cm 3 s -1 at 75 K. An Arrhenius fit to the present data does not provide an adequate description of the rate constant over the entire temperature range, so we choose instead to use a modified Arrhenius type expression (k = A × (T/300) n × exp(-Ea/RT)) where A = 3.58 × 10 -12 cm 3 s -1 , n = 4.45 and (Ea/R) = -379 K. These parameters are valid over the 75 -296 K temperature range. The deviation from Arrhenius behavior at low temperature might be an indication that H atom tunneling plays a small but non-negligible role in the N( 2 D) + C2H6 reaction. Interestingly, photochemical models of Titan's atmosphere 26,27 use the recommended room temperature rate constant 10 𝑘 [( 2 \)?$ 2 ] ^ = 1.9 × 10 -11 cm 3 s -1 as a temperature independent value to describe the reactivity of the N( 2 D)

+ C2H6 reaction at lower temperature. At temperatures relevant to Titan's atmosphere (taking 150 K as a representative temperature) the fit to the present data yields a rate constant of 2.1 × 10 -12 cm 3 s -1 ; almost an order of magnitude lower than the currently recommended value.

Following the joint experimental and theoretical study of this reaction by Balucani et al. 16 the major products are expected to be methanimine (CH2NH) and CH3 with small contributions from ethanimine (CH3CH2NH) and H atoms, NH and C2H5, as well as other products at the 1 % level or less.

The N( 2 D) + C3H8 reaction

The measured rate constants at 296 K were found to be lower than the values obtained previously by Fell et al. 23 and Umemoto et al. 24 of (2.8 ± 0.8) × 10 -11 cm 3 s -1 (corrected value) and (3.1 ± 0.3) × 10 -11 cm 3 s -1 respectively. Given the discrepancy between the present and previous studies, further tests were performed to check for potential errors in the current work. In particular, as N( 2 D) atoms are formed through reaction (1a) rather than prompt formation through direct photolysis, it is possible that N( 2 D) formation and N( 2 D) loss through the N( 2 D) + C3H8 reaction were occurring over the same time period.

Indeed, such an effect is expected to lead to slower pseudo-first-order rate constants for the target reaction. During these test experiments, pseudo-first-order rate constants for the N( 2 D) + C3H8 reaction were measured with three different fixed NO concentrations (2.2, 4.2 and 6.3

× 10 14 cm -3 ) at room temperature. Assuming a rate constant of 1.5 × 10 -10 cm 3 s -1 for reaction

(1), 54 we obtain pseudo-first-order rate constants of 3.4, 6.3, 9.4 × 10 4 s -1 respectively thus allowing us to calculate the half-life times for N( 2 D) formation under these conditions. The delay times when only 13 % of the initial C( 3 P) remains to be converted to N( 2 D) (corresponding to three half-lives) are 64, 33 and 22 µs, indicating that data acquired at delays shorter than these values may not accurately reflect the kinetics of the N( 2 D) + C3H8 reaction.

For the lowest NO concentration, the data acquired after 64 µs was difficult to analyse given the weak signal levels, with a measured rate constant of (1.2 ± 0.2) × 10 -11 cm 3 s -1 . With the intermediate NO concentration, a similar measured rate constant of (1.2 ± 0.2) × 10 -11 cm 3 s -1 was obtained. Only the rate constant value obtained for the N( 2 D) + C3H8 reaction with [NO] = 6.3 × 10 14 cm -3 (yielding 𝑘 [( 2 \)?$ U ] _ (296 K) = (1.7 ± 0.4) × 10 -11 cm 3 s -1 respectively) was retained to derive the temperature dependent fitting parameters. The rate constant decreases as the temperature falls reaching a value of (3.2 ± 0.9) × 10 -12 cm 3 s -1 at 75 K. As the Arrhenius expression does not yield a good fit to the data over the entire temperature range, the modified Arrhenius type expression described above is used instead, where A = 7.36 × 10 -12 cm 3 s -1 , n = 2.99 and (Ea/R) = -252 K. These parameters are valid over the 75 -296 K temperature range. Current photochemical models of Titan's atmosphere 26,27 use the recommended room temperature rate constant 10 𝑘 [( 2 \)?$ U ] _ = 2.9 × 10 -11 cm 3 s -1 as a temperature independent value to describe the reactivity of the N( 2 D) + C3H8 reaction. At 150 K, the fit to the present data predicts a rate constant of 5.0 × 10 -12 cm 3 s -1 ; six times smaller than the currently recommended value.

Influence of the new rate constants on Titan's atmospheric chemistry

We wanted to test the effect of the new rate constants on current models of Titan's atmosphere. To do this, we employed the 1D-photochemical model described by Dobrijevic et al. 26 with more recent updates related to the chemistry of aromatic compounds, [START_REF] Loison | Icarus[END_REF] although a detailed description is not provided here (see Dobrijevic et al. 26 ). This model treats the coupled chemistry of neutral and ionic species (both cations and anions) to provide an accurate description of the chemistry occurring from the lower atmosphere all the way up to the ionosphere, including high resolution absorption cross-sections for an improved description of N2 photodissociation. Earlier models by these authors which did not include ions already updated the neutral chemistry of nitrogen bearing molecules. 8,9,56 Indeed, an uncertainty propagation study by Loison et al. 9 focusing on the chemistry of nitriles, amines and imines clearly identified the N( 2 D) + CH4 and N( 2 D) + C2H6 reactions as key reactions for future study in order to improve the accuracy of model predictions.

Here, we ran two different simulations to examine the effects of the measured rate constants on the species mole fractions as a function of altitude. The first employed the previously recommended rate constants and branching ratios for the N( 2 D) + CH4, C2H6 and C3H8 reactions to yield the standard profiles. The second used the new values as described by the Arrhenius parameters for the respective reactions as presented in section 3. As we do not measure the various product channels for these processes, the relative branching ratios remain unchanged. The differences produced by the two models for various species at 1000 km are summarized in Table 3. As CH3NH2 is formed through the NH2 + CH3 reaction, this species also decreases slightly in the new model due to the less efficient CH2NH2 + + e -dissociative recombination reaction where NH2 is a minor reaction product as described above. Similarly, the NH3 mole fraction also decreases as it forms mostly from reactions of NH2.

The

  photochemical model of Titan's atmosphere. The conclusions are outlined in Section 5.

  vapour pressure and the total pressure and the flow rate through the vessel (𝑃 $%& ' /𝑃 )*++*, ) × 𝐹 )*++*, . The method employed to follow N( 2 D) atoms in the present work was on-resonance pulsed laser induced fluorescence in the vacuum ultraviolet range (VUV LIF) via the 2s 2 2p 3 2 D° -2s 2 2p 2 ( 3 P)3d 2 F transition at 116.745 nm. Tunable radiation at and around this wavelength was produced in an identical manner to study of Nuñez-Reyes & Hickson 6 where a more detailed description of the procedure for its generation can be found. The VUV radiation was collimated by a MgF2 lens at the exit of the tripling cell and steered into the reactor through a 75 cm long sidearm containing baffles to reduce the flux of divergent residual UV radiation. The VUV beam overlapped with the cold supersonic flow at right angles. As the coreagents CH4, C2H6 and C3H8 all absorb strongly in the VUV region, this zone was continuously flushed with either N2 or Ar to minimize the attenuation of the probe beam. The solar blind photomultiplier tube detector (PMT) was placed at right angles to both the probe beam and the supersonic flow, collecting fluorescence from unreacted N( 2 D) atoms.

Figure 1 N

 1 Figure 1 N( 2 D) VUV LIF intensity as a function time recorded at 296 K. (Red solid circles) without C3H8; (blue-solid circles) [C3H8] = 2.3 × 10 15 cm -3 . The signal recorded in the presence of C3H8 has been multiplied by an arbitrary factor to allow a comparison of the two traces on the same intensity scale.

  Decays similar to those displayed in Figure1were recorded at several values of [CxH2x+2] whilst maintaining a fixed value for [NO]. In this way, any change in the observed decay rate of N( 2 D) atoms was entirely due to the change in [CxH2x+2]. Plots of the measured pseudofirst-order rate constant versus [CxH2x+2] thus yielded the second-order rate constant from weighted linear least-squares fits to the data. The decay curves recorded at high [CxH2x+2] were always considerably more scattered than the data recorded at low [CxH2x+2] due to absorption of the N( 2 D) VUV LIF emission by residual saturated hydrocarbon present in the reactor. Two example second-order plots are shown in Figure 2, for the reactions of N( 2 D) with CH4 (purple triangles) and C2H6 (red squares) at 177 K.

Figure 2 Figure 2

 22 Figure 2 Pseudo-first-order rate constant as a function of saturated hydrocarbon concentration recorded at 177 K. (Red solid squares) the N( 2 D) + C2H6 reaction; (purple solid triangles) the N( 2 D) + CH4 reaction. [NO] = 2.7 × 10 14 cm -3 . Solid lines represent weighted linear leastsquares fits to the data. The error bars on individual data points (1s) are derived from exponential fits to temporal profiles such as those displayed in Figure 1.

Figure 3 .

 3 Figure 3. Second order rate constants for the N( 2 D) + CxH2x+2 reactions as a function of temperature. The N( 2 D) + CH4 reaction: (Blue solid circles) this work; (purple open circles)

t1/ 2

 2 = (ln(2)/k¢) = 260 microseconds. Consequently, the loss of N( 2 D) due to the N( 2 D) + CH4 reaction alone over the entire duration of the flow at 177 K is significantly less than the loss due to the competing N( 2 D) + NO reaction (t1/2 = 29 microseconds). Although a clear progression of the pseudo-first-order rate constant as a function of [CH4] is obtained in these experiments (see Figure

5 Conclusions

 5 N( 2 D) + CH4 reaction remains the major loss process for N( 2 D) atoms in the atmosphere of Titan representing approximately 28 % of the total N( 2 D) removal, followed by the N( 2 D) + HCN(22 %) and N( 2 D) + C2H2 (13 %) reactions. Whereas the N( 2 D) + C2H6 reaction was the 9 th most important loss process for N( 2 D) (representing approximately 1% of its total removal), with the new rate constants it now represents only 0.1 % of N( 2 D) loss.While the new model results are similar to those of the standard model, it will be interesting to test the effects of the new and more realistic error bars for these reactions on the mole fraction profiles of a range of atmospheric species obtained through uncertainty propagation studies. Previous models (see Loison et al.[START_REF] Loison | Icarus[END_REF] for example) employed values for the uncertainty factor, F = 1.5 for the N( 2 D) + C2H6, C3H8 reactions. In light of the present experimental results, it is clear that these uncertainties were severely underestimated. While these findings illustrate the difficulty of attributing realistic uncertainties to reactions that have only been investigated at temperatures far from the relevant range, they also reinforce the need for kinetic measurements at appropriate temperatures. Indeed, even though the measured rate constants for the N( 2 D) + CH4 reaction are close to the previously recommended ones, the new uncertainties associated to this process might lead to an entirely different set of key reactions compared to earlier work. As further experimental studies of other important N( 2 D) reactions are currently underway, a more detailed uncertainty propagation study will be performed at a later date once their rate constants and associated uncertainties have been fully characterized.Here we present a temperature dependent kinetic study of the reactions of metastable nitrogen atoms in the 2 D state with the saturated hydrocarbons CH4, C2H6 and C3H8. A supersonic flow (Laval nozzle) reactor was used to perform these measurements at temperatures as low as 75 K. As N( 2 D) atoms are difficult to generate photolytically, a chemical source of these atoms was employed, through the reaction of C( 3 P) atoms with NO molecules as first described by Nuñez-Reyes & Hickson. 6 C( 3 P) atoms themselves were produced by the multiphoton dissociation of CBr4 molecules, while N( 2 D) atoms were detected directly through vacuum ultraviolet laser induced fluorescence. While the measured rate constants for the N( 2 D) + CH4 reaction are in good agreement with earlier work at higher temperature, the values derived for the N( 2 D) + C2H6 and N( 2 D) + C3H8 reactions are significantly smaller at low temperature than previous estimates based on earlier room temperature results. The effects of the new rate constants are tested on a coupled ion-neutral photochemical model of Titan's atmosphere. These simulations confirm that the N( 2 D) + CH4 reaction is the main loss process for N( 2 D) atoms in Titan's atmosphere in addition to being the major source of CH2NH. In contrast, the N( 2 D) + C2H6 and N( 2 D) + C3H8 reactions are shown to be negligibly important processes in the overall photochemistry of Titan's atmosphere. The need for uncertainty propagation studies to test the effects of the new rate constant uncertainties, alongside those of future measurements, is highlighted.

Table 1

 1 Supersonic flow characteristics

	Mach number	1.83± 0.02 a 1.99 ± 0.03	2.97 ± 0.06
	Carrier gas	N2	Ar		Ar
	Density (´10 16 cm -3 )	9.4 ± 0.2	12.6 ± 0.3	14.7 ± 0.6
	Impact pressure (Torr)	8.2 ± 0.1	10.5 ± 0.2	15.3 ± 0.5
	Stagnation pressure (Torr) 10.3	13.9		34.9
	Temperature (K)	177 ± 2	127 ± 2		75 ± 2
	Mean flow velocity (ms -1 ) 496 ± 4	419 ± 3		479 ± 3
	Chamber pressure (Torr)	1.4	1.5		1.2
	a The errors on the Mach number, density, temperature and mean flow velocity (1σ) are
	calculated from separate measurements of the impact pressure using a Pitot tube as a function
	of distance from the Laval nozzle and the stagnation pressure within the reservoir.
	Room temperature measurements were also performed during the present investigation. These
	experiments were conducted at reduced flow velocities without a Laval nozzle, to ensure that
	no pressure gradients existed within the reactor. In a similar manner to the study of Nuñez-
	Reyes & Hickson, 6 N( 2 D) atoms were generated chemically, through the reaction of ground
	state C( 3 P) atoms with nitric oxide
	C( 3 P) + NO ® N( 2 D, 4 S) + CO (1a)
		® O( 3 P) + CN	(1b)

  the pseudo-first-order losses of N( 2 D) with NO (𝑘 1( 2 3)?1@ [NO]), CxH2x+2 (𝑘 1( 2 3)?E F G 2FH2 [C I H KI?K ]), any other excess reagent (such as CBr4 for example (𝑘 1( 2 3)?ELM ' [CBr P ]) and through diffusion (𝑘 1( 2 3),RSTT ). As C( 3 P) atoms do not

	react with CxH2x+2, 𝑘 > < = 𝑘 E( U V)?1@ [NO] + 𝑘 E( U V),RSTT . Nevertheless, as the first 15
	microseconds of the N( 2 D) temporal profiles were not exploitable, it was more convenient to
	use a function of the form
	𝐼 1( 2 3) = Cexp(-𝑘 ; < 𝑡)

2 D) removal, 𝑘 > < is the pseudo-first-order rate constant for N( 2 D) production and t is time. The term 𝑘 ; < comprises several contributions including

Table 2

 2 Measured rate constants for the N( 2 D) + CxH2x+2 reactions

	Reaction	T / K	N a	[CxH2x+2]	[NO]	𝑘 1( 2 3)?E F G 2FH2 /
				/ 10 14 cm -3	/ 10 14 cm -3	10 -13 cm 3 s -1
	N( 2 D) + CH4	296	44	0 -46.4	4.3	(24.8 ± 5.4)
		296	24	0 -45.8	2.2	(21.0 ± 4.3)
		177	34	0 -55.5	2.6	(4.7 ± 2.4)
		127	20	0 -61.9	1.6	(1.7 ± 1.5)
	N( 2 D) + C2H6	296	28	0 -23.6	4.4	(122.0 ± 28.7)
		296	18	0 -24.4	6.3	(126.4 ± 25.7)
		177	39	0 -37.6	2.7	(25.2 ± 5.8)
		127	38	0 -43.2	3.1	(16.8 ± 3.7)
		127	18	0 -43.1	4.6	(17.6 ± 4.6)
		75	19	0 -21.3	2.7	(11.1 ± 3.5)
	N( 2 D) + C3H8	296	21	0 -24.9	6.3	(173.3 ± 37.9)
		177	38	0 -26.7	2.7	(56.1 ± 12.0)
		127	36	0 -30.9	3.2	(42.8 ± 9.0)
		127	18	0 -31.0	4.7	(45.5 ± 9.8)
		75	30	0 -10.5	2.7	(32.5 ± 8.8)
	a Number of individual kinetic profiles recorded for each experiment to derive the second-
	order rate constant.

Table 3

 3 Mole fractions at 1000 km produced by the two model runs. and despite the much larger rate constant). After the new rate constant is included, the N( 2 D) + C2H6 reaction becomes an entirely negligible source of CH2NH (0.3 %), with a decrease of the integrated rate by a factor of 8 between the two models.

	Species	Mole fraction	Difference (with respect to the
		(New Model)	standard model) / %
	C2H6	2.3 ´ 10 -5	+9.3
	C3H8	3.45 ´ 10 -7	+9.4
	C4N2	8.3 ´ 10 -7	+9.4
	NH2	8.7 ´ 10 -7	-10.8
	NH3	8.6 ´ 10 -7	-10.6
	CH2NH	1.75 ´ 10 -5	-11.1
	CH3NH2	9.35 ´ 10 -10	-10.6
	Due to the low mole fraction of C3H8 in Titan's atmosphere, the reaction of N( 2 D) with C3H8
	has a negligible influence on the chemistry in both models. For the other two reactions, the
	new rate constants lead to only small differences in the mole fraction profiles for all species at
	altitudes which correspond to the peak of the N( 2 D) mole fraction profile. Taking CH2NH as
	an example, this species is a major product of the N( 2 D) + CH4 reaction with the channel
	leading to CH2NH + H considered to represent 80 % of the total, while the products CH3 +
	CH2NH are considered to be the sole products of the N( 2 D) + C2H6 reaction in current
	models. By examining the integrated column rates, we note that the N( 2 D) + CH4 reaction
	represents the major source of CH2NH in the standard model (47 %), followed by the NH +
	CH3 ® CH2NH + H (30 %) and CH2NH2 + + e -reactions (although CH2NH2 + is almost
	entirely formed by proton transfer reactions with CH2NH, so this mostly just recycles between
	CH2NH and CH2NH2 + with only a minor loss channel to form NH2 + 3 CH2). This result also
	holds for the new model given the very minor change in the rate constant values for the N( 2 D)

+ CH4 reaction between the two simulations. In the standard model, the N( 2 D) + C2H6 reaction represents the third largest source of CH2NH in Titan's atmosphere, albeit with a contribution of only 2.5 % due to the relatively low mole fraction of C2H6 with respect to CH4 (

Acknowledgements K.M.H. and D.N.R acknowledge support from the French program ''Physique et Chimie du Milieu Interstellaire" (PCMI) of the CNRS/INSU with the INC/INP co-funded by the CEA and CNES as well as preliminary funding from the "Program National de Planétologie" (PNP) of the CNRS/INSU.

One sentence summary for TOC entry:

Low temperature rate constants for the N( 2 D) + C2H6, C3H8 reactions are shown to be much smaller than previously thought.