Mars Methane Sources in Northwestern Gale Crater Inferred From Back Trajectory Modeling.

During its first seven years of operation, the Sample Analysis at Mars Tunable Laser Spectrometer (TLS) on board the Curiosity rover has detected seven methane spikes above a low background abundance in Gale crater. The methane spikes are likely sourced by surface emission within or around Gale crater. Here, we use inverse Lagrangian modeling techniques to identify upstream emission regions on the Martian surface for these methane spikes at an unprecedented spatial resolution. Inside Gale crater, the northwestern crater floor casts the strongest influence on the detections. Outside Gale crater, the upstream regions common to all the methane spikes extend toward the north. The contrasting results from two consecutive TLS methane measurements performed on the same sol point to an active emission site to the west or the southwest of the Curiosity rover on the northwestern crater floor. The observed spike magnitude and frequency also favor emission sites on the northwestern crater floor, unless there are fast methane removal mechanisms at work, or either the methane spikes of TLS or the non-detections of ExoMars Trace Gas Orbiter cannot be trusted.

Introduction

Almost all of the methane in the present‐day Earth's atmosphere can be traced back to biological origins (Cicerone & Oremland, 1988). Extending this observation to Mars, the presence of atmospheric methane could be a biosignature on this seemingly lifeless planet (Yung et al., 2018). Alternative, abiotic methane production mechanisms on Mars invoke past or present geological activity (reviewed in Oehler & Etiope, 2017) such as serpentinization (Oze & Sharma, 2005), which would indicate the presence of liquid water, an indispensable ingredient for life. Abundant methane in the ancient Martian atmosphere could also provide a solution to the conflict between the Faint Young Sun and a warm surface suggested by fluvial and lacustrine features on Mars (e.g., Kite et al., 2017).

In the past two decades, the significance of methane in the Martian atmosphere has motivated a number of remote sensing observations aimed at both retrieving the methane abundance in the Martian atmosphere and mapping out its spatial distribution. These observations have reported inconsistent and highly variable methane concentrations (Aoki et al., 2018; Geminale et al., 2008; Giuranna et al., 2019; Fonti & Marzo, 2010; Formisano et al., 2004; Krasnopolsky, 2012; Krasnopolsky et al., 2004; Mumma et al., 2009). To overcome the technical challenges faced by remote sensing observations, the Tunable Laser Spectrometer (TLS; Mahaffy et al., 2012) on board the Curiosity rover was sent to Gale crater to make in situ measurements. During 7.1 years of operation as of January 2020, twenty direct‐ingest measurements and 16 enrichment measurements (refer to Webster et al., 2015, 2018 for the descriptions for the two measurement types) revealed a baseline level of ∼0.41 parts‐per‐billion‐by‐volume (ppbv) (Webster et al., 2018), with episodic spikes up to ∼21 ppbv (Webster et al., 2021) as summarized in Figure 1. These spikes have been interpreted as discrete, possibly nearby, methane emission events (Webster et al., 2015, 2018). Notably, the latest ∼21 ppbv methane spike was the first spike detected in an enrichment measurement. It has a much higher signal‐to‐noise ratio than the earlier spikes detected in direct‐ingest measurements and therefore has increased the credibility of the previous methane spike observations. Concurrently, the ExoMars Trace Gas Orbiter (TGO) has made solar occultation measurements of methane concentration at mid‐ to high‐altitudes. However, it has reported stringent upper limits down to 0.02 ppbv (Knutsen et al., 2021; Korablev et al., 2019; Montmessin et al., 2021). Assuming methane is a long‐lived species with a 330‐year lifetime as indicated by standard photochemical models (Lefèvre & Forget, 2009), it should be uniformly mixed throughout the Martian atmosphere, so TGO's stringent upper limits have been interpreted as an upper limit for methane concentration in the entire atmosphere, which is contradictory to TLS's significantly more elevated ∼0.41 ppbv background level. But some mechanisms have been proposed to reconcile this inconsistency. For example, TLS performed all its measurements in the near‐surface planetary boundary layer (PBL), and methane, if released from the surface, could accumulate in the shallow nighttime PBL (Moores, Gough, et al., 2019; Moores, King, et al., 2019). Some speculative fast removal mechanisms that can possibly cause temporal and spatial inhomogeneity of methane concentration have also been proposed (Gough et al., 2010; Hu et al., 2016; Knak Jensen et al., 2014), which might also reconcile the inconsistency between the TLS and the TGO results. In this study, we first accept both the results from TLS and TGO, and investigate the circumstances under which their discrepancies can be reconciled. We will then evaluate the likelihood of these circumstances.

Figure 1

Tunable Laser Spectrometer methane signals versus Mars season and local time. The seven data points above 5 ppbv are regarded as “methane spikes” with their indices labeled. The 29 data points below 5 ppbv are regarded as the background abundance. Two background‐level measurements are also marked, one performed immediately before the detection of Spike 6 and the other after Spike 7. Direct‐ingest measurements are shown in circles. Enrichment measurements are shown in squares. Colors show the local time of methane ingestions. Error bars show ±1 σ uncertainty. Adapted from Webster et al. (2018, 2021).

Assuming the existence of methane on Mars is real, its origin will have profound implications for geology and astrobiology. Identification of the methane's origin requires that we first find the surface emission sites, results from which can inform future missions of high priority landing sites and enable them to directly probe the methane source. The results can also guide orbiting instruments to better focus their methane observation strategies.

Inferring the locations of methane emission sites requires correct modeling of complex atmospheric transport processes. An early attempt to do so involved using a diffusion model to represent the spread of observed methane plumes (Mumma et al., 2009), which was shown to be oversimplified by addressing the importance of advection by bulk wind (Mischna et al., 2011). More recently, the Global Environmental Multiscale (GEM)‐Mars general circulation model (GCM) was used to simulate methane transport and then a statistical approach based on the idea of simultaneous satisfaction of multiple observational constraints was used for methane source localization (Giuranna et al., 2019). Results suggested an emission region to the east of Gale crater for TLS's first methane spike (Spike 1 in Figure 1). Later, the Mars Regional Atmospheric Modeling System mesoscale model was used to simulate the transport and dispersion of methane plumes emitted from 10 selected source regions around Gale crater (Pla‐García et al., 2019). Substantial dilution during tracer transport was observed, which demonstrates the importance of incorporating turbulent dispersion into tracer transport modeling. Among all the 10 emission region candidates, the region to the northwest of the crater was favored, different from the prior findings of Giuranna et al. (2019).

The aforementioned studies have all adopted a forward Eulerian approach to identify potential emission sites, in which the model integrates three‐dimensional tracer fields forward in time and quantifies how much tracer released at a specific emission location at a specific time can ultimately reach the detector. However, this trial‐and‐error approach is computationally inefficient, as most of the injected methane does not reach the detector, so usually only a small number of putative emission sites are selected and studied in depth (e.g., Pla‐García et al., 2019). Meanwhile, the spatial resolution of emission regions is limited by the size of GCM grid boxes. In commonly used GCMs that provide a global coverage, the size of grid boxes is typically a few degrees in latitude and longitude, or several hundred kilometers. Such GCMs have difficulty in differentiating emission sites within the 154‐km diameter Gale crater (e.g., Giuranna et al., 2019).

In this work, we adopt an inverse Lagrangian approach (Lin et al., 2003, 2012) of emission site identification to overcome the challenges faced by the forward Eulerian approach. The inverse Lagrangian approach, also known as back trajectory analysis, is widely used in the environmental science community to map out upstream emission regions (e.g., Gerbig et al., 2003; Kort et al., 2008; Lin et al., 2003; Lin et al., 2004; Macatangay et al., 2008; Mallia et al., 2015). An ensemble of computational particles, representing individual air parcels, is released from the detector at the time of detection and transported backwards in time within the model. The particles' transport pathways are determined by the bulk wind, and the particles are dispersed by parameterized subgrid‐scale turbulence. The locations where backward‐traveling particles are found within the PBL and, hence, are potentially affected by surface emission, are identified as potential upstream surface emission regions. The quantitative linkage between measured atmospheric mole fraction at the detector (e.g., TLS) and upstream surface fluxes can be established for any putative emission site via the number density of particles at that emission site (Fasoli et al., 2018; Lin et al., 2003). A single inverse Lagrangian simulation can quantify the influence of all upstream emission regions on a detection, and the spatial resolution of emission regions is not limited by the resolution of the underlying GCM. As such, high‐resolution maps of all upstream emission regions can be produced, which is critical for the search for emission sites within and around a comparatively small crater like Gale.

Methods

GCM Wind Simulations

Since global, high‐quality wind observations on Mars have been lacking to date, we use MarsWRF, a GCM of the Martian atmosphere, to simulate the wind fields necessary for inverse Lagrangian modeling. MarsWRF is derived from the terrestrial Weather Research and Forecasting (WRF) model and is a Mars‐specific implementation of PlanetWRF (Richardson et al., 2007; Toigo et al., 2012). MarsWRF is a finite‐difference grid‐point model projected onto an Arakawa‐C grid with user‐defined horizontal and vertical resolutions. The vertical grid follows a modified‐sigma (terrain‐following) coordinate from the surface to ∼80 km altitude. The total present‐day atmospheric CO2 budget is tuned to fit the Viking Lander annual pressure curves (∼6.1 mbar), and both surface albedo and thermal inertia are matched to Mars Global Surveyor Thermal Emission Spectrometer (MGS‐TES) observations (Christensen et al., 2001; Putzig et al., 2005), while a Mars Orbiter Laser Altimeter topography base map is employed and scaled to the chosen model resolution (Smith et al., 2001).

Multiple studies in the past have validated MarsWRF through comparison of its behavior against data from the Mars Global Surveyor Thermal Emission Spectrometer (Guzewich et al., 2013, 2014; Lee et al., 2011; Toigo et al., 2012), the Mars Reconnaissance Orbiter Mars Climate Sounder (Guzewich et al., 2013), and the weather stations on board Curiosity (Fonseca et al., 2018; Newman et al., 2017) and InSight (Newman et al., 2020), showing MarsWRF reproduces observed atmospheric pressure, atmospheric and ground temperature, near‐surface wind speeds and wind directions well.

MarsWRF permits multiple embedded “nests” with increasing spatial resolutions in a single model run. This allows atmospheric circulations influenced by small‐scale topographic features to be fully resolved in a simulation while the simulation also covers the entire globe at a coarser resolution. In this study, we run MarsWRF at increasing horizontal resolutions around Gale crater. The final model consists of four nested levels, each scaled up in resolution (spatial and temporal) by a factor of three from its “parent” nest. Level 1 provides global coverage with a horizontal resolution of 2° × 2° and a 60‐s timestep. Level 2 encompasses an 80° × 80° domain with a horizontal resolution of 0.67° × 0.67° and a 20‐s timestep. Level 3 encompasses a 26.67° × 26.67° domain with a horizontal resolution of 0.222° × 0.222° and a 6.67‐s timestep. Level 4 encompasses an 8.89° × 8.89° domain with a horizontal resolution of 0.074° × 0.074° (4.4 × 4.4 km) and a 2.22‐s timestep (Figure S1 in Supporting Information S1). Level 4 fully resolves the crater circulation. Two‐way boundary conditions link a nested domain with its parent, with information being passed both up and down between parent and child domains. A description of this process may be found in Richardson et al. (2007). In order to speed up the simulations, we performed test simulations in advance to determine the duration of MarsWRF simulations on each nesting level. A particular nest is no longer necessary after 99% of the initially released backward‐traveling particles have left the domain of that nest.

Given the lack of a global coverage of high‐quality wind observations, it is impossible to precisely reproduce “real” atmospheric circulations on spatial scales smaller than tens of kilometers, as stochastic weather events can significantly impact wind speed and even direction. As a result, we do not intend to reproduce the “real” winds. Rather, we aim to produce “mean” winds that are representative of their respective season and time of sol. For each TLS measurement, we repeat MarsWRF simulations for the corresponding Mars year four or five times, each time starting from a different initial condition. For each Mars year, the different rounds of GCM simulations are all driven by the same seasonally representative dust loadings. Results show slight variations in year‐to‐year conditions as a consequence of stochastic variability in the weather. The variance in results across the four or five times of simulation is, however, small.

On short timescales (<1 week), it is not anticipated there will be a significant change in the mean atmospheric conditions on Mars, so for each of the four or five rounds of wind simulation, we treat the sol of the methane measurement, and one, two, and three sols before and after the measurement as equally representative of the circulation pattern at the time of the TLS measurement, and release particles at the time of sol of each measurement on all of the seven sols. In this way, we form an ensemble of ∼30 back trajectory simulations for each investigated TLS measurement. Then, for each measurement, we average the results of its ∼30 back trajectory simulations. This ensures that discrete weather patterns are smoothed out.

Inverse Lagrangian Analysis

The wind fields from MarsWRF are used to drive the Stochastic Time‐Inverted Lagrangian Transport (STILT) Lagrangian Particle Dispersion Model (Fasoli et al., 2018; Lin et al., 2003) to simulate plume transport and dispersion. STILT is based on the Hybrid Single Particle Lagrangian Integrated Trajectory model (Draxler & Hess, 1998; Stein et al., 2015) that is extensively used in air quality, volcanic ash and industrial plume modeling, and STILT inherits all of the validated components of its predecessor. Although STILT was originally designed for terrestrial use, the Monin–Obukhov similarity theory for the PBL, along with the adherence to the well‐mixed criterion (Thomson, 1987), a manifestation of the second law of thermodynamics, ensures that the physics and fluid dynamics underlying STILT can be applied to all substantial planetary atmospheres, including the atmosphere of Mars, after modifications to some model parameters. The modified parameters include planetary radius, gas constant, angular rotation rate of the planet, surface gravity, dynamic viscosity of air, mean free path of air, molecular weight of air, surface air pressure, specific heat capacity of air, the map of land use, and the map of surface roughness length (Hébrard et al., 2012), etc.

In its application, STILT transports an ensemble of computational particles (10 thousand particles in each simulation in this study) from the site of the detector (here and henceforth, the location of the Curiosity rover) using time‐reversed grid‐scale wind plus a parameterized subgrid‐scale turbulent velocity (Hanna, 1984). The timestep in STILT is determined dynamically based on the wind field, and typically ranges between 1 and 10 min. When combined with a GCM, STILT linearly interpolates the simulated bulk wind from the GCM grid points to the precise positions of each particle at each timestep, and then displaces the particles according to the reversed wind arrow. Meanwhile, STILT adds a random velocity component, determined by a Markov chain process based statistically on the simulated meteorological conditions, to the bulk wind velocity. The random velocity represents turbulent motions that are unresolved by the GCM and results in dispersion of the particle cloud (Figure S3 in Supporting Information S1). Additionally, vertical mixing in the PBL is parameterized by vertically redistributing particles to random altitudes within the PBL (Figure S4 in Supporting Information S1). In the hyper‐near field around the detector, an “effective mixing depth” smaller than the PBL thickness is calculated based on the homogeneous turbulence theory, and particles are redistributed within the effective mixing depth rather than the entire air column within the PBL (Fasoli et al., 2018). This prevents the particles released near the surface from ascending to the top of the PBL instantaneously.

At every timestep in a back trajectory simulation (which corresponds to a putative emission time), STILT tallies the instantaneous particle density in the PBL at all locations and generates a “footprint” map in units of ppbv μmol−1 (Lin et al., 2003), which quantifies the contribution of unit methane emission from a putative emission site to the methane mole fraction at the detector. The footprint value is proportional to the column‐integrated particle number density within the PBL and the molar mass of air, and inversely proportional to the PBL thickness and the average air density within the PBL (Lin et al., 2003). The footprint value at any location (representing a putative emission site) at any time (representing a putative emission time) is equal to the prospective methane mole fraction in units of ppbv above the background level (∼0.41 ppbv) induced by 1 μmol of methane emission at that emission site and emission time. High footprint values indicate emission times and locations when and where emission casts strong influence over the detection. If integrated over all putative emission times, the footprint will measure the influence of a constant‐flux emission on a detection, and its map will show all important upstream regions.

In computing the footprints, the domain is first gridded horizontally (a grid that is separate from that of the GCM) so that STILT can count the number of particles within each horizontal grid and calculate the particle density at all putative emission locations. The resolution of this grid becomes the resolution of the map of the emission regions. We use 2° as the resolution for the domain from 80°S to 80°N and from 60°E eastward to 140°W. For the subdomain from 17.6°S to 8.4°N and from 124.2°E to 150.4°E, we use a resolution of 0.2°, or ∼11.8 km. For the subdomain from 6.64°S to 3.72°S and from 136.24°E to 139.16°E, we use a resolution of 0.02°, or ∼1.18 km. We note that the definition of the STILT footprint in this study is slightly different from its original form in Lin et al. (2003). The new definition has excluded the influence of the grid size and the timestep of footprint calculation on footprint values.

Results

Categorization of Methane Spikes

We focus on the seven methane spikes reported by the TLS instrument during the 7.1 years of the Curiosity mission through January 2020 (Figure 1, Table S1 in Supporting Information S1). The seven spikes can be categorized based on the season and the time of sol of their detections. In terms of seasons, Spikes 1 and 6 were detected in the late northern fall and winter. Spikes 2–5 and 7 were detected in northern spring. In terms of the time of sol, Spikes 1 and 5 were detected in the early afternoon, and Spikes 2–4, 6, and 7 were detected between midnight and early morning. As a result, Spikes 1 and 6 share similar seasonal, regional and global circulation patterns, as do Spikes 2–5 and 7. Spikes 1 and 5 share similar diurnal crater circulation patterns, as do Spikes 2–4, 6, and 7. The similarity in atmospheric circulation patterns also manifests itself in the subsequent STILT footprint maps.

Atmospheric Circulations

MarsWRF simulations show that the circulation at Gale crater consists of three components—a global meridional overturning circulation, a regional circulation, and a crater‐scale circulation. Figure 2 shows an example of near‐surface winds simulated by MarsWRF. In northern winter, the rising branch of the global meridional overturning circulation is centered in the southern hemisphere. Prevailing winds at the topographic dichotomy adjacent to Gale crater are southward and are particularly strong around 270° solar longitude (Ls = 270°), when Spike 6 was detected. In northern spring, the large‐scale prevailing winds at Gale crater are weak. The regional circulation is characterized by upslope northeasterlies along the topographic dichotomy in the afternoon, and downslope southwesterlies in the nighttime. The crater circulation is characterized by upslope winds along the inner crater rim and the slope of Mount Sharp in the afternoon, and downslope winds in the nighttime. The PBL thickness at Gale crater also undergoes a daily cycle between a nighttime minimum thickness of tens of meters, and a daytime maximum thickness of about three kilometers, similar to previous findings by Fonseca et al. (2018).

Figure 2

Simulated winds in the bottom layer of MarsWRF at Ls = 81.84° (Spike 5). The plotted data are an average over the six hours indicated by the time period on the upper left of each panel. Panels (a and b) show the regional circulation, from which one can identify southwesterly downslope winds along the topographic dichotomy from midnight to sunrise, and northeasterly upslope winds from noon to sunset. Panels (c and d) show the Gale crater circulation, from which one can identify downslope winds along the inner wall of the crater rim and along Mount Sharp from midnight to sunrise, and upslope winds from noon to sunset. The crater circulation is well resolved by MarsWRF. Red colors show rising air. Blue colors show sinking air. Contours show surface elevation. Red stars mark the position of Curiosity.

Influential Upstream Regions

Figure 3 shows the time‐integrated footprints of Spikes 1 and 2. Refer to Figure S5 in Supporting Information S1 for the footprints of all seven methane spikes. Within Gale crater, the strongest footprint of Spike 1 lies to the north of the TLS detector (Figure 3a), which is also the case for Spike 5 (Figure S5m in Supporting Information S1). This means that these two early afternoon measurements are both more sensitive to the emission from the north than the emission from other directions. The similarity in the footprints for Spikes 1 and 5 is consistent with the similarity in the early afternoon crater‐scale circulation patterns at the Curiosity site, in which northerlies dominate, although Spikes 1 and 5 were detected in different seasons. For Spike 2, the strongest footprint lies over the entire northwestern crater floor (Figure 3d). This is also the case for Spikes 3, 4, 6, and 7 (Figures S5g, S5j, S5p, S5s in Supporting Information S1), although there are some finer spatial patterns in the footprint map of Spike 6. These five spikes were all detected in the nighttime or in the early morning when the PBL was shallow. The released computational particles are confined within the PBL and only dispersed horizontally; therefore, they imprint almost equally strong footprints over the entire northwestern crater floor. In consequence, a nighttime detection is sensitive to the emission from all the northwestern crater floor.

Figure 3

Maps of time‐integrated STILT footprint, showing the influence of any putative emission site on (a–c) Spike 1 and (d–f) Spike 2. The maps are shown in (a and d) the crater scale, (b and e) the regional scale, and (c and f) the hemispherical scale. The footprint is integrated in putative emission time over a thirty‐sol time window prior to a methane measurement. High footprint values indicate upstream regions. The footprint value at any location (a putative emission site) is equal to the prospective TLS methane signal in ppbv above the ∼0.41 ppbv background after a thirty‐sol constant‐flux methane emission event with an emission flux of 1 μmol s−1 occurs at that location. The choice of integrating over 30 sols is based on the fact that due to dilution, few particles can still imprint footprints after traveling backwards for 30 sols; hence, the maps here show almost all the footprints that can be imprinted, and one can infer almost all the possible upstream regions from these footprint maps. Contours show surface elevation. Stars in (a and d) mark the positions of Curiosity.

Outside Gale crater, the strongest footprint of Spike 1 lies to the north of the crater, as a result of the prevailing regional‐scale northerlies in this season (Figure 3b). This is also true for Spike 6 (Figure S5q in Supporting Information S1). This means that, for these two spikes, if an emission region exists in the neighborhood of Gale crater (but outside the crater), it is most likely located to the north of the crater. The locations of the upstream regions of Spike 2 are, however, less definitive. They are located in the first and third quadrants of Gale crater (Figure 3e). This is also the case for Spikes 3–5, and 7 (Figures S5h, S5k, S5n, S5t in Supporting Information S1). Despite this ambiguity, the strongest footprints of all the seven spikes overlap in a region within 300 km to the north of Gale crater. It is noteworthy that the “E8” and “ESE” regions (142–146°E, 2–10°S), suggested as the most likely emission regions for Spike 1 by Giuranna et al. (2019), do not bear high footprint values in our study and are, hence, not identified as the preferred upstream regions for Spike 1 (Figure 3b).

Further zooming out to the hemispherical scale, the high‐footprint regions of Spike 1 extend from Elysium Planitia into two directions—one heading for the north along the western side of Elysium Mons to Utopia Planitia, and the other heading for the east along the southern side of Elysium Mons to Amazonis Planitia (Figure 3c). This is also the case for Spike 6, although the northern branch appears more prominent (Figure S5r in Supporting Information S1). This suggests that among all the distant large‐scale geographic units, the aforementioned ones are the most likely to be the emission regions for Spikes 1 and 6. For Spikes 2–5 and 7, the high‐footprint regions cover many large‐scale geographic units around Gale crater (Figure 3f, Figures S5i, S5l, S5o, S5u in Supporting Information S1), including the aforementioned Elysium Planitia and Utopia Planitia. It is noteworthy that many thousands of mud‐volcano‐like structures have been found in Utopia Planitia, which could be methane emission hot spots, according to Oehler and Etiope (2017).

Minimum Methane Emission

Based on these footprints, the minimum amount of methane emitted from any putative emission site that could give rise to the observed methane spikes can be calculated. The ∼0.41 ppbv background level from TLS is first subtracted from the seven methane spikes. The remainder of the signals must then be a consequence of recent emission. It is unknown whether the emission was continuous, intermittent, or episodic, but, to put a lower bound on the required methane emission, for each methane spike, we can assume that an instantaneous emission event occurred at the exact moment when a putative emission site had the strongest influence on the methane measurement. Then, for any location (representing a putative emission site), dividing each methane signal (after subtracting the background level) by the maximum footprint value of all putative emission times yields the minimum amount of methane emitted from that putative emission site that would be required to explain the methane signal (Figure 4). Upstream regions, which had the highest footprint values in Figure 3, now bear the smallest values in Figure 4, the latter meaning that they can more easily produce a methane spike. For example, any putative emission site on the northwestern crater floor (the blue region in Figure 4d) is able to produce Spike 2 by emitting only ∼100 kilograms of methane (refer to the left colorbars in Figure 4). In order to directly compare with the results from TGO, under the assumption that methane is a long‐lived species in the atmosphere, this ∼100 kilograms of methane will result in an increase of ∼10−5 ppbv in the global mean methane concentration (refer to the right colorbars in Figure 4). In comparison, if Spike 2 results from an emission event in Utopia Planitia, at least several millions of kilograms of methane must have been emitted, which would result in an increase of several hundreds of pptv (parts‐per‐trillion‐by‐volume, 1 pptv = 10−3 ppbv) in the global mean methane concentration (Figure 4f). For Spike 1, only a small, fan‐shaped area to the north of the Curiosity site is able to produce the spike by emitting ∼100 kilograms of methane (Figure 4a), as a consequence of vigorous mixing in the daytime PBL. Figure S6 in Supporting Information S1 shows the maps of the minimum methane emission for all the seven methane spikes.

Figure 4

The minimum amount of methane emitted from every putative emission location that can produce (a–c) Spike 1 and (d–f) Spike 2. For every putative emission site, an emission event is assumed to occur at the exact moment when the site has the strongest influence on a methane measurement. The left colorbars show the minimum mass of emitted methane as required by the magnitude of the spikes. The right colorbars show the increase in the globally averaged methane concentration after one of the aforementioned smallest emission event occurs. Contours show surface elevation. Stars in (a and d) mark the positions of Curiosity.

The 0.02 ppbv upper limit (Montmessin et al., 2021) of TGO, interpreted as the upper limit on the average methane concentration in the Martian atmosphere, if combined with the 330‐year lifetime estimated from standard photochemical models (e.g., Lefèvre & Forget, 2009), implies that, on average, no more than 6 × 10−5 ppbv of methane (or ∼530 kg of methane) is replenished every year. Then, during the 7.1 years of TLS operation, on average, no more than ∼4.3 × 10−4 ppbv (or ∼3,700 kilograms) could have been emitted into the atmosphere. Assuming the seven methane spikes result from seven emission events, on average, each of them could emit no more than ∼530 kilograms of methane; otherwise, the methane release would have resulted in a significant and potentially observable rise in the background methane concentration. Only the blue areas in Figure 4 qualify as areas able to produce a methane spike with the observed mole fraction but emitting no more than ∼530 kilograms of methane (refer to the left colorbars in Figure 4). More quantitative analysis shows that these “qualified areas” are only 1,560 km2 in total, about 8.4% the total area of Gale crater or 1.1 in 100,000 the total surface area of Mars (Figure S10 in Supporting Information S1). This means that without fast removal mechanisms significantly shortening the lifetime of methane, any methane emission site responsible for the TLS methane spikes must be located within this 1,560 km2 region around the Curiosity site inside Gale crater. If the emission site is located in the interior of this region, there is still some room for other emission sites to exist elsewhere on the planet. If it is located right on the edge of this region, no other emission sites on Mars are allowed to exist. If an emission site has equal chances to reside at any location on the surface of Mars, the probability of it being located within the 1,560 km2 around the site of Curiosity will be equal to the fraction of this region in the total surface area of Mars, only 1.1 in 100,000. Despite this low probability, this is the only way that the TLS spikes and the TGO non‐detections can be reconciled if one accepts the ∼330‐year lifetime of methane. The 1,560 km2 area is still an overestimate, as the assumed situation, where only seven methane emission events occurred during the 7.1 years and all of them were captured by the TLS measurements, is essentially impossible. The actual methane spike frequency at the Curiosity site may be much higher, which will put a much lower upper bound on the amount of methane emitted by a single emission event. In that case, the “qualified emission region” will be confined within even smaller areas that are very close to the location of the Curiosity rover, such as the deep blue areas on the northwestern crater floor in Figure 4. Even for the 1,560 km2 upper limit of the “qualified emission region”, this invokes a coincidence that Curiosity was sent to the immediate vicinity of a methane emission hotspot, almost an impossibility. Although Gale crater was carefully selected as the landing site for Curiosity based on its unique geological context (e.g., Grotzinger et al., 2015), and it has been proposed as a potential methane seepage site because of its location on the heavily faulted topographic dichotomy that is favorable for microseepage (Oehler & Etiope, 2017), it is still very unlikely that one methane emission site on Mars, or even the only emission site on Mars, resides so close to the site of the rover. Therefore, probably at least one of the three assumptions needs to be reevaluated: (a) TLS's methane spikes are real, (b) TGO's upper limits are real and they represent the upper limit of the methane abundance throughout the Martian atmosphere, and (c) the lifetime of methane is ∼330 years. We note that some concerns about contamination by terrestrial methane in the foreoptics chamber and other parts of the TLS instrument have been raised (e.g., Zahnle, 2015), which may have complicated the retrieved in situ methane concentrations, but the possibility of measurement bias was ruled out by the TLS team (Webster et al., 2018). If unknown, rapid methane removal mechanisms are at work, which significantly shorten the lifetime of methane, more methane can be emitted into the atmosphere per year without perturbing the long‐term background methane concentration, and the emission sites will have some freedom to be located at more distant places outside Gale crater, most likely in the upstream regions found in Section 3.3.

Consecutive Methane Measurements

More precise emission site identification is possible when we make use of consecutive methane measurements that report a large difference in methane abundances. At Ls ≈ 266° in Mars Year 33, two measurements were obtained within a few hours. The first measurement started at ∼01:30 local time and detected a 0.332 ppbv signal, close to the background level. Only a few hours later, the second measurement, at ∼06:30 local time, detected Spike 6 with 5.55 ppbv. One possible explanation for the rapid increase in the ambient methane concentration is that an emission event was initiated between the two measurements. Here, we focus on another possibility that the two measurements were both performed in the midst of an emission event, but a change in the wind direction between the two measurements induced the temporal variability of the methane signals. For any emission event, it would produce a methane signal only if the emission site was located in the upstream region of the detector at the time of emission. Figure 5 shows a comparison between the time‐integrated footprints for Spike 6 and for the background level. A significant difference can be found between the upstream regions within Gale crater (Figures 5a and 5d). On the northwestern crater floor, the upstream region of Spike 6, indicated by high footprint values, primarily lies to the west and the southwest of Curiosity rover, whereas the upstream region of the background level primarily lies to the northeast of the rover. Therefore, if one assumes the TLS measurements to be correct, the region to the west and the southwest of Curiosity in northwestern Gale crater is identified as the preferred methane emission site for this spike. There are no significant differences between the upstream regions at larger scales (Figures 5b and 5e, and Figures 5c and 5f).

Figure 5

Comparison between the footprint maps for a background level and for Spike 6, which were measured on the same sol at Ls ≈ 266° in Mars Year 33. Panels (a–c) show the STILT footprint of the background concentration measured at ∼01:30 local time. Panels (d–f) show the STILT footprint of Spike 6, measured at ∼06:30 local time. Contours show surface elevation. The stars in (a and d) show the positions of Curiosity. An emission site with weak influence on the background level and strong influence on Spike 6 would bear a small footprint of the former and a large footprint of the latter. Comparing (a) and (d), regions to the west and the southwest of Curiosity on the northwestern crater floor are such sites. The differences between (b) and (e) and between (c) and (f) at the larger scales are smaller than the uncertainty.

We note that this method based on consecutive methane measurements is able to precisely constrain the location of a nearby emission site, but it has high requirements for the measurements. First, the measurements must be consecutively performed within a short period of time, such that the multiple detections have a good chance to be influenced by the same emission event. We have no information about how long a methane emission event on Mars would typically last, but an initial guess for the time period will be at most a few days and optimally a few hours. Second, the measurements must be performed at different times of sol. As currently only “mean” winds representative of their seasons and time of sol can be simulated with a GCM, the winds at the same time of sol on a few consecutive sols are similar, making it difficult to tell the difference between the upstream regions of different measurements. For example, we attempted to apply this method to Spike 7 detected at 03:51 local time on Sol 2442 (after the landing of Curiosity) and a background level detected at 01:52 on Sol 2446, but, because these two measurements were both made in the middle of night, no significant difference is found between their upstream regions; hence, this method fails to precisely locate the emission region based on these two measurements. The drop in the measured methane concentration during the four sols is probably due to the cessation of an emission event.

Conclusions

If we trust the methane abundances detected by both TLS and TGO and accept the 330‐year methane lifetime from known photochemistry, our back trajectory modeling for atmospheric transport strongly supports surface emission sites in the vicinity of the Curiosity rover in northwestern Gale crater. This includes the special case that TLS, itself, is the methane source. However, an emission site in northwestern Gale crater invokes a coincidence that we selected a landing site for Curiosity so close to an active methane emission site, which is a small probability event. Other possibilities that do not invoke this coincidence include the existence of fast methane removal mechanisms that are unknown to date, false positives of TLS and/or false negatives of TGO. Should future studies confirm the existence of heterogeneous pathways or other unknown photochemical processes for methane destruction, the methane emission sites can be located outside Gale crater, and most likely to the north of the crater.

Our study demonstrates the feasibility and the advantages of applying the inverse Lagrangian modeling technique to emission site identification problems on other planets. The difference between our conclusions about the probable location(s) of the putative emission site(s) and the conclusions in Giuranna et al. (2019) may demonstrate the necessity of small‐scale wind simulations and repetitive simulations for removing weather stochasticity. Methane concentration data from future in situ measurements, especially those collected in consecutive measurements performed within a few hours, could further improve the emission site identification.

Supporting Information S1

Click here for additional data file.