ISOTOPIC RATIOS IN TITAN's METHANE: MEASUREMENTS AND MODELING

The CIRS instrument consists of dual Fourier Transform spectrometers (mid-IR and far-IR) that share a common telescope, foreoptics, scan mechanism, and other devices (Kunde et al. 1996; Flasar et al. 2004). The mid-IR spectrometer is a conventional Michelson arrangement, with two parallel 1 × 10 arrays of photodetectors to sense the radiation. These arrays are known as focal plane 3 (FP3, 600–1100 cm⁻¹, 17–9 μm) and focal plane 4 (FP4, 1100–1400 cm⁻¹, 9–7 μm). Individual detectors in the arrays are square with a 0.273 mrad field of view (FOV), while the spectral resolution is variable (depending on commanded path difference of the movable mirror) from 0.5 to 15.5 cm⁻¹.

In limb-sounding mode the arrays are pointed above the rim of the body's solid surface (Nixon et al. 2009b), thereby sensing a much longer path length through the atmosphere that is especially suitable for increasing the signal-to-noise ratio (S/N) for detection of very low abundance trace species (Vinatier et al. 2007; Teanby et al. 2007; Nixon et al. 2009a) including isotopic variants (Bézard et al. 2007; Vinatier et al. 2007; Nixon et al. 2008a, 2008b; Jennings et al. 2008; Coustenis et al. 2008; Jolly et al. 2010). Usually, the instrument is pointed so that the arrays are perpendicular to the disk edge, thereby sampling 10 altitudes with each side-by-side array without repositioning. However, CIRS has only five receiver channels for the 10 detectors of each focal plane, so the receiver must alternate between odd-numbered and even-numbered detectors on each array to build up the full picture—sacrificing half the incident radiation.

Recently, several observations have been executed in the orthogonal orientation, where the arrays are instead positioned parallel to disk edge, with all 10 detectors of each array sampling only a single altitude. This provides additional S/N advantages because (1) the data from all 10 detectors may later be co-added, if horizontal variations in signal are negligible, and (2) the receivers may be commanded to work in "pair" mode, where the photon counts from detector pairs (1&2, 3&4 ... 9&10 on each array) are added together, enabling the entirety of the input signal to be used. In this mode, each pair of detectors form a "virtual" detector dimensions 0.273 × 0.566 mrad (in the vertical and horizontal directions, respectively).

Previously, several of these limb-parallel, detector pair mode observations were successfully analyzed to place upper limits on abundances of predicted, but undetected, trace species (Nixon et al. 2010). In the present study, we have analyzed data from the MIRLMPAIR002 observation described in Nixon et al. (2010), which was made during the T55 flyby on 2009 May 22, beginning at 02:26:41 UT with a 4 hr duration. The arrays were targeted at 25°S latitude (sampling ± 2° during the observation) and mean altitudes of 107 km (FP3, 7.6 mbar) and 247 km (FP4, 0.27 mbar). We use herein only the spectra recorded by FP4, 941 in all, at the highest spectral resolution of CIRS: 0.5 cm⁻¹ FWHM of instrumental line shape, after Hamming apodization. The range was ∼150,000 km at the mid-point of the observation, and the corresponding projected FOV size for a pair of mid-IR detectors was 41 × 85 km.

After co-adding all 941 FP4 spectra, two spectral regions were isolated for modeling. For the deuterated species, a small region from 1140 to 1160 cm⁻¹ was selected containing only the sharpest (Q-branch) emissions of ¹²CH₃D and ¹³CH₃D, to minimize problems with fitting the underlying baseline, due to haze absorption and the effects of a weak propane band (the ν₇ band of propane centered on 1158 cm⁻¹; see Figure 6(d) of Nixon et al. 2009a) for which no line data currently exist to enable modeling. For the non-deuterated species, a wider range was used, 1260–1350 cm⁻¹, to include the P-Q-R structure of both ¹²CH₄ and ¹³CH₄, giving the best possible S/N.

The radiative transfer calculation and fitting was accomplished using the NEMESIS software (Irwin et al. 2008), which predicts the emerging spectrum based on a one-dimensional vertical atmospheric model and iteratively adjusts this model to fit the spectrum using a minimization technique. The modeling was accomplished in a nearly identical way to that described in Nixon et al. (2009a). The model atmosphere included the major gases ¹²CH₄, N₂, and H₂ using the vertical abundances of Niemann et al. (2010); the minor species C₂H₂, C₂H₄, C₂H₆, and C₃H₈ with profiles as described in Nixon et al. (2009a); plus the methane isotopologues ¹³CH₄, ¹²CH₃D, and ¹³CH₃D with profiles initially set to scaled versions of the ¹²CH₄ at terrestrial isotopic abundances. Gas line data used were as described in Nixon et al. (2009a), except for C₂H₆ where we used the newly available list for the ν₆ and ν₈ bands by Di Lauro et al. (2012). Haze spectral opacity used was that of Vinatier et al. (2010), although we also experimented with gray haze and Khare et al. (1984) haze spectral properties and found negligible effect on our results.

The variable parameters were vertical temperature profile (information governed by ¹²CH₄ at fixed profile abundances), and scaling factors (one parameter for all levels) applied to the vertical haze profile, and the profiles of C₂H₂, C₂H₄, C₂H₆, C₃H₈, ¹³CH₄, ¹²CH₃D, and ¹³CH₃D. We tried retrieving in two methods (1) temperature first, then gas and haze abundances, or (2) temperature and gas/haze scale factors simultaneously. A slightly better fit resulted from method (2), which we retained for our final results. Both spectral regions were simultaneously fitted. Note that the finite footprint of the detectors, which convolves the limb radiation, was explicitly accounted for using eight weighted rays spaced 10 km apart, as described in Nixon et al. (2009b). Also, we determined that a small, sub-gridscale correction to the wavenumber scale provided in the CIRS database was required of +0.08 cm⁻¹, found by oversampling both model and data spectra to a grid of 0.01 cm⁻¹, and subsequent cross-correlation to obtain the optimum wavenumber shift. This correction is necessary due to slight drifts in the frequency of the reference laser that controls the interferogram sampling.

Figure 1 shows the spectral data (upper plots, black) and final model fit (upper plots, red). The lower plots show the spectral signature of the isotopologues, calculated by removing the minor CH₄ species from the model once the final fitting parameters were calculated, and then computing a synthetic spectrum using these retrieved parameters. The abundances of the three minor methane isotopologues (at 247 km) are found from the fitting (and using ¹²CH₄ = (1.48 ± 0.09) × 10⁻² above the cold trap from Niemann et al. 2010), while the error on the retrieved abundances is computed by NEMESIS according to the formalism given in Equation (22) of Irwin et al. (2008). Figure 2 shows how four isotopic ratio estimates can be derived from the abundances of the four methane isotopologues. The ratios and their errors are given in Table 1, along with the error-weighted mean ¹²C/¹³C and D/H.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 2 and Figure 3 compare these results to other recent estimates of D/H and ¹²C/¹³C, considering only measurements of the more abundant species (i.e., CH₄ and H₂). Our D/H from this work is indistinguishable from the value measured by Abbas et al. (2010), and somewhat greater though still compatible with the other measurements by CIRS and GCMS, and also the ground-based telescopes (error bars overlap). Our measured value of (1.59 ± 0.33) × 10⁻⁴ is fully consistent with the terrestrial value (1.56 × 10⁻⁴), while exceeding by one order of magnitude the value measured for Saturn ((1.6 ± 0.2) × 10⁻⁵; Fletcher et al. 2009). Also, our ¹²C/¹³C (86.5 ± 8.2) is in agreement with the recent re-measurements by Cassini INMS (Mandt et al. 2012) and Huygens GCMS (Niemann et al. 2010)—revised upward from earlier values due to improvements in the calibration and modeling—and also compatible the earlier CIRS measurements by Bézard et al. (2007) based on the pair (¹²CH₃D, ¹³CH₃D).

Zoom In Zoom Out Reset image size

With regard to the other previous CIRS measurement of ¹²C/¹³C, based on ¹²CH₄ and ¹³CH₄: ¹²C/¹³C = 76.6 ± 2.7 (Nixon et al. 2008a), our present finding is marginally compatible by virtue of overlapping error bars. However, the former conclusion that ¹²C/¹³C in methane is enriched with respect to the terrestrial/giant planet values, and those of heavier hydrocarbons on Titan (ethane and acetylene), is no longer secure. We believe that the change in the central value from ∼77 to ∼87 is due to different systematic errors in the nadir-sounding retrievals (Nixon et al. 2008a) and the limb-sounding retrievals (present work), with limb-sounding having the advantage that the weighting functions for each species are sharply peaked at the same (tangent) altitude. There may be some additional systematic error due to the non-overlapping line approximation made in the correlated-k radiative transfer calculation. Finally, the small error bar of the previous work was predicated on combination of purely random errors, allowing for infinite error reduction, without regard to error floors imposed by systematic errors (e.g., absolute transition strengths in line atlases). From these considerations, we recommend that the present result (86.5 ± 8.2) which has a more conservative error and is also in agreement with the recent GCMS and INMS re-analysis is considered to supersede the previous result.

Once the breakup of methane is initiated by photolysis, a series of reactions (e.g., Yung et al. 1984) leads to the formation of acetylene (C₂H₂), which is itself easily photolyzed to ethynyl (C₂H). Once this radical has been produced, it quickly attacks methane removing an H-atom (hydrogen abstraction) resulting in C₂H₂ and CH₃ (see reaction series S5 of Wilson & Atreya 2009). However, this reaction is catalytic in the sense that C₂H is recycled many times from C₂H₂, and because the photolysis of C₂H₂ can proceed at longer wavelengths than the photolysis of CH₄, this process effectively photosensitizes methane to longer-wavelength photolysis at altitudes below 700 km. This rapid reaction accounts for ∼46% of total methane destruction (photolysis plus chemical loss) and is the single largest source of methane loss (Wilson & Atreya 2009) overall.

Taking into account the possible isotopic varieties of methane, the possible reactions are (see Figure 4)

$$\begin{eqnarray*} {\rm C_{2} H} + {\rm ^{12}CH_{3}D} \rightarrow & {\rm C_{2} H_{2}} + {\rm ^{12}CH_{2} D} \quad [{\rm k1a}] \quad& {\rm R1A} \\ {\rm or} \rightarrow & {\rm C_{2} HD} + {\rm ^{12}CH_{3}} \quad [{\rm k1b}] \quad& {\rm R1B} \\ {\rm C_{2} H} + {\rm ^{12} CH_{4}} \rightarrow & {\rm C_{2} H_{2}} + {\rm ^{12}CH_{3}} \quad [{\rm k12}] \quad& {\rm R12} \\ {\rm C_{2} H} + {\rm ^{13} CH_{4}} \rightarrow & {\rm C_{2} H_{2}} + {\rm ^{13}CH_{3}} \quad [{\rm k13}] \quad& {\rm R13.} \end{eqnarray*}$$

Zoom In Zoom Out Reset image size

We may define the Kinetic Isotope Effect (KIE) ratios for this reaction as follows: q₁₃ = k₁₂/k₁₃ and q_D = k₁₂/(k_1a + k_1b). These fractions are >1 if destruction of ¹²CH₄ proceeds faster than the alternative (the so-called normal KIE) and <1 if the rare isotopologue is depleted faster. However, due to the difference in ground state zero-point vibrational energies (ZPEs) for the different methane isotopologues, the reaction rates k should differ by some degree and so q ≠ 1.0. The ZPE shifts for CH₃D and ¹³CH₄ have been estimated by Nair et al. (2005) to be −599 ± 1 and −25.5 ± 0.5 cm⁻¹, respectively, relative to ¹²CH₄ (9834 cm⁻¹), therefore the reaction rate for ¹³CH₄ should more closely resemble that of ¹²CH₄ than will CH₃D, and the KIE q should be closer to unity.

The KIEs for methane in a variety of hydrogen abstraction reactions have been extensively measured for terrestrial atmospheric applications, focusing on CH₄ + Cl and CH₄ + OH, but CH₄ + CN has also been investigated (see Table 3). These show that hydrogen KIE is typically tens of percent (q_D ∼ 1.1–2.0) while carbon KIE is a few percent (q₁₃ ∼ 1.01–1.08). The KIE for CH₄ + C₂H has been studied only once at cold temperatures, and only for the larger hydrogen effect (Opansky & Leone 1996)—although a recent paper by Matsugi et al. (2011) added data points at warm temperatures (>290 K). The Opansky & Leone (1996) work, undertaken for application to Titan, measured the reaction rates of both CH₄ and CD₄ with C₂H down to 190 K. We have followed the proposal of Lunine et al. (1999) that an effective reaction rate for CH₃D can be assigned by adding

$$\begin{equation} k_{{\rm CH}_3{\rm D}} = (0.25 \times k_{{\rm CD}_4}) + (0.75 \times k_{{\rm CH}_4}). \end{equation} \tag{ 1 }$$

The KIE for hydrogen in this reaction is then $q_D = k_{{\rm CH}_4} /k_{{\rm CH}_3{\rm D}}$ and is given in Table 3.

However, the carbon KIE for the reaction CH₄ + C₂H is also of great interest because the ratio ¹²C/¹³C in methane may potentially be more informative about the history of Titan's methane than D/H. This is because we have a larger degree of certainty in the starting value of ¹²C/¹³C—which varies very little throughout the solar system—compared to D/H, which varies widely. In the absence of laboratory measurements of q₁₃, we have therefore attempted to constrain the value through ab initio calculations, as described in the following subsection.

The success of ab initio calculations of chemical kinetics in hydrogen abstraction reactions, including kinetic isotope effects, has previously been demonstrated (Sellevåg et al. 2006; Temelso et al. 2006; Vandeputte et al. 2007). Here we apply high-level ab initio methods and conventional transition state theory (TST; Eyring 1935) coupled with several one-dimensional tunneling corrections of successively increasing complexity to estimate the reaction rates (R12, R13). The reaction barrier and thermodynamics were determined using the UCCSD(T) (unrestricted coupled-cluster singles, doubles, and perturbative triples) level of theory employing correlation consistent basis sets extrapolated to the complete basis set (CBS) limit. ZPE and thermal corrections were included using UCCSD(T) with the cc-pVTZ (correlation-consistent polarized valence triple-ζ) basis set while a first-order correction to the adiabatic approximation was included using the diagonal Born-Oppenheimer Correction (DBOC) at the UCCSD/cc-pVTZ level of theory.

The rate enhancement due to tunneling was determined by calculating the transmission coefficient through three different one-dimensional reaction potentials. The Wigner correction (Wigner 1932) models the reaction barrier as an untruncated parabola of a certain width. The Skodje–Truhlar approach (Skodje & Truhlar 1981) considers the reaction barrier as a truncated parabola of a defined width and height while the Eckart method (Eckart 1930) fits the zero-point-corrected energy of the reactants, transition state, and products to an asymmetrical Eckart function. A detailed description is given in the Appendix.

Our first objective was to attempt to replicate the hydrogen KIE for the reaction C₂H + CD₄ → C₂HD + CD₃ as measured by Opansky & Leone (1996), and which can be used to infer R1A + R1B. The results of our calculations are shown in Figure 5 using the various tunneling corrections, and generally show good agreement with the laboratory data. This validation was extremely important, as it showed that the model was working well and therefore we should feel confident to extrapolate our scheme to the unmeasured abstraction of hydrogen from ¹³CH₄. Using the method of (Lunine et al. 1999; Equation (1)) we find q_D = 1.148 ± 0.001 at 298 K and q_D = 1.212 ± 0.002 at 175 K, from the mean and standard deviation of our TST calculations.

Zoom In Zoom Out Reset image size

We therefore proceeded to use TST to compute the rate for R13 and also the q₁₃, which is plotted in Figure 6, and values listed in Table 4. As expected, the KIE for carbon is much less than that for hydrogen, close to 2% across the range of calculated temperatures. Using the standard deviation between the curves as the error, we find q₁₃ ≃ 1.018 ± 0.001 at 298 K and ≃ 1.019 ± 0.002 at 175 K.

Zoom In Zoom Out Reset image size

In the next section, we will show the application of these q_D and q₁₃ estimates to modeling the evolution of Titan's atmosphere.

In the case of no methane resupply, we will expand on the model of Cordier et al. (2008), considering possible fractionation that could occur during the five methane depletion/creation mechanisms identified by Wilson & Atreya (2009). These five processes are included in a generalization of the two-process model of Cordier et al. (2008), which we will use for both D/H evolution, but also ¹²C/¹³C evolution. We replace their Equation (3) (the Rayleigh distillation relation) with

$$\begin{eqnarray} \log R(t) & = & \left(1 - \frac{\sum _{i} \left(F_i / q_i \right)}{\sum _{i} F_i} \right)^{-1} \log E(t) \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} \ & = & \left(1 - \frac{1}{q_{\rm eff}} \right)^{-1} \log E(t), \end{eqnarray} \tag{ 3 }$$

where E(t) is the enrichment of the heavy isotope in methane (E₁₃(t) = f₁₃(t)/f₁₃(0) and E_D(t) = f_D(t)/f_D(0)); f₁₃(t) and f_D(t) are the (heavy/light) methane mixing ratios (¹³CH₄/¹²CH₄ and ¹²CH₃D/¹²CH₄) at time t; R(t) is the reservoir ratio, the ratio of initial methane (¹²CH₄ ≃ total methane) in the atmosphere to the amount of methane today; q_i is the fractionation due to the ith process, defined such that q > 1 in the normal case where the lighter isotope is lost faster (inverse definition to Cordier et al. 2008); and F_i is the rate of methane destruction due to the ith process. Also note that we have defined 1/q_eff as the single value that replaces the mean of (1/q_i) weighted by the importance of each process (F_i): 1/q_eff = ∑_i(F_i/q_i)/∑_iF_i which simplifies both the formula and the concept. For convenience, we will hereafter cease to explicitly write out the time dependence of E and R in each case.

Within this general model, we allow six sub-models for different assumptions about the fractionation processes: (1A) no escape, chemical loss rates from Lavvas et al. (2008; hereafter L08); (1B) as in 1A, but using chemical loss rates from Wilson & Atreya (2009; hereafter WA09); (1C) as in 1A, but using instead chemical loss rates from Krasnopolsky (2009; hereafter K09). (1D), (1E), and (1F) are identical to 1A, 1B, and 1C, respectively, but with additional fractionation by hydrodynamic escape (Yelle et al. 2008; Strobel et al. 2010). In this paper, we do not compute models with sputtering loss of methane, which depends on methane density and is therefore not so easily included in the Rayleigh formula. In any case, the present day sputtering loss rate of 2.8 × 10⁷ cm⁻² s⁻¹ (De La Haye et al. 2007) is much smaller than the rates due to photochemistry and possible hydrodynamic loss 10⁹–10¹⁰ cm⁻² s⁻¹, therefore adds only a small perturbation to the results obtained with these assumptions. Fractionation due to sputtering is discussed further in a companion paper (Mandt et al. 2012).

Table 5 gives the values for the F_i and q_i used in each model. The F_i for each photochemical model were in general obtained from the respective papers, although for L08 the requisite numbers were not given in the paper, but kindly supplied by P. Lavvas upon request. We note that there is significant variation in the values for F_i between the three models (see Figure 7), especially in the proportion of chemical depletion attributed to the reaction of CH₄ with C₂H. In particular, only 26% of the chemical loss in K09 is attributed to the reaction of ethynyl with methane—the particular chemical reaction for which we have a measured KIE—compared to 60% in WA09 and 73% in L08, which therefore produces a very substantial change in the inferred fractionation times, as discussed later. In addition, all three photochemical models (K08, WA09, and K09) use 0.3 as the quantum yield for C₂H formation following acetylene photolysis, while recent laboratory measurements (Läuter et al. 2002; Kovács et al. 2010) indicate a yield closer to unity. Thus, all these models may significantly underestimate the C₂H production rate on Titan and similarly the importance of the subsequent C₂H + CH₄ loss process.

Zoom In Zoom Out Reset image size

We do not consider the time evolution of the solar UV flux, which is expected to decrease over time as discussed in Mandt et al. (2012). This introduces an error (about 6% decrease in flux over 200 Myr) that affects the chemical loss rate by a similar amount. This error becomes significant when the KIE is very small (q_eff ≃ 1) and/or E becomes significantly larger than unity, and is discussed further in Section 5.

All six models use the same q_i values. The q from photolysis follows estimates given by Nair et al. (2005), while q due to the KIE of hydrogen abstraction in methane by ethynyl (or "photosensitized" destruction, in the parlance of Wilson & Atreya 2009) was estimated earlier in this paper for both hydrogen and carbon KIE. The q for other chemical processes is unknown, and set to unity, as is that for chemical production of methane. This model deficiency due to unknown KIE magnitudes causes q_eff to be underestimated, and so the enrichment timescales are overestimated by some undetermined amount and should be considered upper limits in this sense.

In models 1D, 1E, and 1F we add an additional term for methane loss and fractionation by hydrodynamic escape (F_i): 1.35, 2.90, and 3.48 (× 10⁹ cm⁻² s⁻¹) in L08, WA09, and K09, respectively, and is a significant fraction of the total methane loss. This loss is also assumed constant over time. For q we use the ratio (¹²C/¹³C)_exobase/(¹²C/¹³C)_surface derived from a new analysis of the INMS data set by Mandt et al. (2012) and is set identically for both heavy isotopologues CH₃D and ¹³CH₄. Table 5 also gives both the total F from all processes, calculated separately for the six sub-models, and the q_eff in each case.

To actually proceed with obtaining predictions from the model, we need to impose some constraints. Our starting point will be to presume that we have knowledge of the enrichment in ¹³C at the present day, i.e., E₁₃, from which we then have sufficient information ((q₁₃)_eff) to calculate R using Equation (3). We will justify this assumption, by arguing that ¹²C/¹³C has proved to be surprisingly constant throughout the solar system. Outer solar system planetary values have been derived for the atmospheres of Jupiter (92.6^+4.6_{− 4.1}; Niemann et al. 1998) and Saturn (91.8^+8.4_{− 7.8}; Fletcher et al. 2009), in good agreement with a mean ¹²C/¹³C derived for eight comets (90 ± 4; Hutsemékers et al. 2005) and even with the terrestrial inorganic standard value (Pee Dee Belemnite, 89.4; Werner & Brand 2001).

In this work, we use the error-weighted mean of the Jovian, Saturnian, and cometary ratios: 91.3^+2.8_{− 2.7} as the range of the initial ¹²C/¹³C. Similarly, we derived a mean ratio for present-day Titan from the observations: 89.7 ± 1.0 (Table 2). Combining these gives us a range for the present-day enrichment in E₁₃ of 1.017^+0.034_{− 0.032}, or effectively 1.000 < E₁₃ < 1.051 as we do not have any processes in our model that increase the carbon isotopic ratio (i.e., cause E to decrease). The alternative option, beginning with a presumed D/H value, is much riskier, since D/H varies by over an order of magnitude from the giant planets to terrestrial planets and comets. Therefore, we will compute the E_D afterward based on timescales derived from E₁₃.

The assumed range of E₁₃ values are shown in Column 2 of Table 6. The corresponding R, given in Column 3, gives us the amount of methane implied in the original methane atmosphere in units of the current methane amount (see also Figure 8). We can then apply Equation (5) of Cordier et al. (2008), with slight modification, to compute the time for this change (depletion) to occur:

$$\begin{equation} t = \frac{M_{\rm CH_4} (R-1 ) }{m_{\rm CH_4} \sum _i F_i }. \end{equation} \tag{ 4 }$$

In this formula $M_{{\rm CH_4}}$ is the present column mass of CH₄ =189.1 g cm⁻², calculated by us by direct quadrature of our atmosphere model; $m_{{\rm CH_4}}$ is the mass of a CH₄ molecule =(16.04/N_A) g, where N_A = 6.022 × 10²³; and F_i has been previously defined. (R − 1) is used rather than R, so that the present inventory is not included in the calculation. The results are given in Column 4 of Table 6.

Zoom In Zoom Out Reset image size

We then use R to compute E_D using Equation (3) and replacing the (q₁₃)_eff with (q_D)_eff. The estimate of the deuterium enrichment can be used to compute the original D/H in Titan's atmosphere, if a present day-value is assumed, e.g., from Table 2. We note that the maximum computed enrichment in D/H is ∼1.65, much less than required to enrich a Saturnian D/H ((1.6 ± 0.2) × 10⁻⁵; Fletcher et al. 2009) to any realistic Titan value, which is typically one order of magnitude greater (see Table 2).

Finally, we have computed the mass of carbon that has been processed during the methane depletion from M = M_atm(R − 1) × (12.00/16.04), where M_atm is the mass of methane in the present atmosphere, calculated by us to be 1.576 × 10⁵ GT (gigatons). The implications of the D/H and carbon mass ranges are discussed further in Section 5.

In the second scenario, we assume that methane is being replenished to the atmosphere at a rate equivalent to its net destruction by all mechanisms. In the absence of resupply all Titan's methane is predicted to vanish in the next ∼20 Myr, therefore this model avoids the problem posed by our observing a rather unique and unsustainable point in Titan's history.

It is possible to show that in this case, a simple relation holds between the isotopic ratio in the atmosphere and the isotopic ratio in the supply flux as follows. We consider two isotopes of methane: "light" (L, i.e., ¹²CH₄) and "heavy" (H, i.e., either ¹³CH₄ or ¹²CH₃D), with atmospheric molecular populations N_L(t) and N_H(t), and total loss rates from all mechanisms l_L(t) and l_H(t). The general time-evolution equations for the populations of each molecular species are

$$\begin{eqnarray} \frac{dN_L}{dt} & = & -l_L(t) N_L(t) + \frac{f_0}{1+f_0}\phi (t) \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \frac{dN_H}{dt} & = & -l_H(t) N_H(t) + \frac{1}{1+f_0}\phi (t), \end{eqnarray} \tag{ 6 }$$

where ϕ(t) is the total methane supply flux (e.g., from an internal, or perhaps external source) with constant isotopic ratio f₀. Therefore when supply balances loss, the relation ϕ(t) = l_L(t)N_L(t) + l_H(t)N_H(t) must hold. If in addition, the loss processes are unchanging with time (or changing on long timescales compared to equilibration of molecular species) and the supply constant then nothing perturbs the populations of both isotopologues and dN_L/dt = dN_H/dt = 0. Therefore,

$$\begin{eqnarray} l_L N_L & = & \frac{f_0}{1+f_0} \phi (t) \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} l_H N_H & = & \frac{1}{1+f_0} \phi (t). \end{eqnarray} \tag{ 8 }$$

Dividing these last equations yields

$$\begin{equation} \frac{N_L}{N_H} = f_0 \frac{l_H}{l_L} \Rightarrow f = f_0 / q_{\rm eff}, \end{equation} \tag{ 9 }$$

as q_eff is just the total effective fractionating factor we met earlier. In other words, the (light/heavy) isotopic ratio ((D/H)⁻¹ or ¹²C/¹³C) stabilizes under steady resupply conditions at a value equal to the isotopic ratio in the source flux, divided by the effective (weighted) fractionation factor, also assumed constant. If, however, the loss conditions change over time (e.g., changing solar flux), then the supply must also change to maintain balance.

Assuming the ¹²C/¹³C of supply is the outer solar system value we estimated earlier (f₀ = 91.3) we can use the six (q₁₃)_eff derived for Model 1 (see Table 5) and Equation (9) to derive six steady state predictions for ¹²C/¹³C today: 90.2 (2A: L08, no escape); 90.4 (2B: WA09, no escape); 90.9 (2C: K09, no escape); 88.0 (2D: L08, no escape); 85.1 (2E: WA09, escape), and 84.5 (2F: K09, escape). The values for the first three (non-escape) models are very compatible with our estimated Titan mean value (89.7 ± 1.0) while the values for the escape models are compatible with individual measurements in Table 2 other than GCMS. Similarly, (CH₃D/CH₄)_atm = (q_D)_eff (CH₃D/CH₄)_supply, so if the present day value is close to terrestrial, then the source must be some ∼4%–15% lower in D/H (depending on model) to allow for subsequent enrichment. In this model scenario, it is not possible to put any constraint on the amount of time for which methane has been resident in the atmosphere because once a steady state has been reached there is no evolution over time.

Note that this steady state condition has previously been advocated by Jennings et al. (2009), who measured a zero-enrichment (terrestrial) value for ¹²C/¹³C in ethane and contrasted it with the enrichment seen in ¹³CH₄. They suggested that a carbon KIE with q₁₃ ≃ 1.1 could produce a stable enhancement in methane's ¹²C/¹³C consistent with the (preliminary) GCMS value of 82.3 ± 1.0 (Niemann et al. 2005). Since this work, both the enhancement in the present-day Titan ¹³CH₄ and the q₁₃ have been substantially revised closer to unity, allowing the scenario to remain self-consistent.

The third possibility is that atmospheric methane is in a non-steady state, but increasing, rather than decreasing (as in Model 1). In this case, the supply of methane occurs faster than depletion, and therefore any effect of isotopic enrichment is washed out by the greater supply of fresh methane. In this model, the isotopic ratios are maintained arbitrarily close to the supply values, and must exceed (for ¹³C, or be lesser than, for D) the steady state isotopic ratios: (¹²CH₄/¹³CH₄)_supply/(q₁₃)_eff < (¹²CH₄/¹³CH₄)_atm < (¹²CH₄/¹³CH₄)_supply and (CH₃D/CH₄)_supply < (CH₃D/CH₄)_atm < (q_D)_eff (CH₃D/CH₄)_supply.

This model would be strongly favored if, for example, Titan turned out to have no present endogenic activity, as argued by Moore & Pappalardo (2011). Contrariwise, this creates the Copernican problem mentioned previously that we are observing Titan during a unique period in its history, when the final remnants of a (possibly much larger) methane inventory are finally vanishing. Using the probabilistic version of the Copernican argument advocated by Gott (1993) for continuous events, there is a 50% probability that Titan had no more than twice the present methane present in the atmosphere, and 90% probability that the initial inventory was no more than 10 × the present amount, prior to the beginning of the current depletion cycle. This initial input may have come from a phase of episodic outgassing as envisioned by Tobie et al. (2006), with the most recent event occurring 400–1400 Myr ago.

The results from our Model 1 indicate a very large range of possible enrichment times based on the measured present and assumed initial values for ¹²C/¹³C. Using the mean values for enrichment, we find timescales of 4–9 Myr if hydrodynamic escape is occurring, increasing to 56, 112, or 1646 Myr based on the no-escape photochemistry of L08, WA09, and K09, respectively. Using maximum possible enrichments increases these times dramatically: up to 47 Myr for hydrodynamic escape, and 1.1, 4.3, and 6600 Gyr for the chemistry of L08, WA09, and K09—the last of which is clearly implausible. Note that while the models of L08 and WA09 agree more closely on the chemical depletion rates compared to K09, the models of WA09 and K09 agree more closely on the escape rate than does L08. This reinforces the earlier observation that our model is extremely sensitive to input parameters, and better agreement between the photochemical models for F_i due to C₂H + CH₄ in particular would be helpful.

A drawback exists in our simple model, namely, the lack of time dependence in the solar flux, and hence destruction rate. Calculations using a time-dependent model of methane depletion (Mandt et al. 2012) indicate much shorter timescales: 208 and 586 Myr, respectively, for E = 1.017 and E = 1.051 using the chemistry of K09 and no escape (cf. 1.6 and 6600 Gyr above). Therefore, we conclude that using the simple analytical Rayleigh distillation relationship in the presence of time-dependent processes can result in large errors, and should only be used for order-of-magnitude estimates when the enrichment times are short compared to the timescales of time-dependent processes.

Although many uncertainties exist, it does seem possible to conclude that models that include escape can achieve the inferred present levels of carbon fractionation in short times, no more than one present methane lifetime (∼20 Myr). On the other hand, fractionation by chemistry alone will take much longer, hundreds to thousands of Myr, based on our model results for most likely E₁₃ (56–1646 Myr) and the 586 Myr upper limit derived using the time-dependent model. We note that Hörst et al. (2008) have estimated that ∼300 Myr is required to create the present inventory of atmospheric CO, which requires both the presence of carbon (from methane) and oxygen (O⁺, deposited at the top of the atmosphere). If this model is interpreted as a minimum time for the existence of methane in the atmosphere, then models that allow fractionation during hydrodynamic escape (1D–1F) must be excluded.

Additional considerations may include the age of the surface, if the crust underwent significant upheaval and resurfacing during a major outgassing event and has subsequently aged only due to cratering and mechanical erosion. Neish & Lorenz (2012) recently estimated the surface age in two slightly different impact models as 200 Myr or 1000 Myr: incompatible once again with our hydrodynamic escape models, which fractionate too rapidly. We therefore arrive at an estimated range for the methane duration of 300–600 Myr, bounded by the CO chemistry and the time-dependent model for methane depletion, and contained inside the surface age range. We note that the probabilistic Copernican argument (Gott 1993), however, implies that there are only 6% and 3% likelihoods respectively of observing the last ∼20 Myr of a 300 or 600 Myr continuous process. This is a quantitative way of saying that it appears unlikely that we are observing a narrow period in Titan's history when the last methane of a much larger initial inventory is disappearing.

Regarding the implied amount of carbon processed during the methane depletion, we find most-likely amounts (in GT) of 370,000 (1A), 580,000 (1B); 8,200,000 (1C); 70,000 (1D); 34,000 (1E); 30,000 (1F). We can infer that this total carbon is ultimately deposited on the surface, although we must correct the latter three totals for carbon escape, arriving at: 63,000 (1D); 26,000 (1E); 22,000 (1F). These may be compared with the surface deposits of carbon estimated by Lorenz et al. (2008): 16,000–160,000 GT (lakes) and 160,000–640,000 GT (dunes). Therefore, considering these six cases where methane is not resupplied to the initial atmosphere, cases 1D–1F with hydrodynamic escape do not produce sufficient methane to explain the minimum surface inventory (176,000 GT). Case 1C also probably produces too much, while 1A and 1B fall inside the nominal range.

As we have previously stated, Model 1 presents a significant problem in that it requires us to be observing a unique period of Titan's history (the Copernican Problem). Therefore, much attention has been focused on mechanisms for replenishing the methane, in particular a search for endogenic eruptive features ("cryovolcanic flows") that could supply methane, trapped originally in Titan's interior as clathrate hydrates, to the atmosphere (e.g., Fortes et al. 2007; Choukroun et al. 2010). Some early putative volcanic features from Cassini mapping (Sotin et al. 2005; Lopes et al. 2007) have later been disproven, while more recent evidence (e.g., Nelson et al. 2009; Wall et al. 2009) has been disputed (Soderblom et al. 2009; Moore & Pappalardo 2011). A full discussion of the arguments for and against endogenic activity on Titan is beyond the scope of this paper, and changes rapidly due to discoveries made by the ongoing Cassini mission. Suffice it to say that while such phenomena present significant advantages for explaining the present atmosphere, the evidence as yet is insufficient to sustain or disprove the hypothesis.

The case of the total surface carbon deposits (176,000–800,000 GT as discussed in previous section) is potentially more serious, as the steady state scenario processes about 120,000 GT of carbon (the current atmospheric inventory) every 20 Myr, or 6000 GT Myr⁻¹. Therefore, durations of this state beyond ∼130 Myr produce in excess of 800,000 GT and require explanations for the missing surface carbon, possibly reburied or otherwise incorporated into the crust, as advocated by Wilson & Atreya (2009). In summary, while Model 2 is compatible with available isotopic evidence from Titan's methane, resolving this scenario with the other existing evidence, especially the surface hydrocarbon deposits, remains to be proven.

It is unlikely that Titan's methane is in an exact steady state with supply equaling depletion (Model 2). If methane resupply is ongoing, it is likely sporadic, or else either decreasing (variant of Model 1) or increasing (Model 3) at the present time. Steadily increasing atmospheric methane could occur due to the gradual brightening of the Sun, slowly thawing out Titan's crust. In this case, the surface would also likely be active and young, and similar evidence for endogenic activity might be required as in Model 2. However, this scenario could help with the issue of carbon buildup at the surface, if methane was less abundant in the past. At the present time, isotopic information from methane can neither confirm not deny such a hypothesis.

Modeling hydrogen abstraction barriers involving radicals is challenging because of spin contamination in many ab initio wavefunction methods (Andrews et al. 1991), and large self-interaction errors in density functional methods (Zhang & Yang 1998), as demonstrated for the case of ethynyl (Temelso et al. 2006) and other radical species (Chuang et al. 2000; Vandeputte et al. 2007). A proven way to overcome spin contamination is to use the robust CCSD(T) (Raghavachari et al. 1989) method with a large basis set. Here, we fully optimized the geometry of the stationary points (reactants, transition state, and products) and determined their harmonic vibrational frequencies using the unrestricted CCSD(T) with cc-pVTZ basis set. We utilized the analytic gradient and Hessian code implemented in the CFOUR package (Stanton et al. 2010).

Employing the cc-pVXZ (X = T, Q, 5, ... , ∞) (Dunning 1989; Kendall et al. 1992) basis sets has the great advantage that the UCCSD(T)/cc-pVXZ electronic energies systematically converge to their infinite (complete) basis set (CBS) limit with increasing cardinal number X. The UCCSD(T)/cc-pVXZ energy is decomposed into its reference Hartree–Fock energy (E^HF_X) which converges exponentially to the HF CBS limit (E^HF_CBS) and the correlation energy (E^C_X) which converges as X⁻³ to the E^C_CBS. We extrapolated E^HF_X (X = T, Q, 5) to the E^HF_CBS by fitting to Equation (A1) (Feller 1993) and the E^C_X (X = T, Q) using Equation (A2) (Helgaker et al. 1997):

$$\begin{equation} E{^{{\rm HF}}_{{\rm CBS}}}=E{^{{\rm HF}}_{X}} + {\rm Be}^{-CX} \end{equation} \tag{ A1 }$$

$$\begin{equation} E{^{C}_{{\rm CBS}}}=\frac{X^3E{^{C}_{X}}-(X+1)^3E{^{C}_{X+1}}}{X^3-(X+1)^3} \end{equation} \tag{ A2 }$$

$$\begin{equation} E{^{}_{{\rm CBS}}}=E{^{{\rm HF}}_{{\rm CBS}}}+E{^{C}_{{\rm CBS}}}. \end{equation} \tag{ A3 }$$

The extrapolated E_CBS (A3) was used to calculate benchmark quality reaction barrier and energy within the framework of the Born-Oppenheimer (BO) approximation. Table 7 shows the convergence of the UCCSD(T) reaction barrier and energy toward the CBS. A first-order correction to the BO energy was included using the DBOC using UCCSD/cc-pVTZ at the UCCSD(T)/cc-pVTZ optimized geometry. DBOC is a first-order adiabatic correction to the BO approximation, and instead of assuming that nuclei are infinitely heavy, it takes into account the finite mass of the nuclei (Handy et al. 1986). While it is a small correction (Valeev & Sherrill 2003), its mass dependence is important in calculating reaction rates and KIEs. See Table 8.

Table 8. Computed Thermochemical Reaction Profile for R2 and R3

	12CH4+C2H	13CH4+C2H	CD4+C2H
Imaginary mode
$||\tilde{\nu }^\ddagger ||$(cm−1)	238.954	235.358	216.919
Reaction barriera,b,c
ΔE‡b	1.288	1.288	1.288
ΔE‡[+DBOC]	1.319	1.320	1.293
ΔE‡(0 K)	0.564	0.563	0.924
ΔE‡(0 K)[+DBOC]	0.596	0.595	0.930
ΔG‡(175 K)[+DBOC]	2.130	2.136	2.612
ΔG‡(298.15 K)[+DBOC]	3.184	3.194	3.694
Reaction energya,b,c
ΔEb	−27.602	−27.602	−27.602
ΔE[+DBOC]	−27.622	−27.623	−27.629
ΔE(0 K)	−29.552	−29.528	−28.732
ΔE(0 K)[+DBOC]	−29.572	−29.549	−28.760
ΔG(175 K)[+DBOC]	−29.600	−29.578	−28.801
ΔG(298.15 K)[+DBOC]	−29.624	−29.606	−28.894

Notes. ^aIn kcal mol⁻¹. ^bUCCSD(T)/CBS electronic energy. ^c[+DBOC] included at the UCCSD/cc-pVTZ level; (0 K) implies the inclusion of ZPVE corrections using UCCSD(T)/cc-pVTZ vibrational frequencies; (175 K) and (298.15) include thermal corrections on top of the DBOC and ZPVE corrections.

Download table as:  ASCIITypeset image

The above scheme provides very accurate electronic energies for the reactants, transition state, and products and one can easily calculate the classical barrier (ΔE‡) to the hydrogen abstraction from methane and its isotopologues by ethynyl radical. To get further thermodynamic information, we need to calculate the partition function for the reactants and transition state and employ TST (Eyring 1935) to get reaction rates. First, the zero-point vibrational energy correction (ΔZPVE), is estimated simply as one-half of the sum of the harmonic vibrational frequencies and added to the classical barrier to give ΔE‡(0). We obtain finite temperature corrections to the enthalpy ΔH(T) ΔS(T) using the vibrational, rotational, and translational partition functions in conjunction with the harmonic oscillator, rigid rotator, and particle-in-a-box models (Mcquarrie 2000; Cramer 2004). Adding these thermodynamic corrections to the classical barrier gives the enthalpy (ΔH‡) and Gibbs free energy (ΔG‡) of activation. Given the these quantities, one can use TST to get the reaction rate [k(T)]:

$$\begin{equation} k(T)=\kappa (T)\left(\frac{k_bT}{h}e{^\frac{-\Delta H^\ddagger + T\Delta S^\ddagger }{RT}}\right)=\kappa (T)\left(\frac{k_bT}{h}e{^\frac{-\Delta G^\ddagger }{RT}}\right). \end{equation} \tag{ A4 }$$

Here κ(T) is the transmission coefficient, k_B is the Boltzmann constant, h is Planck's constant, and R is the universal gas constant. κ(T) is unity in the classical (no tunneling) case; however, in reactions involving light atoms like hydrogens, quantum mechanical tunneling through the reaction barrier can increase the transmission coefficient, and hence enhance the reaction rates. The effect of tunneling is typically most pronounced at low temperatures where it can be the dominant reaction mechanism.

Ignoring all tunneling pathways but the one dimension along the reaction coordinate, we have estimated the tunneling transmission coefficient, κ(T). The simplest way to account for the effect of quantum mechanical tunneling on reaction rates is to use the Wigner approximation (Wigner 1932). It assumes that the reaction proceeds along a parabolic potential whose width is defined by the imaginary frequency ($\tilde{\nu } ^\ddagger$) of the transition state. The Wigner transmission coefficient κ_w(T) is

$$\begin{equation} \kappa _w(T)=1+ \frac{1}{24}\left(\frac{hc\tilde{\nu }^\ddagger }{k_bT}\right)^2, \end{equation} \tag{ A5 }$$

where $\tilde{\nu }^\ddagger$ is the magnitude of the imaginary frequency along the reaction coordinate in wavenumbers and c is the speed of light. This simplification is valid only when k_BT ≫ hc$\tilde{\nu }^\ddagger$. Using our UCCSD(T)/cc-pVTZ $\tilde{\nu }^\ddagger$ of 239 and 235 cm⁻¹ for the ¹²CH₄–C₂H and ¹³CH₄–C₂H transition states, we can see that the Wigner tunneling correction will work well only for T ≫ 340 K.

The Skodje–Truhlar (Skodje & Truhlar 1981) tunneling correction is a more robust approximation that is valid over a large temperature range because it factors in not only the width of the barrier but also its height. The Skodje–Truhlar transmission coefficient for an exothermic reaction with α ⩾ β like ours is given by

$$\begin{equation} \kappa _{{\rm ST}}(T)=\frac{\beta \pi /\alpha }{\sin (\beta \pi /\alpha)}-\frac{\beta }{\alpha -\beta }e^{\left[(\beta -\alpha)\Delta E^\ddagger (0)\right]}, \end{equation} \tag{ A6 }$$

where α = $\frac{2\pi }{h\tilde{\nu }^{\ddagger }}$, β = $\frac{1}{k_bT}$, and ΔE‡(0) is the zero-point energy corrected barrier.

Unlike the Wigner and Skodje–Truhlar corrections which assume a parabolic form for the reaction potential, the Eckart tunneling method (Eckart 1930) fits the zero-point-corrected energies of the reactants, transition state, and products to an asymmetric Eckart function. It models the reaction potential more realistically and it has analytical solutions. For a full description of the Eckart tunneling correction, the reader is encouraged to see Johnston & Heicklen (1962) or Vandeputte et al. (2007). A summary of the q₁₃ KIE using the different tunneling approximations is shown in Table 4.

The three zero-curvature tunneling methods described above assume tunneling occurs in one dimension along the reaction path. There exist other more complicated ways of calculating the reaction rate (such as variational transition state theory) and accounting for multi-dimensional tunneling. Truhlar et al. (1996) provide a good summary of these approaches.