© ESO, 2015

1. Background and motivation

In recent years, the presence of methane in the Martian atmosphere has been proposed by several authors (Krasnopolsky et al. 2004; Formisano et al. 2004; Mumma et al. 2009; Fonti & Marzo 2010; Krasnopolsky 2012; Villanueva et al. 2013), but some of the previous claims have been questioned and doubts have been expressed for each of them (Lefèvre & Forget 2009; Zahnle et al. 2011; Kerr 2012). In this context, the latest results of the Tunable Laser Spectrometer (TLS) of the Sample Analysis at Mars (SAM; Mahaffy et al. 2012; Webster et al. 2013a,b, 2015) are of great relevance. From in situ measurements at Gale crater, an upper limit of 1.3 ppbv for the Mars atmospheric methane has initially been reported Webster et al. (2013b). After the acquisition of further data and a reprocessing of the entire dataset, spanning 605 sols, the MSL Science Team has found a mean methane value of 0.69 ± 0.25 ppbv. They have also found a strong and rapid variability with a high methane abundance period, lasting at least 60 sols, during which the mean value was 7.19 ± 2.06 ppbv, decreasing to the mean value again only 47 sols later (Webster et al. 2015).

However, the comparison between these values and the derived methane content reported in the literature (10 to 15 ppbv on average, with a reported variability between 0 and 60 ppbv) is not straightforward. In fact, TLS has sampled a single location (roughly 1 km in size) and only the very lowest part of the Martian atmosphere (about 1 m), while in all the other cases (from both orbiting instruments and Earth-based observations) methane has been detected by integrating the entire thickness of the atmosphere and over broad regions of the surface, providing the ability to analyze and compare almost simultaneously data from places that are located far apart on the planet surface. Clearly, assessing the distribution and variability of methane in the Martian atmosphere might be very important to constrain the possible production and destruction mechanisms. It is therefore mandatory to understand any possible weakness in the measurement techniques that have led to methane detection on a planetary scale.

Several methane production mechanisms have been discussed in the literature, including 1) volcanism (Krasnopolsky et al. 2004); 2) meteoritic and/or cometary infall (Formisano et al. 2004); 3) H₂O-CO photolysis (Bar-Nun & Dimitrov 2007); 4) serpentinization (Oze & Sharma 2005); and 5) biogenic origin (Formisano et al. 2004; Krasnopolsky et al. 2004). The first three of these hypotheses have been excluded (Formisano et al. 2004; Ryan et al. 2006; Atreya et al. 2007) while the last two are still a matter of debate (Atreya et al. 2011). In hypotheses 4 and 5, liquid water must be involved.

An important point stressed by Geminale et al. (2008), Mumma et al. (2009), Fonti & Marzo (2010) and Geminale et al. (2011) was that the methane mixing ratio in the Martian atmosphere varies spatially and temporally. It is remarkable that the maximum content of methane is related to regions such as Nili Fossae, where serpentinization and carbonation have occurred (Ehlmann et al. 2009). Since these regions are interpreted as being among the oldest of the Martian surface (≈4 Ga old), a possible temporary storage reservoir has been invoked in the form of methane clathrate particles (Chastain & Chevrier 2007; Chassefière 2009). Since these clathrates are suggested to be present in the subsurface and at the base of the cryosphere, it is possible that they could be transported to the surface, releasing methane due to some destabilization process related to the Martian seasonal cycle (Atreya et al. 2011; Mousis et al. 2013).

The temporal variability points out the presence of a very efficient loss mechanism for Martian methane. UV photolysis and oxidation destruction processes of CH₄, by O(¹D) and OH, are extremely efficient and imply a CH₄ lifetime of 300−600 years (Atreya et al. 2007) in contrast with the rapid reduction of CH₄ associated with observations. To explain the more rapid loss of CH₄ implied by the observations, several mechanisms have been analyzed, and the most probable currently seems to be the heterogeneous oxidation of CH₄ by the Martian surface (Atreya et al. 2007; Lefèvre & Forget 2009; Mumma et al. 2009). This implies a reactive surface that could remove methane much more efficiently than the photochemical loss (200 days, Lefèvre & Forget 2009). The most probable oxidant that could render the Martian surface reactive is hydrogen peroxide (H₂O₂), which was detected in 2003−2004 with a mixing ratio of 20−40 ppbv (Encrenaz et al. 2004; Clancy et al. 2004). Another strong oxidant is perchlorate salt, which was discovered at the Phoenix landing site (Hecht et al. 2009) and has been proposed to cause the oxidation of methane, but laboratory measurements have shown that neither ClO⁴⁻ salts nor H₂O₂ alone could explain the observed temporal and spatial variability of atmospheric methane (Gough et al. 2011).

In this context, the work of Fonti & Marzo (2010) is particularly important. Using a statistical clustering technique, they analyzed spectra acquired with the Thermal Emission Spectrometer (TES), spanning three Martian years (MY). The spectral range covered by TES (200–1700 cm^-1; 50–6 μm) allowed them to study, for the first time after the attempt of Maguire (1977), the second strongest methane band, the ν₄ Q-branch at 1304 cm^-1 (7.7 μm), while all the other analyses concentrated on the stronger 3018 cm^-1 (3.3 μm) ν₃ Q-branch band system. The disadvantages due to the rather low resolution of the instrument (12.5 or 6.25 cm^-1) have been compensated for by averaging a large number of spectra (about three million, in total). One of the most important result found by Fonti & Marzo (2010) is that at the solar longitude (Ls) values investigated (i.e., 0°, 90°, 180° , and 270°), the methane content exhibits a seasonal variability cycle and a recurrent spatial distribution with higher values over Tharsis, Arabia Terrae, and Elysium.

The detected variability, both temporal and spatial, is definitely intriguing, but none of the attempts to explain the observations have been fully satisfactory. As an example, we report the preliminary analysis we made to understand whether the methane temporal variability, found by Fonti & Marzo (2010), follows a regular cycle, with possible correlations with other atmospheric cycles, already well understood. To this purpose, the global methane content was evaluated as a function of time and the trend was fit with a second-order truncated Fourier series, for example, Eq. (1), to reflect the fact that water vapor (and maybe methane) have two maxima and two minima in each year: (1)where the values of the five parameters a₀, a₁, a₂, b₁, and b₂ were adjusted using the χ² technique to obtain the best fit. A similar analysis was performed for the aerosol opacity and the water vapor cycle on the data obtained with the same TES spectrometer in the same period of time. These same data have been analyzed by Smith (2004), who reported the temporal and spatial variation of the dust optical depth and the water vapor column abundance at 14:00 local time, as a function of Ls, for the same three MYs considered by Fonti & Marzo (2010). The original data have been averaged in latitude, therefore only the temporal variations of the global value of the two quantities were considered. The results of the three fits are shown in the three panels of Fig. 1. This figure, as well as Figs. 2 and 4, were modified from Mancarella et al. (2014).

Fig. 1 Seasonal variation of the global methane content in the Martian atmosphere (top panel, points from Fonti & Marzo 2010) compared to the data (points) analyzed by Smith (2004) for atmospheric dust (central panel; the dust optical depth is evaluated at 1075 cm-1) and water vapor (bottom panel). The data, fit with a double sinusoidal function (solid lines), seem to suggest a possible methane cycle, with possible phase correlations to dust and water vapor cycles. The two red vertical lines mark the maxima (solid) and minima (dashed) of the fit methane. The determination coefficient R2 is given for each fit, R2 = 1 corresponding to the perfect fit.

The phase correlations shown in Fig. 1 between the methane-interpolated function and the two other trends might suggest possible interpretations of the methane variability, but, unfortunately, the information retrievable from this analysis might be premature because of the low temporal sampling of the global methane content. In turn, the effort needed to improve the temporal coverage would not be justified without a validation, leading to an improved confidence in the result obtained by Fonti & Marzo (2010). We here discuss our effort to improve the confidence in the analysis reported by Fonti & Marzo (2010) of the TES data, and thus enable higher temporal sampling of methane content.

In the following section we review the statistical clustering used by Fonti & Marzo (2010), focusing on the weak points of the method. At the same time, we describe and discuss the improvements we implemented in their procedure in the attempt to address the weak points in the analyses of Fonti & Marzo (2010). At the end of the section we also introduce an alternative method for the analysis of the data, without using the clustering technique. In the next section we describe and discuss the result obtained after implementing all the improvements. Finally, in the last section we draw some conclusions.

2. Methodology

2.1. Statistical clustering approach

The use of a statistical clustering technique by Fonti & Marzo (2010) to analyze TES data was justified by the large number of TES spectra used to obtain a good signal-to-noise ratio by averaging data at each Ls investigated. The unsupervised clustering technique assigns each spectrum of the dataset to one of the clusters, which eliminates bias from the procedure. More detailed information on this technique can be found in Marzo et al. (2006) and references therein, while a successful application to planetary remote sensing spectral data is described in Marzo et al. (2008).

The clustering procedure requires defining the selection criterion used for sorting the spectra into various clusters. This is a critical step, since the criterion allows the unsupervised sorting procedure to create clusters of spectra sharing similar behavior. Fonti & Marzo (2010) restricted the clustering to a narrow spectral region where the methane band is expected in the TES data at 1304 cm^-1. However, the resulting clusters might not be completely unambiguous because a range of unassociated spectral behavior may be mixed within a single cluster.

This ambiguity may arise from the unavoidable presence of anomalous spectra, such as those shown in Fig. 2. The peculiar behavior of these spectra is usually linked to temporary instrumental artifacts, and if their number is relatively large, they could affect the clustering procedure. Other spectra that could affect the clustering are those obtained at high latitudes, with a noise level much higher than the average due to the lower surface temperatures. It is therefore very important to remove anomalous and noisy spectra, so that the clustering can be applied to a clean, even if reduced, dataset.

Fig. 2 Two examples of anomalous spectra belonging to the dataset Ls = 180° ± 5° of MY24. One of the spectra (black line) has emissivity values increasing constantly well above 1 in the high wavenumber region, while the other (red line) shows a sudden drop in emissivity at about 1200 cm-1 and a systematic ripple (inset: see text) between 800 and 1200 cm-1.

2.2. Spectra selection

Fonti & Marzo (2010) attempted to avoid noisy and problematic spectra by considering only those taken in the best possible conditions. They imposed the following constraints on the original selection from the TES database: 1) latitude between −60° and +60°; 2) local time: 11:00–15:00; 3) emission angle ≤5°; and 4) ±2° around a specified Ls. We have adopted the same constraints in selecting the spectra, except for increasing the Ls interval to ±5°. This implies that we are able, at best, to provide qualitative comparisons with Fonti & Marzo (2010), but have a larger number of spectra in our original dataset.

2.3. Preprocessing of spectra

In addition, Fonti & Marzo (2010) used a simple two step preprocessing to remove further possibly problematic spectra from the selected TES dataset. First, they discarded all the spectra with at least one emissivity value outside the interval 0–2. Second, they selected 80 spectral channels, including both the CO₂ band and the methane band spectral regions, and applied an initial clustering to the data in the range 560–1410 cm^-1. The result was a large number of clusters, most of them containing only a few spectra, with a single cluster containing the vast majority of the spectra of the original dataset. This last cluster was their preprocessed dataset, on which they applied further analysis.

However, this approach may not be particularly effective in removing the noisy spectra and might even not completely remove anomalous spectra. Their presence may result in a large uncertainty in the derived methane content. For this reason, we have developed a more elaborate three-step preprocessing procedure aimed at removing the vast majority of the anomalous and noisy spectra, while retaining the greatest number of spectra in the dataset to be analyzed for CH₄. The preprocessing procedure was only applied to the data between 200 and 1600 cm^-1, excluding the first five channels of each spectrum because they are null, and the last ten channels because of their great variability. Below we provide details of the preprocessing procedure, showing examples, all of which, including those in Fig. 2, refer to the dataset of Ls = 180° ± 5° of MY24, containing 977 561 spectra, not quite doubling the 516 701 spectra of Fonti & Marzo (2010), where they used ±2° of Ls.

2.3.1. Upper and lower limit

Similar to Fonti & Marzo (2010), but more restrictive, in the first step we removed all spectra with at least one emissivity value in any spectral channel outside the range 0.05−1.10. The lower limit was defined after confirming that the minimum of the CO₂ main band at 670 cm^-1 was not affected. The upper limit was defined after confirming that only anomalous spectra were removed from the dataset. We chose an emissivity value greater than unity, since the calibration procedure of the TES radiance data affects the spectral region around the methane band, where some emissivity values might slightly exceed unity, even for spectra without any anomaly. The number of spectra remaining after this step is 843 625 (about 86% of the initial dataset).

Fig. 3 Histogram showing the number of spectra vs. the RP value. The values of the mean and of the standard deviation of the distribution are reported.

2.3.2. Ripple removal

As shown in the inset of Fig. 2, some spectra exhibit a regular peak-to-peak variability, hereafter ripple, with a periodicity of two spectral channels (in the example shown in Fig. 2 in the 800−1200 cm^-1 range). This could be very dangerous for our analysis because the methane band has exactly the same periodicity, although limited to a single period. We tried to remove the spectra with the ripple in the methane band spectral region. To address this problem, we defined a ripple parameter (RP), given by (2)where ∑ _evenϵ and ∑ _oddϵ are the sum of the emissivity values in the even and odd channels, respectively, of the spectral range 1200−1400 cm^-1. If the emissivity values are randomly distributed, without a two-channel periodicity, then the value of the RP would tend to one. If a ripple with a large enough amplitude is present, then the RP value will be higher or lower than one by an amount depending on the ripple amplitude. Figure 3 shows a histogram of the number of spectra vs. the RP value. The distribution is generally symmetrical around the peak value of about 1.002. For further analysis we retained all the spectra within ±0.01 of the mean value (0.992 < RP < 1.012). The acceptance interval corresponds to roughly ±2σ and only removes the spectra in those bins with <10% of the number of spectra in the most populated bin. This is a good compromise between removing most problematic spectra without significantly reducing the number of selected spectra. After this second step 782 478 spectra remain (80% of the initial dataset).

The result of steps 1 and 2 of the pre-processing is shown in Fig. 4, where the average of the spectra before (top panel) and after (bottom panel) eliminating data are shown along with their associated standard deviations. The difference between the standard deviation of the two spectra is particularly significant between 1200 cm^-1 and 1600 cm^-1 and quite consistent between 600 cm^-1 and 1200 cm^-1.

Fig. 4 Top panel: average and standard deviation of the 977 561 TES spectra belonging to Ls = 180° ± 5° of MY24 without any preprocessing. Bottom panel: average and standard deviation of the remaining 782 478 TES spectra after applying the two initial preprocessing steps described in the text.

Fig. 5 Cartoon illustrating how the noise parameter is defined. a) Three notional noisy spectra are shown together with the notional grand average of all the entire hypothetical dataset, assumed to be noiseless. Example 1 (circles) is basically superimposed on the grand average; example 2 (triangles) is displaced, but with roughly the same trend of the grand average; and example 3 (squares) has a trend different from the grand average. b) The difference between each of the three spectra, shown in panel a) and the grand average; the three straight lines represent the linear fit of each of the three differential spectra. On the right side are shown the three differences between the linear fits and the data shown in panel b), relative to example 1 (panel c)), example 2 (panel d)), and example 3 (panel e)), respectively.

2.3.3. Noise and temperature cutoff

In spite of the dramatic improvement, particularly around 1300 cm^-1, where the methane band is located, we have certainly not yet removed most of the noisy spectra from the original dataset. To select only the lowest noise spectra, it would be beneficial to associate a noise parameter with each spectrum and select the spectra as a function of this parameter. Unfortunately, defining a noise parameter for TES emissivity spectra is not straightforward. This is because TES spectra are collected in radiance and are later converted in emissivity to facilitate the comparison with each other, independently of the surface temperature of Mars during the observation.

The noise-equivalent radiance (NER) of TES is stated in Christensen et al. (2001): it corresponds to 1-sigma radiance noise level of an individual spectral sample and is 2.5 × 10^-8 W cm^-2 str^-1 (cm^-1)^-1 in the spectral range between 300 and 1400 cm^-1. Deriving from this radiance noise value the corresponding emissivity noise value, to be associated with each of the emissivity spectra in our dataset, requires taking into account 1) the wavenumber range; 2) the observed target temperature (the surface temperature of Mars, in this case); and 3) the theoretical relation between radiance and temperature given by Planck’s law. The wavenumber range of interest is 1200−1400 cm^-1, around the methane band at 1304 cm^-1 and within the spectral range 300–1400 cm^-1 cited by Christensen et al. (2001). The surface temperature of Mars has been evaluated by the TES team for each spectrum and can be found in the TES dataset¹.

It is therefore possible to associate an expected emissivity noise value with each emissivity spectrum, based on the NER and the surface temperature of Mars associated with that spectrum. In principle, it would be possible to select the lower noise spectra by simply rejecting all the spectra with a surface temperature below a chosen value, but that would not take into account anomalous spectra affected by a noise higher than that expected for the associated surface temperature. Hence the necessity to define a noise parameter for each TES spectrum, based only on the spectral behavior in the spectral region between 1200 and 1400 cm^-1.

To do this, we first calculated the average spectrum for the 782 478 left after the two pre-processing steps described above. A cartoon representation of this grand average is shown as the thick solid line in Fig. 5a. Then, for each of the 20 spectral channels between 1200 and 1400 cm^-1, we calculated the difference, difference 1, between each spectrum and the average spectrum, as illustrated in Fig. 5b. Since the noise of the grand average is much lower than that of any single spectrum, we assumed it to be essentially noiseless and chose as the noise parameter of a given spectrum the standard deviation of the 20 difference values evaluated for that spectrum. However, a slight offset, or, more important, a slope difference, examples 2 and 3 in Fig. 5a, would result in an overestimation of the noise parameter associated with each spectrum, as shown by the data points in Fig. 5b. In other words, the noise parameter to be associated with each spectrum would also have included the effect of the variability among different spectra and not only that within each spectrum.

Fig. 6 Noise parameter vs. temperature scatter plot: each black dot represents one spectrum. The red curve represents the expected noise as a function of temperature, obtained with the theoretical relation between the surface temperature and the given radiance noise, while the green crosses mark the average of the noise parameter distribution. The two blue lines denote the noise and temperature limits used for spectra selection.

To more accurately represent the noise of a single spectrum, we calculated the linear fit, between 1200 and 1400 cm^-1, of the 20 difference values, as shown in Fig. 5b. To take into account the effect of the inherent variability among different spectra, for instance, examples 2 and 3 in Fig. 5a, we calculated the difference, difference 2, between each of the linear fits, shown in Fig. 5b , and the corresponding difference of the spectra to the grand average (difference 1, also shown in Fig. 5b). Difference 2 values are comparable for the three examples shown in Figs. 5c−e, in spite of the different slope of fit of example 3 in Fig. 5b. The standard deviation of the 20 values of difference 2 is representative of the variability within each spectrum, in the region of interest, and was chosen as the noise parameter associated with each spectrum. A noise parameter defined in this way and associated with each TES emissivity spectrum, is consistent with the NER definition given in Christensen et al. (2001) for the radiance spectra.

The result of this procedure is summarized in Fig. 6, where the distribution of the 782 478 spectra of our dataset is shown as a function of noise vs. temperature. Each black dot corresponds to a spectrum. The distribution of the spectra in this plot can be easily compared with the expected emissivity noise values, calculated using the theoretical relation between NER and surface temperature, as represented by the red curve. The green crosses mark the trend in the scatter plot of the average noise parameter values, obtained by binning the data in 32 temperature intervals, each 2.5 K wide. As can be seen, the statistical distribution of our derived noise values for each temperature are remarkably close to the theoretically expected value for that temperature, at least where the number of spectra is sufficiently high to allow meaningful statistics.

Now we defined and applied the third step of the preprocessing procedure. If the distribution of the spectra in Fig. 6 were much narrower (i.e., basically coincident with the red curve), it would be possible to choose for the cutoff either a noise value or the corresponding temperature value. Since the real spectra are distributed in a rather wide strip around the red curve, it is better to select the spectra using the cutoff values in both noise and temperature. As a result of the relation between the two quantities, we can set the value of each of the two variables and evaluate the corresponding value of the other variable.

To retain the maximum number of spectra and at the same time minimize the number of noisy spectra, we decided to select only those spectra with a surface temperature ≥250 K and an emissivity noise parameter ≤0.017. The two blue lines in Fig. 6 indicate these two values, and it can be seen that the green crosses begin to deviate from the theoretical trend just left of their crossing point, which is obviously exactly on the red curve. The two lines divide the scatter plot into four quadrants, and we selected for the subsequent analyses only the spectra lying in the bottom right quadrant: their number is 710 874 (about 73% of the initial dataset).

The effects of the various steps of the preprocessing procedure are summarized in Table 1.

2.4. Clustering criterion

As pointed out at the beginning of this section, the selection criterion is crucial for a successful application of the clustering analysis. In Fonti & Marzo (2010) the selection criterion used for statistical clustering was based on the emissivity values in the methane channel (TES channel #110 at 1304 cm^-1) and in the adjacent channel at lower wavenumbers (TES channel #109 at 1294 cm^-1). This criterion, however, might sometimes produce unexpected selections, as in the example shown in Fig. 7.

Fig. 7 Example of the application of clustering, using the criterion of Fonti & Marzo (2010), on the 782 478 spectra of the Ls = 180° ± 5°, MY24, dataset, after the first two steps of the preprocessing: the two spectra represent the average of the spectra contained in each cluster. Inset: enlargement of the spectral range between 1200−1400 cm-1.

In Fig. 7 the different behavior of the two spectra in the whole spectral range is clearly visible. In particular, there is a clear difference in the emissivity level at low wavenumber and in the depth of the wide band, around 1050 cm^-1, which is due to atmospheric dust. These two effects might be due to a possible connection between the two emissivity values used for the clustering procedure and some spectral characteristic, influenced by the spectral behavior at these wavenumbers, such as the surface temperature and/or emissivity or the load of atmospheric dust. This is confirmed by the trend of the two spectra in the inset of Fig. 7, where we show an enlargement of the spectral region near the methane band. The two spectra cross each other very close to the spectral channels whose emissivity values were chosen for the selection criterion. The anomalous behavior may also be linked to the fact that this is the spectral region where the TES radiance spectra are normalized to obtain the corresponding emissivity spectra.

To address this problem, we decided to adopt a new selection criterion based on the value of the methane band depth (BD), defined as in Clark & Roush (1984), (3)where b is the emissivity value in the channel corresponding to the methane band (TES channel #110 at 1304 cm^-1) and c is the average of the two emissivity values on each side of the band (TES channel #109 at 1294 cm^-1 and TES channel #111 at 1314 cm^-1).

The result of the clustering procedure, applied to the same dataset as in Fig. 7, using the BD is shown in Fig. 8. The spectra of the two clusters have almost the same overall behavior for the entire spectral range and differ chiefly in the region of the methane band (inset), where the effects of the crossing of the two spectra are greatly reduced. The residual difference might be due to the presence of remaining noisy spectra, to be removed with the third step of the preprocessing.

Fig. 8 Example of the application of clustering for the same dataset as in Fig. 7, using the BD as clustering criterion. Inset: enlargement of the spectral range between 1200–1400 cm-1.

The BD clustering criterion, although it basically solves the problem of selecting spectra by some characteristic(s) other than the methane content, apparently does not reduce the spectral variability in the three selected spectral channels at all. This is probably because the BD clustering criterion cannot distinguish between a low (or high) emissivity value in each of the three channels that is due to the presence of methane and the same behavior that is due to the effect of spectral noise. In other words, the spectra with an emissivity value in the central channel lower than those in the other two channels will be assigned to the same cluster, whether this is due to the methane band or to the random distribution of the emissivity values linked to the residual spectral noise. The net result is that the contribution of methane to the variability is dwarfed by the effect of this random noise.

2.5. Alternative to clustering

An alternative solution to the clustering procedure would be trying to detect the presence of methane by comparing each average spectrum with a tailored synthetic spectrum, where the amount of methane can be varied. This direct approach, however, has already been unsuccessfully tried by Fonti & Marzo (2010). This failure was due to the difficulty of reproducing, with the required accuracy, the average of a very large number of spectra having a wide range of variability both of the surface conditions and of the atmospheric state. The possibility of a correct reproduction of the resulting average spectrum depends both on the capability of developing a suitable theoretical model for synthetic spectrum and on the absence of too many anomalous and/or noisy spectra in the observational data. Hence the need of the pre-processing procedure able to supply the best possible set of spectra to the final step of the methane detection procedure. In the next section, after briefly describing some important features of the clustering procedure adopted by Fonti & Marzo (2010) , we report and discuss the result obtained using all the improvements described and discussed in this section.

3. Results and discussion

Before discussing the results obtained by applying the clustering procedure to the preprocessed dataset of TES data (last row of Table 1), it is important to review the procedure used by Fonti & Marzo (2010) to derive the methane amount at each Ls examined.

The first step in estimating a CH₄ abundance at each Ls was assessing the natural number of clusters, using the Caliński & Harabasz (1974) criterion (Marzo et al. 2006). For all the datasets the value of the Caliński & Harabasz (1974) parameter reached a maximum for two clusters, suggesting the natural presence of two clusters. Fonti & Marzo (2010) associated the cluster with the lower emissivity value in the methane band channel with the presence of methane. Then Fonti & Marzo (2010) derived a methane parameter, for each Ls from the ratio between the average of all spectra and the average of the cluster associated with the methane-free spectra. These methane parameters were then calibrated, using the methane upper limit of Maguire (1977; Fonti & Marzo 2010), and the quantity of methane present at each Ls was eventually retrieved.

When we applied the criterion of Caliński & Harabasz (1974) to the pre-processed dataset in the current study, we found the natural number of clusters less clearly defined. The result of this analysis is shown in Fig. 9 where the almost monotonic increase of the Caliński & Harabasz (1974) parameter seems to suggest that no cluster structure is present.

Fig. 9 Value of the Caliński & Harabasz (1974) parameter as a function of the number of clusters. The trend generally increases monotonically, suggesting no natural number of clusters. A slight decrease of the Caliński & Harabasz (1974) parameter is present after clusters 7 and 11.

Fig. 10 Result of the clustering in two clusters of the preprocessed dataset: each spectrum is the average of all the spectra assigned to that cluster by the clustering procedure and is identified by the number of spectra present in the corresponding cluster. Inset: ratio of each of the two spectra with the average spectrum of the entire dataset, in the spectral range 1200−1400 cm-1. The red straight line corresponds to the unity value.

We explored all the possible alternative choices at local maxima, dividing the pre-processed dataset into 7 and 11 clusters. We finally also forced the division into two clusters to provide a comparison to Fonti & Marzo (2010). The discussion of the results is better done by focusing on Fig. 10, where we show the output of the clustering of the preprocessed dataset into two clusters, as an example representative of every clustering alternative we have attempted.

Each of the two spectra shown in Fig. 10 is the average of all the spectra assigned to that cluster. The number of spectra present in each cluster is given in the legend, and it is almost the same: the difference between the two clusters is smaller than 5% of the total number of spectra in the entire dataset, suggesting a rather homogeneous distribution.

In the inset of Fig. 10 we show in the spectral region between 1200 and 1400 cm^-1 the ratio between each of the two spectra and the average of the entire preprocessed dataset. Showing the ratio can simplify the discussion of the results shown in Fig. 10, since Fonti & Marzo (2010) used a similar ratio as evidence for the presence of methane and from this derived its amount in the Martian atmosphere. The problem with the spectral ratios, shown in the inset of Fig. 10, is that they are fully compatible with the selection criterion imposed, but the amplitude of the oscillations in the region of the methane band is orders of magnitude larger than that expected for the possible variation of the methane content in the Martian atmosphere. In addition, they are almost symmetrical with respect to the unity value of the ratio (red line of Fig. 10), suggesting that the clustering might have been driven by some spectral parameter with a continuous variability much wider than that reasonably associable with the methane amount. On the other hand, the symmetrical distribution of the ratios, as well as the constantly increasing trend of the Caliński & Harabasz (1974) parameter in Fig. 9, seems to point to a residual strong influence of the spectral noise.

Another important consideration in interpreting the inset of Fig. 10 is that the two spectral ratios cross each other very close to the spectral channels where the TES radiance spectra are normalized to obtain the corresponding emissivity spectra, as already mentioned when discussing the spectra shown in Fig. 7. This occurrence can be linked to the residual temperature differences among the spectra and/or differences in surface emissivity or atmospheric opacity. The clustering procedure could have been driven more by such differences than by the different methane content. This explanation could also be linked to the variability shown by the spectra in the dust band (1000–1200 cm^-1), clearly visible in Fig. 10.

Unfortunately, the conclusion of this effort is that, in spite of all attempts, we have been unable to produce distinct methane and methane-free clusters and, as a result, we cannot use the same method as Fonti & Marzo (2010) to estimate the methane abundance. This led us to consider comparing the average spectrum of the preprocessed dataset with a synthetic spectrum, and to this purpose, we modified an existing procedure used to model the terrestrial atmosphere (Amato et al. 2002). Details of this procedure are described in Liuzzi et al. (2013) and Liuzzi et al. (2015), and we summarize here only the main features and the strengths of this approach.

The main characteristic of the procedure is the capability of retrieving simultaneously all the main geophysical parameters affecting the planetary radiance as observed by the orbiting TES instrument. The parameters are surface temperature and emissivity, atmospheric temperature profile, dust and ice mass-mixing ratio vertical profile. The water vapor columnar load is the only parameter retrieved separately from the others, since it produces only minor, though visible, features in the spectra. In this application we also tuned the methane columnar amount separately, aiming for a better fit of the radiance in the range 1280−1320 cm^-1, where the methane band is located. The procedure is a fully integrated retrieval framework, consisting of two modules: a forward one, the radiative transfer module, and an inverse module, which iteratively solves the radiative transfer equation.

The forward model is a monochromatic radiative transfer code solving the discrete radiative transfer equation (including reflected and solar radiance contributions) on a fixed pressure grid of 41 atmospheric layers, producing an analytical parametrization of the variable that allows accurate and fast calculations of the radiance, and of its Jacobians, which are the first derivatives, with respect to every variable. Radiance and Jacobians are calculated at a resolution of 10^-2 cm^-1 and then convolved to TES resolution.

Fig. 11 Top panel: comparison of the average TES spectrum (red line), of the pre-processed dataset, with the best-fit synthetic spectrum (black line). Bottom panel: comparison between the residual brightness temperature differences between the two spectra shown in the top panel and the expected TES NEdT for a brightness temperature of 280 K.

The inverse retrieval procedure is an adaptation to Mars of a method already extensively validated for Earth remote sensing (Carissimo et al. 2005). The inverse solution of the radiative transfer equation is achieved through a Gauss-Newton iterative scheme, which unfortunately is not straightforward to apply to the TES data. As an example, we have not been able to consider the actual variability of the surface emissivity, and we expressed it as a linear combination of three different spectral library members (Bandfield & Smith 2003) that are representative of the most common soil types in the portion of the planet we selected for our analysis (±60° in latitude).

The result obtained by comparing the average spectrum of the preprocessed dataset with the best-fitting synthetic spectrum is shown in Fig. 11, and the corresponding best-fit parameters are reported in Table 2. The very good agreement achieved in most of the spectral range is visually evident in the top panel of Fig. 11. It is quantified in the bottom panel, where the residual brightness temperature differences between the two spectra are compared to the expected TES noise-equivalent delta temperature (NEdT). The NEdT was evaluated, starting from the values of TES NER in the whole spectral range (Christensen et al. 2001), and using a procedure very similar to that described in Sect. 2 to calculate the emissivity noise. The black line in the lower panel of Fig. 11, representing the temperature difference, is generally within the two red lines, corresponding to the expected spectral noise, for most of the wavenumber range. Even in the very critical spectral region around the main CO₂ absorption band, at 660 cm^-1, the agreement is acceptable for our purposes and provides confidence in the robustness of the modeling procedure adopted.

More important is the result shown in Fig. 12, which is an enlargement of Fig. 11 in the spectral region 1200−1400 cm^-1. In Fig. 12, the average TES spectrum (red line) of Fig. 11 is compared with two synthetic spectra. One has an atmospheric methane content of 0 ppbv (black line), and the second of 33 ppbv (green line), the average amount of methane associated by Fonti & Marzo (2010) with their MY24 Ls 180° dataset.

Fig. 12 Comparison in the methane band spectral region of the average TES spectrum of the preprocessed dataset (red line) with two synthetic spectra that differ only in the amount of atmospheric methane: 0 ppbv (black line) and 33 ppbv (green line).

In Fig. 12 the two synthetic spectra are perfectly superimposed on each other with the exception of the channel corresponding to the methane band (#110 at 1304 cm^-1, indicated by the arrow in Fig. 12). The spectrum containing zero methane (black line) is slightly higher than the synthetic spectrum containing 33 ppbv methane (green line) and the average TES spectrum (red line). Unfortunately, in spite of the intriguing fact that the average spectrum seems to agree better with the synthetic spectrum containing 33 ppbv of methane, nothing definitive can be concluded due to the discrepancy between the synthetic spectra and the TES average over most of the spectral region shown in Fig. 12. The difference between the synthetic and measured spectra in the methane band, ≤0.1 K, is much smaller than those outside the methane region, ≈0.5 K. This probably implies that the degree of precision we need to convincingly assess the presence of methane in the Martian atmosphere is extremely hard to reach, if even possible, at present. None of the techniques and procedures discussed here (clustering, preprocessing, modeling of synthetic spectra) can be improved to the required level, at least without a major effort and a long and demanding commitment.

4. Conclusions

We have thoroughly revisited the procedure used by Fonti & Marzo (2010) in their statistical analysis of TES spectra, with the aim of better understanding whether the reported seasonal variability of the global methane content follows a regular cycle and if this cycle has any correlation with the seasonal cycles of other atmospheric components, such as water vapor and/or dust. We started with a revision of the clustering procedure used by Fonti & Marzo (2010) in the attempt of removing possible sources of ambiguity in the data selection, and we focused on two directions: 1) reducing the number of anomalous spectra in the dataset without significantly reducing the total number of spectra; and 2) revising the clustering selection criterion.

We applied our approach to the dataset of Ls = 180° ± 5° of MY24, containing initially 977 561 spectra, then reduced to 710 874. For this purpose, we used a three-step preprocessing protocol with the aim of defining a standard procedure for the selection of the less noisy spectra of each dataset. Unfortunately, the application of the new selection criterion, based on the definition of band depth and applied to the reduced dataset, did not result in an unambiguous identification of one or more clusters, characterized by the presence of the methane adsorption band at 1304 cm^-1.

This result can be due to a combination of several reasons, including, but not limited to 1) the inadequacy of the selection criterion; 2) the inadequacy of the preprocessing; 3) the inadequacy of the dataset used for the test; 4) the low spectral resolution of the TES spectra; and 5) the difficulty of extracting the signal variations linked to the methane band from the much higher noise level affecting the single spectrum, by simply averaging on a large number of spectra.

Of these, 2) and 3) most probably marginally influence the result, while the remaining reasons together, probably with different weights, prevent obtaining a clearly interpretable output from the clustering procedure. In any case, none of the above reasons is inherent to the possible presence and variability of methane in the Martian atmosphere, and therefore, at least in principle, our results neither contradict nor support the work of Fonti & Marzo (2010). This can be linked to the fact that, apparently, the disadvantages due to the change of the selection criterion outweigh the advantages obtained with the preprocessing, namely the removing of the most noisy spectra and the reduction of the overall noise level of the preprocessed dataset. The new selection criterion has indeed considerably reduced, although not eliminated, the ambiguity present in the previous selection criterion, exemplified in Fig. 7. Unfortunately, it has at the same time increased that same ambiguity by strongly enhancing the contribution of spectral noise and, as a consequence, made the process of assigning each spectrum to the appropriate cluster almost impossible.

We directly compared the average spectrum of the preprocessed dataset with a best-fit synthetic spectrum, which has been unsuccessfully attempted by Fonti & Marzo (2010). The result was encouraging, with the difference between the two spectra, on average, below the instrumental noise of a single TES spectrum. Unfortunately, the result is not sufficient to unambiguously identify the presence of methane. This is due to the inherent variability of the TES spectra and the extreme difficulty of reducing the variability below a level comparable with the effect of methane on the spectrum. The presence of a tiny amount of methane is visible when comparing the two synthetic spectra, which differ only for the methane content (Fig. 12). Obviously, this difference cannot be considered a confirmation of the presence of methane in the Martian atmosphere because of the much larger differences between the observed and the synthetic spectrum in most of the spectral range shown in Fig. 12.

Our findings confirm what was found, but not explicitly mentioned, by Fonti & Marzo (2010), which is that no synthetic spectrum can reproduce an average Martian spectrum with the precision required to assess the presence of methane.

It is clear that there is ample space for future work on the problems of using TES data to assess the Martian methane question. These range from changing the selection criterion to using a partially supervised clustering procedure, from limiting the spectral variability by reducing the range of the chosen spectra, to using only TES data obtained with a 6.25 cm^-1 spectral resolution. It is even possible, although not easy, to improve the modeling procedure aimed at reducing the uncertainty level of the fit. The probability of reaching definite results about the presence and variability of the Martian methane by implementing the mentioned improvements into the procedure developed so far, is probably rather low, compared to the required effort.

This consideration about the opportunity to continue the work described in this paper is considerably strengthened by the very recent detection of methane by TLS of SAM on the Curiosity rover (Webster et al. 2015). The near future holds the promise for acquiring more definitive evidence about the Martian methane question from measurements of orbiting instruments designed with the declared purpose to detect methane in the Martian atmosphere with very high accuracy and reliability. The Methane Sensor for Mars (MSM) onboard the Mars Orbiter Mission (known as Mangalyaan, an Indian mission) has reached its final orbit around Mars in late September 2014 (Sundararajan 2013). After this, the ESA mission ExoMars, scheduled to be launched in January 2016, will carry onboard the Trace Gas Orbiter (TGO). Among its tasks are observations of the Sun through the Martian atmosphere using two instruments, Nadir and Occultation for Mars Discovery (NOMAD) and the Atmospheric Chemistry Suite (ACS), which together cover a spectral range from the near-UV to thermal IR (Zurek et al. 2011).

As a final remark, we wish to stress that even though we have not reached our aim of improving the temporal sampling of the TES data beyond that of Fonti & Marzo (2010), the extensive work we have done in the attempt to improve their procedure has produced some interesting byproducts. In addition to extending to the Martian atmosphere the technique used to model the terrestrial atmosphere (Liuzzi et al. 2013, 2015), we developed new methods for the preprocessing of a dataset of planetary spectra, aimed at eliminating spurious data and selecting only high-quality spectra. These methods, summarized in Figs. 3, 5, and 6, can be applied in a wide range of different situations and may prove highly beneficial to researchers working not only in planetary spectroscopy, but also in other fields that require a preliminary selection of the best spectra from an extensive dataset.

Acknowledgments

S.F. and T.L.R. thank NASAs Visiting Research Associate program, at NASA Headquarters, for enabling SF’s visits to NASA Ames that were invaluable for the completion of this effort. This research has been partially supported by the Italian Space Agency (ASI) and the Italian Ministry of University and Research.

References

All Tables

All Figures

Fig. 1 Seasonal variation of the global methane content in the Martian atmosphere (top panel, points from Fonti & Marzo 2010) compared to the data (points) analyzed by Smith (2004) for atmospheric dust (central panel; the dust optical depth is evaluated at 1075 cm-1) and water vapor (bottom panel). The data, fit with a double sinusoidal function (solid lines), seem to suggest a possible methane cycle, with possible phase correlations to dust and water vapor cycles. The two red vertical lines mark the maxima (solid) and minima (dashed) of the fit methane. The determination coefficient R2 is given for each fit, R2 = 1 corresponding to the perfect fit.

In the text

	Fig. 2 Two examples of anomalous spectra belonging to the dataset Ls = 180° ± 5° of MY24. One of the spectra (black line) has emissivity values increasing constantly well above 1 in the high wavenumber region, while the other (red line) shows a sudden drop in emissivity at about 1200 cm-1 and a systematic ripple (inset: see text) between 800 and 1200 cm-1.
In the text

	Fig. 3 Histogram showing the number of spectra vs. the RP value. The values of the mean and of the standard deviation of the distribution are reported.
In the text

	Fig. 4 Top panel: average and standard deviation of the 977 561 TES spectra belonging to Ls = 180° ± 5° of MY24 without any preprocessing. Bottom panel: average and standard deviation of the remaining 782 478 TES spectra after applying the two initial preprocessing steps described in the text.
In the text

Fig. 5 Cartoon illustrating how the noise parameter is defined. a) Three notional noisy spectra are shown together with the notional grand average of all the entire hypothetical dataset, assumed to be noiseless. Example 1 (circles) is basically superimposed on the grand average; example 2 (triangles) is displaced, but with roughly the same trend of the grand average; and example 3 (squares) has a trend different from the grand average. b) The difference between each of the three spectra, shown in panel a) and the grand average; the three straight lines represent the linear fit of each of the three differential spectra. On the right side are shown the three differences between the linear fits and the data shown in panel b), relative to example 1 (panel c)), example 2 (panel d)), and example 3 (panel e)), respectively.

In the text

Fig. 6 Noise parameter vs. temperature scatter plot: each black dot represents one spectrum. The red curve represents the expected noise as a function of temperature, obtained with the theoretical relation between the surface temperature and the given radiance noise, while the green crosses mark the average of the noise parameter distribution. The two blue lines denote the noise and temperature limits used for spectra selection.

In the text

	Fig. 7 Example of the application of clustering, using the criterion of Fonti & Marzo (2010), on the 782 478 spectra of the Ls = 180° ± 5°, MY24, dataset, after the first two steps of the preprocessing: the two spectra represent the average of the spectra contained in each cluster. Inset: enlargement of the spectral range between 1200−1400 cm-1.
In the text

	Fig. 8 Example of the application of clustering for the same dataset as in Fig. 7, using the BD as clustering criterion. Inset: enlargement of the spectral range between 1200–1400 cm-1.
In the text

	Fig. 9 Value of the Caliński & Harabasz (1974) parameter as a function of the number of clusters. The trend generally increases monotonically, suggesting no natural number of clusters. A slight decrease of the Caliński & Harabasz (1974) parameter is present after clusters 7 and 11.
In the text

Fig. 10 Result of the clustering in two clusters of the preprocessed dataset: each spectrum is the average of all the spectra assigned to that cluster by the clustering procedure and is identified by the number of spectra present in the corresponding cluster. Inset: ratio of each of the two spectra with the average spectrum of the entire dataset, in the spectral range 1200−1400 cm-1. The red straight line corresponds to the unity value.

In the text

	Fig. 11 Top panel: comparison of the average TES spectrum (red line), of the pre-processed dataset, with the best-fit synthetic spectrum (black line). Bottom panel: comparison between the residual brightness temperature differences between the two spectra shown in the top panel and the expected TES NEdT for a brightness temperature of 280 K.
In the text

	Fig. 12 Comparison in the methane band spectral region of the average TES spectrum of the preprocessed dataset (red line) with two synthetic spectra that differ only in the amount of atmospheric methane: 0 ppbv (black line) and 33 ppbv (green line).
In the text