A three-dimensional atmospheric dispersion model for Mars

Atmospheric local-to-regional dispersion models are widely used on Earth to predict and study the effects of chemical species emitted into the atmosphere and to contextualize sparse data acquired at particular locations and/or times. However, to date, no local-to-regional dispersion models for Mars have been developed; only mesoscale/microscale meteorological models have some dispersion and chemical capabilities, but they do not offer the versatility of a dedicated atmospheric dispersion model when studying the dispersion of chemical species in the atmosphere, as it is performed on Earth. Here, a new three-dimensional local-to-regional-scale Eulerian atmospheric dispersion model for Mars (DISVERMAR) that can simulate emissions to the Martian atmosphere from particular locations or regions including chemical loss and predefined deposition rates, is presented. The model can deal with topography and non-uniform grids. As a case study, the model is applied to the simulation of methane spikes as detected by NASA’s Mars Science Laboratory (MSL); this choice is made given the strong interest in and controversy regarding the detection and variability of this chemical species on Mars.


Introduction
Chemical species emitted into the atmosphere can travel hundreds or even thousands of kilometers from their source location. Their distributions can have a strong impact on the availability of species in the context of atmospheric chemistry. General circulation models (GCMs) on Mars, which typically have resolutions of the order of 1º of latitude and longitude (e.g., Richardson et al. 2007), usually include routines for computing the distribution of passive tracers. In some cases, such models have been extended to deal with chemistry (e.g., Lefèvre and Forget 2009); these models have proven to be an essential tool for studying the Martian atmosphere, particularly the atmospheric dynamics and chemistry, dust, and water cycles on Mars. However, the resolution constraints in GCMs limit their use in the local-to-regional domain. Mesoscale/microscale models, or GCMs with zoom capabilities, have successfully operated in this scale for weather forecasting on Earth.
They have only recently been applied to other planetary atmospheres (e.g., Rafkin et al. 2001;Tyler et al. 2002;Richardson et al. 2007;Spiga and Forget 2009), which has expanded our knowledge in several fields, such as atmospheric circulation, aeolian processes, the water cycle, and methane studies (Spiga and Forget 2009 and references therein;Steele et al. 2017;Pla-Garcia et al., 2019).
These models calculate the dispersion of species in the atmosphere using an online approach; that is, the governing equations of species are embedded in a meteorological model and are solved simultaneously with the meteorological predictions. This online modeling of species dispersion has the advantage of representing the atmosphere more realistically in cases where species distribution (e.g., aerosols) influences meteorology, but it results unnecessary and has some disadvantages in other cases. One such disadvantage arises when many simulations with different chemical assumptions are required; in such cases, longer computational times are needed to produce meteorological forecasts with online air quality predictions. Offline chemical transport simulations only require a single meteorological dataset to produce many chemical transport simulations for a particular study  (Grell et al. 2004;Leelossy et al. 2014), such as in analysis of emission scenarios; that is, it is not required to run the complete meteorological prediction with each dispersion simulation. When studying Mars, researchers usually have to deal with sparse observations or data that involve insufficient temporal and spatial resolution; therefore, the use of models to assist in the interpretation of observations and to test several scenarios is mandatory. While specific offline local-to-regional atmospheric dispersion models (hereafter referred to as dispersion models) are widely used on Earth to predict and study the effects of both local and regional emissions (Leelossy et al. 2014 and references therein), to date, such models have not been developed for Mars. This paper introduces the three-dimensional atmospheric dispersion versatile model for Mars (DISVERMAR), as a complementary tool to the modeling efforts performed to date in the Mars atmosphere. The new model is able to easily simulate both point and regional emissions into the Martian atmosphere. First-order chemical loss can also be included in the simulations, as well as predefined deposition rates for the chemical species of interest. The model takes meteorological data from existing and well-tested regional circulation models such as MarsWRF (Richardson et al. 2007;Toigo et al. 2012;Newman et al. 2017). The model can operate in different temporal and spatial domains from local to regional, being constrained by the resolution of the input meteorological data (wind, temperature, and pressure fields). The model can also deal with topography and non-uniform grid spacing, which is particularly useful for studying surface emissions and boundary layer processes, in which higher resolution is needed close to the surface than at higher levels.
This paper demonstrates how DISVERMAR could assist in the interpretation of data acquired by a surface mission on Mars. In particular, this work presents a case study involving the simulation of methane spikes detected by NASA's Mars Science Laboratory (MSL), which have been attributed to surface emissions in the vicinity of the MSL (Webster et al. 2015;Webster et al. 2018). Methane on Mars represents a puzzle in that there is not a satisfactory explanation of its origin, source locations, and detected variability (e.g., Yung et al. 2018;Moores et al. 2019;Viúdez-Moreiras et al. 2020a), including apparent mismatches among observations (e.g., Webster et al. 2015;Webster et al. 2018;Korablev et al. 2019), and it is therefore a potential application for DISVERMAR.

Methods
The developed three-dimensional numerical model solves the advection-diffusion equation in an Eulerian framework, which is suitable for the local-to-regional scale in which DISVERMAR will operate, including planetary boundary layer (PBL) studies and chemistry (Leelossy et al. 2014 and references therein): where N is the species number density, v the threedimensional wind field, D e the eddy diffusivity matrix, and S , R , and D the source, loss, and deposition terms, respectively. The discretization of this PDE is performed by high-order central differences for the diffusive term and an upwind scheme for the advective term (e.g., Durran 2010;Hirsch 2007). The order can be selected by the user, for each particular application. The implemented code in the model allows the option to integrate the PDE with the forward Euler scheme or by using a predictorcorrector scheme; the selection will depend on the application. The model makes use of a standard Arakawa-C grid ( Fig. 1 left), wherein the chemical species is staggered by one-half a vertical grid below the vertical velocity w, and the u and v wind components are at the same level as the chemical species but staggered horizontally by one-half of the horizontal grid. A Mercator projection is applied to the horizontal domain. In the vertical domain, the model can operate in z-coordinates or in terrain-following σ-z-coordinates. The last vertical coordinates system is used to deal with topography, where σ is defined according to the equation: where z srf = z top − z srf , z top is the model top domain coordinate and z srf the ground level. In that case, the gradient operator present in Eq. 1 is modified to take into account the topography in accordance with Hirsch (2007). Furthermore, the vertical grid can deal, if needed, with higher resolution close to the surface, decreasing progressively as altitude increases in order to refine the mesh near the surface, e.g., in case surface emissions are being modeled. The vertical grid spacing and decreasing rate is user specified.
The upwind scheme employed in DISVERMAR is a class of a non-oscillatory Lax-Wendroff scheme based on Smolarkiewicz and Margolin (1998), which corrects the truncation error in an iterative manner. This scheme can also correct the truncation error term proportional to flow divergence, and consider diffusion and source terms, ensuring a fully second-order scheme. Thus, for a basic presentation of the scheme and the advection in two dimensions, the first pass, or pass k = 0, of the algorithm computes the uncorrected species number density for the iteration n + 1, N where the flux function F is defined in terms of the local Courant number, C x and C y , for the x-and y-axes, respectively. Thus, for the x-axis, the flux function F is defined by: where where v x is the wind speed in the x-axis, t the timestep, and x the mesh resolution in the x-axis. The y-axis presents similar equations. Then, the next pass k of the algorithm computes iteratively a corrected field of the species number density for the iteration n + 1 by: and B (n+1) k y matrices can be computed by the following equation on the right edge of a cell: The higher-order approximation of the species number density for the iteration n + 1 is computed for M algorithm passes. However, the two-pass scheme already offers second-order accuracy (Smolarkiewicz and Margolin 1998). The extension of the scheme from two dimensions to three dimensions is straightforward.
The three-dimensional wind field as a function of local time is obtained from a meteorological model. Although using the same grid in both models would be recommended, this situation is not usual given the grid optimization for the particular processes being simulated in each model. Furthermore, the terrain-following coordinates could differ between DISVERMAR and the mesoscale model used to retrieve the meteorological fields. Thus, such fields are interpolated to the DISVERMAR grid prior to the simulation to ensure mass conservation, as is usual in dispersion models (e.g., Ishikawa 1993). This capability makes DISVERMAR versatile in terms of using inputs derived from mesoscale meteorological simulations. The current implementation includes an interface with the MarsWRF model (Richardson et al. 2007;Toigo et al. 2012), although other meteorological models could be used to input the meteorological fields to DISVERMAR.
The transport of chemical species in the lower atmosphere of Mars, as is the case in the Earth's atmosphere, is usually dominated by advection in the horizontal domain. In the vertical, however, the strongly unstable conditions present during daytime play a central role in the resulting transport of species, and it must be taken into account whether the grid resolution is large in the sense that it does not resolve such eddies.
The closure turbulence parameterizations implemented in Earth atmospheric models are based on numerous field campaigns and data analyses that have been refined over the years. However, only a few surface missions have landed on Mars to date, and consequently, few data are available, with most of them having been obtained at surface level. Therefore, it is not possible to develop validated closure turbulence schemes for Mars based on these data, and even several datasets recently acquired on the Martian surface are challenging to interpret due to  (2021), the approach of current Mars atmospheric models in their closure schemes is to use parameterizations typically adapted from those developed for Earth, with schemes and their parameters being modified if necessary to minimize mismatches with the few available empirical data. The subgrid-scale transport (caused by mechanical shear and buoyancy) is parameterized in DISVERMAR by the eddy diffusivity matrix, D e , assuming the coefficients for heat (e.g., Draxler and Hess 1998). Zero values are considered for the non-diagonal elements (e.g., Smagorinsky 1963;Louis 1979;Draxler and Hess 1998). However, while a classical local scheme using this parameterization in Eq. 1 could be valid for stable and weakly unstable conditions, it may lead to mismatches if the simulation domain covers the Martian PBL at daytime, when strong unstable conditions are at work, and the model grid does not fully resolve the turbulence. The closure turbulence parameterization implemented in DISVERMAR integrates a hybrid scheme (Hong and Pan 1996) that modifies the diffusive term ∇(D e * ∇N ) in Eq. 1 by adding a countergradient flux γ to the vertical domain ∂ ∂z D zz ∂N ∂z − γ , thus accounting for non-local effects that could be significant within the PBL, such as large eddies. The countergradient term is particularly relevant for properly computing the thermal structure of the PBL and thus the height of the PBL. In fact, most of the differences using this approach usually come from the non-local effects of temperature (Hong et al. 2006). The countergradient term γ is given by: where b is a coefficient taken as 7.8 (Hong and Pan 1996), S f is the species surface flux, and ω s is the mixed-layer velocity scale taken as ω s = u * /φ m , where u * is the friction velocity and φ m is a wind profile function dependent on the stability conditions, which also depends on the height of the PBL, h (Hong and Pan 1996). The vertical diffusivity coefficient, D zz , is computed by: where k is the Von Karman constant, Pr is the Prandtl number, and the coefficient p is taken to be 2 (Troen and Mahrt 1986;Hong and Pan 1996). Some parameters of the closure scheme, such as the height of planetary boundary layer h and the friction velocity u * , can be taken from the meteorological model outputs if available, as is usually the case in dispersion models for Earth. The Hong and Pan (1996) scheme has been widely validated on Earth and is also included in MarsWRF (MRF scheme) (see, e.g., Richardson et al. 2007). This scheme presents good performance against the few data available for Mars, although it should be highlighted that comparisons with MSL REMS surface data show lower predicted surface wind speeds probably due to a mismatch in the mixing of momentum down to the near-surface from higher layers (Newman et al. 2017).
The eddy diffusivity coefficients in the horizontal domain are computed from the wind gradient by: Fig. 1 The DISVERMAR horizontal and vertical computational grid (A, B panels, respectively) where h is the horizontal grid spacing and c is a factor to be taken as 2 × 10 −3 (Draxler and Hess 1998). The loss and deposition terms are parameterized in the current version of DISVERMAR by a characteristic time constant, ζ l and ζ d , respectively, leading to The process removal rate in a timestep, J , is obtained by solving the aforementioned relationship: where ζ is its characteristic time constant. Then, species number densities are updated accordingly. For a deposition process, ζ = ζ d in the first prognosed layer and ζ equals zero in the levels above.

Validation of the numerical model
DISVERMAR has been verified against scenarios with analytical solution. This subsection presents some selected test cases performed with DISVERMAR in order to verify the model implementation, based on the test cases presented in Sankaranarayanan et al. (1998).
Here, the PDE is discretized for the diffusive term by second-order central differences, and the integration in time is performed, for simplicity, with the forward Euler scheme. The vertical domain operates in z-coordinates using a uniform grid. The high-order upwind scheme considers two passes in these examples.
The first test (test 3.1.1) corresponds to the simulation of the one-dimensional advection and diffusion of a Gaussian pulse. The analytical solution is given by: where x is the spatial coordinate, t is the time, N (x, t) is the concentration of the species being dispersed as a function of location and time, x 0 is the coordinate of the center of the Gaussian pulse of unitary height at t = 0, v x is the wind velocity, and D x is the species diffusion coefficient.
The parameters for this test are based on Sankaranarayanan et al. (1998). Thus, D x = 0.005 m 2 /s, v x = 0.8 m/s, and the Gaussian pulse of unitary height is initialized at x 0 = 1 m. The space step, Δx, and timestep, Δt, are taken to be 0.025 m and 0.0125 s, respectively. The Courant number (v x Δt/Δx) is 0.4, and the Pèclet number (v x Δx/D x ) equals 4. Results are shown in Fig. 2. As can be seen, the numerical model predicts with excellent performance the dispersion of the Gaussian pulse, presenting very low discrepancies with the analytical model. The second test (test 3.1.2) enlarges test 3.1.1 to a twodimensional domain. The analytical solution is given in this case by: where x and y are the spatial coordinates, t is the time, N x, y, t is the concentration of the species being dispersed as a function of location and time, (x 0 , y 0 ) is the location of the center of the Gaussian pulse at t = 0, v x and v y are the wind components in the x-and y-axes, respectively, and D x and D y are the species diffusion coefficients in the x-and y-axes, respectively.
The parameters for this test are D x = D y = 0.01 m 2 /s, v x = v y = 0.8 m/s, and the Gaussian pulse of unitary height is initialized at (0.5, 0.5) m. Δx (= Δy) and Δt are taken to be 0.025 m and 0.00625 s, respectively, with a Courant number that is equal to 0.4. Results are shown in Fig. 3. As in the previous test, the numerical model performance is very good, again showing very low discrepancies with the analytical model. Thus, the analytical model predicts the center of the Gaussian pulse after 1.25 s to be located at (1.5, 1.5) with N(1.5,1.5) = 0.17, while the numerical model predicts the center of the pulse at the same location but with N(1.5,1.5) = 0.16.
The third test (test 3.1.3) considers the advection and diffusion of an instantaneous point source in a threedimensional domain. The analytical solution is given by: where where x, y, and z are the spatial coordinates, t is the time, N x, y, z, t is the concentration of the species being dispersed as a function of location and time, (x 0 , y 0 , z 0 ) is the location of the center of the point source at t = 0, v x , v y , and v z are the wind components in the x-, y-, and z-axes, respectively, and D x , D y , and D z are the species diffusion coefficients in the x-, y-, and z-axes, respectively.
The parameters for this test are D x = D y = 500 m 2 /s, D z = 0.01 m 2 /s, v x = v y = 0.2 m/s, v z = 0 m/s. The domain covers 40,000 × 40,000 m in both the x-and y-axes, and L z = 100 m in the z-axis. The spatial resolutions Δx (= Δy) and Δz are taken to be 2000 m and 5 m, respectively, and Δt = 100 s. The Courant number is equal to 0.02.
The instantaneous point source at unitary height is injected at (12,000, 12,000) m. For a smooth initial condition, the numerical simulation is initialized at t 0 = 5000 s and it is extended during 3 h, considering the species distribution derived from Eqs. 20-23 for t = t 0 as initial conditions. In addition, Dirichlet boundary conditions are taken from the analytical solution for each timestep, although the bulk of the initial mass injected in the domain at t 0 remains in the domain for the simulation time considered; thus, the concentration is close to zero at the boundaries.
Results are shown in Fig. 4 for two altitudes: (a) altitude z 1 , at L z /2, the center of the vertical domain (z 1 = 50 m) and (b) altitude z 2 , at L z /4 of the vertical domain (z 2 = 25 m). The left column corresponds to the initial time of the numerical simulation. (The instantaneous point source was injected at t = − t o in that reference system.) As can be seen, the species is dispersed through the domain, reaching z 2 , at lower altitudes with higher concentrations a couple of hours after the initialization, and decreasing again at those altitudes as the species is more dispersed. The model performance is good, again showing very low discrepancies with the analytical model. The maximum relative error between the analytical and numerical solutions is 7.12%, computed over the maximum concentration reached at each timestep. In absolute terms, the maximum deviation is always less than 4.0 × 10 −3 .

Case study application of DISVERMAR to a real problem on Mars
On Earth, methane is mainly of biogenic origin. However, understanding the origin of the methane detected by the MSL has proven to be a complex issue. This is exacerbated by the limited data that have been acquired from the Martian surface. Thus, the complexity of the methane measurements and the fact that the MSL has to distribute the electrical energy generated by its radioisotope thermal generator to various instruments, results in the limited number of methane measurements that have been taken to date. In situ observations by the MSL have suggested the detection of methane spikes (i.e., variable methane abundance in a short timescale). Webster et al. concerning the methane origin or its source location relative to the MSL. As MSL is a rover, if the methane source were relatively close to the rover, the MSL team could consider moving the MSL to the hypothetical region where methane is being emitted and investigating about its origin. The lack of extensive modeling and simulation supporting the sparse methane observations (e.g., a systematic dispersion analysis of methane from hypothetical source locations) and the lag between observations and computations, prevents to move forward such a possibility. The sparse data available from MSL measurements do not allow ascertaining the exact time scale of such events detected by the MSL, but the few data available suggest that they could occur on a timescale from a few minutes to less than a few sols (Webster et al. 2018;Moores et al. 2019;Giuranna et al. 2019).
Given that the origin of the methane observed on Gale Crater is currently unknown, the emission flux into the atmosphere could be either episodic or continuous. While the potential production of methane spikes by an episodic release of methane in the vicinity of the MSL should be expected, favoring these kinds of emissions as an explanation for the spikes (e.g., Lefèvre 2019), the production of methane spikes by a continuous release is an alternative that may deserve attention. Methane spikes produced by a continuous release cannot be explained plainly, and their occurrence may mostly depend on the location of the emission source, the suitability of the local meteorology in the region around MSL, and the variability of the emission flux (Viúdez-Moreiras et al. 2020a). The MSL is located inside Gale Crater, a region that, according to MSL's observations, is characterized by both complex topography (Fig. 5) and wind patterns (e.g., Newman et al. 2017;Viúdez-Moreiras et al. 2019). In order to demonstrate how DISVERMAR could assist in the interpretation of data acquired during a surface mission on Mars, the goal of this case study is to determine whether a continuous release of methane in the vicinity of MSL in Gale Crater could produce methane spikes. Being more specific, the objective is to determine whether the local meteorology around the MSL would aid in the production of methane spikes at MSL's location over a time scale of less than a sol, under the constraints imposed by MSL, i.e., sudden increases of as much as ~ 5-10 ppbv from base values less than 1 ppbv (Webster et al. 2015(Webster et al. , 2018. To accomplish this goal, a grid of possible source locations in a spatial domain covering a subset of Gale Crater (Fig. 5) around the MSL and centered on that MSL's general location, was established. Figure 6 shows the spatial domain of the simulations, including the distribution and emission area of the simulated source locations. MarsWRF mesoscale simulations were performed using the same setup as presented in Newman et al. (2017) for the Gale Crater region and during an illustrative season close to the southern winter to obtain the meteorological data. A single simulated sol after the MarsWRF simulation spin-up was used to input the meteorological fields to DISVERMAR, which was sufficient for the purposes of this study. The boundary-layer height and the friction velocity were obtained from MarsWRF model outputs (see Sect. 2). As in the previous tests, the diffusive term of the PDE is discretized using second-order central differences, and the high-order upwind scheme considers two passes. The integration in time is performed, for simplicity, with the forward Euler scheme. The grid is comprised of 50 × 50 horizontal nodes, with a spatial resolution of ~ 900 m, and 25 vertical levels. The vertical domain operates in σ-coordinates using variable grid spacing. The first vertical level is prognosed at ~ 10 m from the surface, with the resolution decreasing progressively as far as up the model top at 20 km. Dirichlet boundary conditions were used for the spatial domain, except on the surface, where a Neumann boundary condition is at work representing an emission flux of methane in a particular location. The methane is emitted in a region of ~ 7 km 2 (which covers a quadrant of nine grid elements which are represented by the gray squares shown in Fig. 6 (totaling 49 simulations). The worst-case scenario regarding the potential production of spikes is considered, that is, a continuous emission flux of methane without any variability in emissions (i.e., a constant emission flux) is included in the simulations. The emission flux was adjusted to match the methane mixing ratio observed by the MSL between 0 and 2 LTST at MSL's general location (Webster et al. 2018).
The theoretical gas-phase methane lifetime in the Martian atmosphere (Lefèvre and Forget, 2009) is included in the R term (Eq. 1), although negligible in this particular simulation given the time period coverage of a sol and the large methane lifetime expected in the Martian atmosphere. No deposition is included for these simulations.
A simulation of a continuous emission flux should be performed with care given that a spin-up period is also needed in dispersion simulations in order to avoid initialization issues in the distribution of methane in the domain. Simulations were initialized without methane and spun-up in order to allow a well-established diurnal cycle for the methane abundance in the domain. The spin-up period for these simulations was less than a sol (Additional file 1: Fig. S1). A timestep of 30 s was used, with the wind field being updated every 5 min. The results presented in this section were obtained from the simulation period after the spin-up of the DISVERMAR simulations.
An initial analysis of the simulation results showed two different groups of source locations. The first group, corresponding to locations to north and northwest of the MSL (15-17, 22-24 and 29-49, see Fig. 6), produced persistent methane abundances of more than 10-20 ppbv at MSL's location during the daytime. The methane abundance even exceeded 100 ppbv in several cases. The estimated emission fluxes were of the order of 10 3 kg/ sol (Additional file 1: Fig. S2). Therefore, these emission source locations could be incompatible with MSL's observations (Webster et al. 2015(Webster et al. , 2018. The behavior of the first group of emission sources at MSL's location is largely the result of a change in the wind regime over the slopes of Aeolis Mons (the mountain at the center of Gale Crater) and the crater floor during the diurnal cycle. The MSL is mainly subjected to southeasterly katabatic winds blowing over Aeolis Mons, which are reinforced by the regional and global circulation present at that season (e.g., Tyler and Barnes 2013;Rafkin et al. 2016;Newman et al. 2017;Viúdez-Moreiras et al. 2019), during the timeslot in which the MSL is measuring the background methane. This promotes air masses coming from southern locations but makes it difficult for air masses to come from the north. Given that the emission source must match MSL's observations, the ~ northern emission sources must have higher emission fluxes than the southern ones. During daytime, however, when the wind regime changes and upslope anabatic winds dominate this region on the slopes of Aeolis Mons, the northern emission sources may produce high and persistent methane abundance at MSL's general location. Figure 7 illustrates this mechanism. This figure depicts the surface methane distribution produced by two emission sources at different locations (locations 18 and 47) in the vicinity of the MSL, during both nighttime (0:00 h LTST) and daytime (12:00 h LTST). Both emission sources present similar methane abundances, in accordance with MSL's observations during nighttime; however, the unfavorable air mass transport during nighttime for the northwestern location (location 47) requires a much higher emission flux to match the MSL's observations. This implies   (Webster et al. 2015;. Figure S3 shows the simulated methane abundance during the diurnal cycle at MSL's general location in the first prognosed layer, for this group of emission locations, whose daytime abundance could reach even higher abundances, therefore incompatible with MSL. In any case, methane spikes were also observed in this group of simulations. As can be seen, each emission source produces a different pattern in the methane levels as a function of the local meteorology. Thus, northeastern locations will tend to produce present higher abundances during the afternoon, while eastern locations tend to result in higher abundances during the morning. Conversely, the second group of simulations (i.e., the locations to the south and southeast of MSL, 1-14, 18-21, and 25-28, see Fig. 6) produced suitable levels at MSL's general location during daytime and, in almost all cases (73%), methane spikes at MSL's location, even without taking into account other mechanisms that are not resolved in the model simulations, such as microscale flows, which could significantly increase the variability of methane abundance as detected by the MSL. Figure 8 shows the simulated methane spikes along the diurnal cycle in the first prognosed vertical level of the model (~ 10 m) at MSL's general location, for the emission locations corresponding to this group. As can be seen in the figure, source locations in the vicinity of the MSL can produce methane spikes of to 25 ppbv, both during the nighttime and daytime, above base levels of less than 1 ppbv between 0-2 LTST (in accordance with the background observations of MSL). Forty-five percentage of the emission locations produced methane spikes between 5 and 30 ppbv, and 27% between 1 and 5 ppbv. The remaining 27% did not produce significant spikes. The emission fluxes were highly variable for this group, ranging from 10 −1 to 10 3 kg/sol depending on the emission location (Additional file 1: Fig. S2). Western locations (Fig. 8, blue color) presented the highest emission fluxes due to the unfavorable winds in most of sol (during both nighttime and daytime), whereas the eastern locations within the southern-southeastern group of emission sources (Fig. 8, red color) presented the lowest due to the favorable ~ southeasterly flows descending on the slopes of Aeolis Mons during nighttime (Fig. 7). The timing and intensity of these spikes will also vary depending on many factors, such as the season and meteorological conditions. Figure 9 shows the evolution of the near-surface methane abundance in the region around the MSL considering an example of a methane spike produced during nighttime by a southern emission source (location 4). Figure 10 shows the volume mixing ratio, wind speed, and wind direction in the first prognosed layer at MSL's location. Given the proximity between the emission source and the rover, such winds can be considered as a proxy of the wind magnitudes along the trajectory of air masses between both locations. As shown in Fig. 9, near-surface east-southeasterly winds dominate around the source location at 0.5 h-2.5 h LTST, as a result of the katabatic downslope winds at Aeolis Mons. MSL's general location is not in the trajectory of the main air masses coming from the methane source, and therefore, the methane abundance is low at the rover's location. However, winds rotate to the south in a double-peak structure between 3 h and 4.5 h LTST (Figs. 9 and 10), transporting significant air masses to MSL's location, which causes a prompt increase in the methane abundance ( Figs. 9 and 10). In addition, this region experiences a calm period during which the dispersal of methane is mostly suppressed, thus favoring the spikes, particularly the ~ 10 ppbv methane spike. Downslope winds rotate to ~ south-southeasterly and strengthen again after 4.5 h LTST (Figs. 9 and 10) significantly lowering the methane abundance; wind direction remains roughly constant in the vicinity of the source location until 6.5 h LTST (Fig. 9). At 7.5 h LTST, winds rotate to ~ northeasterly, resulting in advective flows in the opposite direction to MSL's location and thus lowering methane abundance to negligible levels ( Fig. 9). Fulfilling the goal of this case study (to evaluate whether a continuous release of methane in the vicinity of MSL in Gale Crater could produce methane spikes), the DISVER-MAR simulations strongly suggest that the variability of winds in the vicinity of the MSL is sufficient to produce strong methane spikes at MSL's general location from a nearby emission source, even without invoking sol-tosol variability in winds, variable emission fluxes or eventual emissions into the atmosphere, or microscale flows not resolved in the meteorological data. The invocation of variable emission fluxes, eventual emissions, and/ or microscale flows could also dramatically increase the observed variability in the methane abundance at MSL's landing site. Therefore, methane spikes observed by the MSL could be the result of a nearby emission source and could develop in the worst-case scenario of a constant emission source, as the DISVERMAR simulations suggest. Further analysis should be performed to fully test this hypothesis, however.
In addition, the DISVERMAR simulations suggest that, in the absence of a more complex emission pattern, north-northwest source locations are incompatible with MSL's observations, favoring a hypothetical emission source at southern and southeastern locations (i.e. in the northwest slopes of Aeolis Mons).
As discussed above, the methane spikes observed at MSL's location would be the result of air masses with high methane content reaching the MSL before being dispersed in the atmosphere. Other locations in the nearsurface or even at other altitudes could reach higher abundances than those observed at MSL's location, which will depend on the complex three-dimensional advection and diffusion transport that occur in the atmosphere. Therefore, the methane abundance during spike events detected at MSL's location should not be considered as indicative of an abundance of methane in the near-surface atmosphere and beyond. The maximum levels of methane in the atmosphere will strongly depend on the location of the emission source. Thus, for example, methane levels reached more than ~ 100 ppbv in locations near the emission source, for the spike presented previously (Fig. 9), which would be much higher than the maximum Fig. 10 Volume mixing ratio, horizontal wind speed (vh), and wind direction (Ψ d ) in the first prognosed layer at MSL's location, covering the methane spike presented in Fig. 8 and produced by a source location at 4.7 deg of latitude, 137.4 deg of longitude. Wind direction is shown from which the wind is blowing, following the standard meteorological convention (0° < Ψ d ≤ 360°, with the origin pointing to the north)