(2021) 8:53
Viúdez‑Moreiras ﻿Prog Earth Planet Sci
https://doi.org/10.1186/s40645-021-00445-4

RESEARCH ARTICLE

Progress in Earth andPlanetary Science

Open AccessA three‑dimensional atmospheric dispersion
model for Mars
D. Viúdez‑Moreiras*

Abstract: Atmospheric local-to-regional dispersion models are widely used on Earth to predict and study theeffects of chemical species emitted into the atmosphere and to contextualize sparse data acquired at particularlocations 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 Eule‑rian atmospheric dispersion model for Mars (DISVERMAR) that can simulate emissions to the Martian atmospherefrom 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 inter‑
est in and controversy regarding the detection and variability of this chemical species on Mars.
Keywords: Martian atmosphere, Atmospheric dispersion modeling, Atmospheric chemistry
1 Introduction
Chemical species emitted into the atmosphere can travel
hundreds or even thousands of kilometers from theirsource 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 andchemistry, dust, and water cycles on Mars. However, theresolution constraints in GCMs limit their use in thelocal-to-regional domain. Mesoscale/microscale models, or GCMs with zoom capabilities, have successfully
operated in this scale for weather forecasting on Earth.
*Correspondence: viudezmd@inta.es
Centro de Astrobiología (CSIC-INTA) & National Institute for Aerospace
Technology (INTA), Torrejón de Ardoz, Madrid, Spain

They have only recently been applied to other planetaryatmospheres (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, andmethane 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, thegoverning equations of species are embedded in a meteorological model and are solved simultaneously withthe meteorological predictions. This online modeling ofspecies dispersion has the advantage of representing theatmosphere more realistically in cases where species distribution (e.g., aerosols) influences meteorology, but it
results unnecessary and has some disadvantages in othercases. 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© The Author(s) 2021. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://​creat​iveco​mmons.​org/​licen​ses/​by/4.​0/.

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

(Grell et al. 2004; Leelossy et al. 2014), such as in analysis
of emission scenarios; that is, it is not required to run thecomplete 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 useof models to assist in the interpretation of observations
and to test several scenarios is mandatory. While specificoffline local-to-regional atmospheric dispersion models(hereafter referred to as dispersion models) are widelyused on Earth to predict and study the effects of both
local and regional emissions (Leelossy et al. 2014 andreferences therein), to date, such models have not beendeveloped 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 toeasily simulate both point and regional emissions into theMartian 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 modeltakes 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 neededclose 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 casestudy involving the simulation of methane spikes detectedby NASA’s Mars Science Laboratory (MSL), which havebeen 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, anddetected 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.

2 Methods
The developed three-dimensional numerical model
solves the advection–diffusion equation in an Eulerianframework, which is suitable for the local-to-regional

Page 2 of 16

scale in which DISVERMAR will operate, including
planetary boundary layer (PBL) studies and chemistry(Leelossy et al. 2014 and references therein):∂N− ∇(N v) + ∇(De ∗ ∇N ) = S + R + D,∂t(1)

where N is the species number density, v the threedimensional wind field, De the eddy diffusivity matrix,
and S , R , and D the source, loss, and deposition terms,
respectively. The discretization of this PDE is performedby 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 theuser, 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 predictor–
corrector 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 samelevel as the chemical species but staggered horizontallyby one-half of the horizontal grid. A Mercator projection is applied to the horizontal domain. In the verticaldomain, the model can operate in z-coordinates or interrain-following σ–z-coordinates. The last vertical coordinates system is used to deal with topography, where σ isdefined according to the equation:

σ =ztop − z
,
�zsrf

(2)

where zsrf = ztop − zsrf , ztop is the model top domain
coordinate and zsrf the ground level. In that case, thegradient 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 themesh near the surface, e.g., in case surface emissions arebeing modeled. The vertical grid spacing and decreasingrate is user specified.
The upwind scheme employed in DISVERMAR is aclass of a non-oscillatory Lax–Wendroff scheme basedon 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 andsource terms, ensuring a fully second-order scheme.
Thus, for a basic presentation of the scheme and theadvection in two dimensions, the first pass, or passk = 0, of the algorithm computes the uncorrected spe(n+1)
cies number density for the iteration n + 1, Ni,j 1, by:

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

Page 3 of 16


 
n
n
n
F Ni,j , Ni+1,j , Cx,i+ 1 ,j
2


n
n
−F Ni−1,j, Ni,jn , Cx,i−
1
2 ,j

 
n
n
, Cy,i,j+− F Ni,jn , Ni,j+1
1
2


n
n
, Ni,jn , Cy,i,j−,−F Ni,j−1
1

(n+1)

2

(n+1)
(n+1)
Ni+1,j k − Ni,j k(n+1)k
,
1 =
(n+1)
(n+1)
x,i+ 2 ,jNi+1,j k + Ni,j kA(3)

where the flux function F is defined in terms of the local
Courant number, Cx and Cy, for the x- and y-axes, respectively. Thus, for the x-axis, the flux function F is definedby:

F (NL , NR , Cx ) = CxA NL + CxB NR ,

(4)

where

Cx =

vx t
,
x

(5)

CxA =

1
(Cx + |Cx |),
2

(6)

CxB =

1
(Cx − |Cx |),2(7)where vx 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: (n+1)k+1(n+1)k(n+1)k(n+1)k(n+1)kNi,j= Ni,j− F Ni,j, Ni+1,j , C
x,i+ 21 ,j


(n+1)
(n+1)
(n+1)k−F Ni−1,j k , Ni,j k , Cx,i− 21 ,j (n+1)k(n+1)k(n+1)k, Ni,j+1 , C− F Ni,j
y,i,j+ 21(n+1)k
(n+1)
(n+1)
,
−F Ni,j−1 k , Ni,j k , C
1
y,i,j− 2

(8)

where
(n+1)kCx(n+1)k−1
 (n+1)k−1 (n+1)k= Cx
Ax
 1 − Cx

(9)




(n+1)k−1
 (n+1)k−1 (n+1)k= CyBx
 1 − Cy(10)(n+1)k−1− Cx(n+1)kCy(n+1)k−1− Cx(n+1)k−1 (n+1)k
By
,

Cy

(n+1)k−1

Cy(n+1)k

Ax

.

(n+1)

k
k
The Axand Bymatrices can be computed by thefollowing equation on the right edge of a cell:(n+1)
Ni,j 1 = Ni,jn −(11)(n+1)k(n+1)k(n+1)k(n+1)k
1 Ni+1,j+1 − Ni,j+1 − Ni+1,j−1 − Ni,j−1(n+1)k.
B 1 =y,i+ 2 ,j
2 N (n+1)k − N (n+1)k + N (n+1)k + N (n+1)k
i+1,j+1i,j+1
i+1,j−1i,j−1(12)
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. Althoughusing the same grid in both models would be recommended, this situation is not usual given the grid optimization for the particular processes being simulated ineach 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 meteorologicalsimulations. The current implementation includes aninterface with the MarsWRF model (Richardson et al.
2007; Toigo et al. 2012), although other meteorologicalmodels 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 unstableconditions 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 thesense 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 refinedover 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 obtainedat surface level. Therefore, it is not possible to develop
validated closure turbulence schemes for Mars based onthese data, and even several datasets recently acquired on
the Martian surface are challenging to interpret due toViúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

Page 4 of 16

Fig. 1 The DISVERMAR horizontal and vertical computational grid (A, B panels, respectively)

instrumentation issues (see, e.g., Viúdez-Moreiras et al.
2019, 2020b). To develop fully validated Martian models, it will be necessary to wait for more sophisticatedmissions to be sent to Mars. Thus, as noted in Newman
et al. (2021), the approach of current Mars atmosphericmodels 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 mechanicalshear and buoyancy) is parameterized in DISVERMAR
by the eddy diffusivity matrix, De , 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 thisparameterization in Eq. 1 could be valid for stable andweakly 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∇(De ∗ ∇N ) in Eq. 1 by adding
flux γ  a countergradient∂N∂to the vertical domain ∂z Dzz ∂z − γ , thus account-ing 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:

γ =b

Sf
,
ωs

(13)

where b is a coefficient taken as 7.8 (Hong and Pan 1996),
Sf 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 theheight of the PBL, h (Hong and Pan 1996). The vertical
diffusivity coefficient, Dzz , is computed by:
1 z p
Dzz = kωs z 1 −
,
(14)
Pr
h
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 betaken 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 forMars, although it should be highlighted that comparisonswith 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:





∂u ∂v 21∂u ∂v 2
2
√
+
−
Dh =
+
,
c�h
∂y
∂x
∂x
∂y
2
(15)

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

Page 5 of 16

Fig. 2 Analytical (red color) and numerical (blue color) results for the test case 3.1.1

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 thecurrent version of DISVERMAR by a characteristic time
constant,
leading to  ζl and ζd , respectively,
N
N
∂N
=
−
=
−
R = ∂N
D
=
and
.
The
process
∂t
ζ
∂t
ζ
l

l

d

d

removal rate in a timestep, J , is obtained by solving the
aforementioned relationship:J = N 1 − e−�t/ζ ,(16)
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.

3 Results and discussion
3.1 Validation of the numerical model

DISVERMAR has been verified against scenarios withanalytical solution. This subsection presents someselected test cases performed with DISVERMAR in orderto verify the model implementation, based on the test
cases presented in Sankaranarayanan et al. (1998).
Here, the PDE is discretized for the diffusive term bysecond-order central differences, and the integration intime 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:(x−x0 −vx t )2
1
N (x, t) = √e− Dx (4t+1) ,4t + 1(17)

where x is the spatial coordinate, t is the time, N (x, t)
is the concentration of the species being dispersed as afunction of location and time, x0 is the coordinate of thecenter of the Gaussian pulse of unitary height at t = 0,vx is the wind velocity, and Dx is the species diffusion
coefficient.
The parameters for this test are based on Sankaranarayanan et al. (1998). Thus, Dx = 0.005 ­m2/s, vx = 0.8 m/s,
and the Gaussian pulse of unitary height is initialized
at x0 = 1 m. The space step, Δx, and timestep, Δt, aretaken to be 0.025 m and 0.0125 s, respectively. The Courant number (vxΔt/Δx) is 0.4, and the Pèclet number(vxΔx/Dx) equals 4. Results are shown in Fig. 2. As can beseen, the numerical model predicts with excellent performance the dispersion of the Gaussian pulse, presenting very low discrepancies with the analytical model.

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

Fig. 3 Analytical (A) and numerical (B) results for the test case 3.1.2, for three timesteps (t = 0.01 s, t = 0.25 s, and t = 1.25 s)

Page 6 of 16

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

Page 7 of 16

Fig. 4 Analytical (A) and numerical (B) results for the test case 3.1.3. Two altitudes are considered: (a) altitude z1, at Lz/2, the center of the verticaldomain (z1 = 50 m) and (b) altitude z2, at Lz/4 of the vertical domain (z2 = 25 m). Three timesteps are shown in columns (t = 0 h, t = 2 h and t = 4 h).
Note that the time domain shown in this figure corresponds with the numerical simulation, i.e., the instantaneous point source is injected at
t =  − to (see text)

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

Page 8 of 16

The second test (test 3.1.2) enlarges test 3.1.1 to a twodimensional domain. The analytical solution is given inthis case by:
N x, y, t =−1e4t + 1(x−x0 −vx t )2
Dx (4t+1)(y−y0 −vy t )2+
Dy (4t+1)

, (18)

wherex 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, (x0, y0) is thelocation of the center of the Gaussian pulse at t = 0, vx
and vy are the wind components in the x- and y-axes,
respectively, and Dx and Dy are the species diffusion coefficients in the x- and y-axes, respectively.
The parameters for this test are Dx = Dy = 0.01 ­m2/s,
vx = vy = 0.8 m/s, and the Gaussian pulse of unitary heightis 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 thenumerical model predicts the center of the pulse at thesame location but with N(1.5,1.5) = 0.16.
The third test (test 3.1.3) considers the advection anddiffusion of an instantaneous point source in a threedimensional domain. The analytical solution is given by:N x, y, z, t = Nx (x, ux , t)Ny y, uy , t Nz (z, uz , t),
(19)
where
−
1
Nx (x, ux , t) = √e
4π Dx t



4Dx t



(y−y0 −vy t )2





4Dy t

(z−z0 −vz t )2
4Dz t

(20)

,





−


1
e
Ny y, uy , t = 
4π Dy t−1
Nz (z, uz , t) = √e
4π Dz t(x−x0 −vx t )2

,

,

(21)

(22)

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, (x0, y0, z0) is thelocation of the center of the point source at t = 0, vx, vy,
and vz are the wind components in the x-, y-, and z-axes,
respectively, and Dx, Dy, and Dz are the species diffusion
coefficients in the x-, y-, and z-axes, respectively.

The parameters for this test are Dx = Dy = 500 ­m2/s,
Dz = 0.01 ­m2/s, vx = vy = 0.2 m/s, vz = 0 m/s. The domain
covers 40,000 × 40,000 m in both the x- and y-axes, andLz = 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 isinjected at (12,000, 12,000) m. For a smooth initial condition, the numerical simulation is initialized at t0 = 5000 s
and it is extended during 3 h, considering the species
distribution derived from Eqs. 20–23 for t = t0 as initialconditions. In addition, Dirichlet boundary conditions
are taken from the analytical solution for each timestep,
although the bulk of the initial mass injected in thedomain at t0 remains in the domain for the simulation
time considered; thus, the concentration is close to zeroat the boundaries.
Results are shown in Fig. 4 for two altitudes: (a) altitude
z1, at Lz/2, the center of the vertical domain (z1 = 50 m) and(b) altitude z2, at Lz/4 of the vertical domain (z2 = 25 m).
The left column corresponds to the initial time of thenumerical simulation. (The instantaneous point source
was injected at t =  − to in that reference system.) As can beseen, the species is dispersed through the domain, reaching z2, 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 ateach timestep. In absolute terms, the maximum deviationis always less than 4.0 × ­10−3.
3.2 Case study application of DISVERMAR to a real
problem on MarsOn 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 thelimited number of methane measurements that havebeen taken to date. In situ observations by the MSL havesuggested the detection of methane spikes (i.e., variablemethane abundance in a short timescale). Webster et al.
(2015) stated that the methane spikes could be the resultof a nearby source location. Previous modeling efforts
with a mesoscale model points in this direction (e.g., PlaGarcía et al. 2019). However, there is no solid evidence

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

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 hypotheticalsource locations) and the lag between observations andcomputations, 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 suggestthat 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 theatmosphere could be either episodic or continuous.
While the potential production of methane spikes by anepisodic 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 analternative that may deserve attention. Methane spikes
produced by a continuous release cannot be explained
plainly, and their occurrence may mostly depend on thelocation of the emission source, the suitability of thelocal meteorology in the region around MSL, and thevariability of the emission flux (Viúdez-Moreiras et al.

Page 9 of 16

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 thevicinity 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 asmuch as ~ 5–10 ppbv from base values less than 1 ppbv
(Webster et al. 2015, 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 thespatial 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 purposesof this study. The boundary-layer height and the friction

Fig. 5 A Topographic region around the Mars Science Laboratory (Gale Crater). Ground altitude is shown in kilometers above the Mars geoid. B, C
Topography along the latitude and longitude axes, respectively, passing along the MSL. MSL’s general location is also indicated by a dot

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

Page 10 of 16

Fig. 6 Spatial domain of the simulations, including the set of simulated source locations. The emission area is represented in gray squares for eachsimulated location. Methane emission flux is uniformly distributed in the emission area, which covers nine elements of the grid. Isolines depict the
ground altitude in kilometers. The MSL’s general location is highlighted by a red dot. Each emission location is labeled by a number. Locations within
φ ± 0.05º, where φ is the longitude of the MSL (~ 137.4°), are enveloped in фR2. Western locations to φ-0.05º are enveloped in фR1 and eastern ones
to φ + 0.05° in фR3velocity 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, andthe high-order upwind scheme considers two passes. The
integration in time is performed, for simplicity, with theforward 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, withthe resolution decreasing progressively as far as up themodel top at 20 km. Dirichlet boundary conditions wereused for the spatial domain, except on the surface, wherea 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 ­km2 (which covers aquadrant 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 particularsimulation given the time period coverage of a sol and thelarge 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 alsoneeded in dispersion simulations in order to avoid initialization issues in the distribution of methane in thedomain. Simulations were initialized without methaneand 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 thana sol (Additional file 1: Fig. S1). A timestep of 30 s wasused, with the wind field being updated every 5 min. The
results presented in this section were obtained from thesimulation 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

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

estimated emission fluxes were of the order of ­103 kg/
sol (Additional file 1: Fig. S2). Therefore, these emission source locations could be incompatible with MSL’s
observations (Webster et al. 2015, 2018). The behaviorof the first group of emission sources at MSL’s locationis largely the result of a change in the wind regime overthe slopes of Aeolis Mons (the mountain at the centerof Gale Crater) and the crater floor during the diurnal cycle. The MSL is mainly subjected to southeasterly
katabatic winds blowing over Aeolis Mons, which arereinforced 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 thebackground 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

Page 11 of 16

source must match MSL’s observations, the ~ northernemission sources must have higher emission fluxes than
the southern ones. During daytime, however, when thewind regime changes and upslope anabatic winds dominate this region on the slopes of Aeolis Mons, the northern emission sources may produce high and persistentmethane 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) anddaytime (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

Fig. 7 Surface methane mixing ratio (output from the first vertical model level at 10 m) in the vicinity of the MSL in logarithmic scale, both during
nighttime (0:00 h LTST) and daytime (12:00 h LTST) produced by two different sources at different locations (locations 18 and 47) relative to theMSL (represented by a black dot). The emission area is represented by a square. Both emission sources present similar methane abundances, inaccordance with MSL’s observations (~ 1 ppbv or less, Webster et al. 2018) during nighttime; however, the unfavorable air mass transport during
the night for the northwestern emission source (location 47, bottom panels) requires a much higher emission flux to match MSL’s observations. It
implies that the extremely favorable anabatic flows blowing over the slopes of Aeolis Mons, and the development of the planetary boundary layer,
during daytime, produce a sustainable increase in the methane abundance at MSL’s location to values of the order of 100 ppbv, which is not what
has been observed at all (Webster et al. 2015; 2018)

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

Page 12 of 16

Fig. 8 Simulated methane spikes during the diurnal cycle at MSL’s general location, in the first prognosed layer (~ 10 m), for the south andsoutheast emission sources (see text). Numbers around the spikes are referred to the labels established for each emission location as presented inFig. 6. Almost every source location in the vicinity of MSL produces methane spikes at MSL’s general location. Locations within фR2 are represented
in green color. Western locations (фR1)are represented in blue color and eastern ones (фR3)in red colorthat the extremely favorable anabatic flows over theslopes of Aeolis Mons, together with the development of
the planetary boundary layer during daytime, produce asustainable high abundance of methane at MSL’s generallocation, to values of the order of 100 ppbv, which is notwhat has been observed at all (Webster et al. 2015; 2018).
Figure S3 shows the simulated methane abundance during the diurnal cycle at MSL’s general location in the firstprognosed 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 willtend to produce present higher abundances during theafternoon, 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 variabilityof 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 lessthan 1 ppbv between 0–2 LTST (in accordance with thebackground 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 ­103 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 theslopes of Aeolis Mons during nighttime (Fig. 7). The timing and intensity of these spikes will also vary dependingon many factors, such as the season and meteorologicalconditions.
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 winddirection in the first prognosed layer at MSL’s location.

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

Page 13 of 16

Fig. 9 Surface methane distribution, in logarithmic scale, in the vicinity of the MSL in several time instants covering a methane spike presented inFig. 8 and produced by a source location centered at 4.7 deg of latitude, 137.4 deg of longitude, in time intervals of 1 h starting on local true solar
time (LTST) = 0.5 h. MSL’s general location is also indicated by a dot

Given the proximity between the emission source andthe 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 comingfrom 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,

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

Page 14 of 16

Fig. 10 Volume mixing ratio, horizontal wind speed (vh), and wind direction (Ψd) in the first prognosed layer at MSL’s location, covering themethane 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 fromwhich the wind is blowing, following the standard meteorological convention (0° < Ψd ≤ 360°, with the origin pointing to the north)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 vicinityof 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 andthus 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 DISVERMAR simulations strongly suggest that the variability ofwinds in the vicinity of the MSL is sufficient to produce
strong methane spikes at MSL’s general location from anearby 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 invocationof variable emission fluxes, eventual emissions, and/or microscale flows could also dramatically increase theobserved 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 andcould develop in the worst-case scenario of a constantemission 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 thenorthwest slopes of Aeolis Mons).
As discussed above, the methane spikes observed at
MSL’s location would be the result of air masses withhigh 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 asindicative of an abundance of methane in the near-surface atmosphere and beyond. The maximum levels of
methane in the atmosphere will strongly depend on thelocation 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

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

abundance detected at MSL’s location (Fig. 10). However, the DISVERMAR simulations show that the mixing ratio decreases for that event (which occurs during
nighttime) to less than 1­ 0−3 ppbv above the PBL. When
analyzing other spikes, including those that occurred
during daytime, it can be noted that the methane levels
never exceed 0.01 ppbv. This reduction in methane abundance above the PBL, together with the fact that spikes
are quickly dispersed in the atmosphere in a few hours
or less (Fig. 8), could significantly hinder their potential
detection from orbit.

4 Conclusions
A new atmospheric dispersion model for the Martian
atmosphere has been presented. This model, DISVERMAR, can be classified as a three-dimensional local-toregional-scale Eulerian atmospheric dispersion model. It
has been developed to simulate emissions to the atmosphere from particular locations or regions, including
chemical loss and predefined deposition rates. The model
can deal with complex topography and non-uniform
grids.
The model could be used for several applications
related to atmospheric transport and chemistry on Mars,
in a context for which limited empirical data are available due to difficulties in delivering and operating any
instrument on the Martian surface. This model could bea valuable tool in terms of assisting in the interpretationof such observations, complementary to the modeling
efforts performed to date in the Mars atmosphere.
As a case study, the model is applied to the methane
abundance detected by NASA’s Mars Science Laboratory in Gale Crater, Mars, to determine whether acontinuous release of methane in the vicinity of MSL
in Gale Crater could produce methane spikes. The DISVERMAR simulations strongly suggest that the variability of winds in the vicinity of the MSL is sufficient
to produce strong methane spikes at MSL’s generallocation from a nearby emission source, even without
invoking sol-to-sol variability in winds, variable emission fluxes or eventual emissions into the atmosphere,or microscale flows not resolved in the meteorological
data. In addition, the DISVERMAR simulations suggest
that, in the absence of a more complex emission pattern, north–northwest source locations are incompatible with the MSL’s observations, favoring a hypotheticalemission source at southern–southeastern locations.
Abbreviations
DISVERMAR: Three-dimensional atmospheric dispersion versatile model for
Mars; MSL: Mars Science Laboratory; PBL: Planetary boundary layer; PDS:
NASA’s Planetary Data System.

Page 15 of 16

Supplementary Information
The online version contains supplementary material available at https://​doi.​
org/​10.​1186/​s40645-​021-​00445-4.
Additional file 1. Supplementary material.
Acknowledgements
The author thanks C.E. Newman and M.I. Richardson for retrieving theMarsWRF model code and their valuable inputs to perform the mesoscale
simulations necessary to generate the meteorological fields for this study. The
author also thanks three anonymous reviewers for their helpful comments,
which have improved this manuscript.
Authors’ contributions
The author read and approved the final manuscript.
Availability of data and materials
The topographic maps used to input the model are obtained from the NASA’s
Planetary Data System (PDS) (https://​pds-​geosc​iences.​wustl.​edu/​mgs/​mgs-m-​
mola-5-​megdr-​l3-​v1/​mgsl_​300x/). Model outputs are available upon request.

Declarations
Competing interests
The authors declare that they have no competing interests.
Received: 1 April 2021 Accepted: 27 August 2021

References
Draxler RR, Hess GD (1998) An overview of the HYSPLIT_4 modellingsystem for trajectories, dispersion, and deposition. Austr Meteorol Mag
47:295–308
Durran DR (2010) Numerical methods for fluid dynamics. Springer, New York
Giuranna M et al (2019) Independent confirmation of a methane spike on
Mars and a source region east of Gale Crater. Nat Geosci 12:326–332
Grell GA, Knoche R, Peckham SE, McKeen SA (2004) Online versus offline air
quality modeling on cloud-resolving scales. GRL 31:L16117
Hirsch C (2007) Numerical computation of internal and external flows.
fundamentals of computational fluid dynamics, vol 1. JohnWiley & Sons,
HobokenHong SY, Pan HL (1996) Nonlocal boundary layer vertical diffusion in a
medium range forecast model. Mon Weather Rev 124:2322–2339
Hong SY et al (2006) A new vertical diffusion package with an explicit treat‑ment of entrainment processes. Montly Weather Rev 134:2318–2341
Ishikawa H (1993) Mass-consistent wind model as a meteorological preproces‑
sor for tracer transport models. J Appl Meteorol 33:733–743
Korablev O, Vandaele AC, Montmessin F et al (2019) No detection of methane
on Mars from early ExoMars Trace Gas Orbiter observations. Nature
568:517–520. https://​doi.​org/​10.​1038/​s41586-​019-​1096-4
Leelossy A, Molnár F, Izsák F, Havasi A, Lagzi I, Mészáros R (2014) Dispersion
modeling of air pollutants in the atmosphere: a review. Cent Eur J Geosci
6(3):257–278
Lefèvre F (2019) The enigma of methane on mars. In: Cavalazzi B, Westall F(eds), Biosignatures for astrobiology. Springer
Lefèvre F, Forget F (2009) Observed variations of methane on Mars unex‑
plained by known atmospheric chemistry and physics. Nature. https://​doi.​org/​10.​1038/​natur​e08228
Louis JF (1979) A parametric model of vertical eddy fluxes in the atmosphere.
Bound-Layer Meteorol 17:187–202
Moores JE, King PL, Smith CL, Martinez GM, Newman CE, Guzewich SD et al(2019) The methane diurnal variation and microseepage flux at Gale Cra‑ter, Mars as constrained by the ExoMars Trace Gas Orbiter and Curiosity
observations. Geophys Res Lett 46:9430–9438
Newman CE, Gómez-Elvira J, Marín M, Navarro S, Torres J, Richardson MI,
Battalio JM, Guzewich SD, Sullivan R, de la Torre-Juárez M, Vasavada AR,

Viúdez‑Moreiras ﻿Prog Earth Planet Sci(2021) 8:53

Bridges NT (2017) Winds measured by the Rover Environmental Monitor‑
ing Station (REMS) during the Mars Science Laboratory (MSL) rover’s
Bagnold Dunes Campaign and comparison with numerical modelingusing MarsWRF. Icarus 291:203–231
Newman CE et al (2021) Multi-model meteorological and aeolian predictionsfor Mars 2020 and the Jezero Crater Region. Space Sci Rev 217(20)Pla-Garcia J, Rafkin SCR, Karatekin Ö, Gloesener E (2019) Comparing MSL
Curiosity rover TLS-SAM methane measurements with Mars Regional
Atmospheric Modeling System atmospheric transport experiments. J
Geophys Res Planets 124:2141–2167
Rafkin SC, Haberle RM, Michaels TI (2001) The Mars regional atmosphericmodeling system: model description and selected simulations. Icarus
151:228–256
Rafkin SC, Pla-Garcia J, Kahre M, Gomez-Elvira J, Hamilton VE, Marín M, Navarro
S, Torres J, Vasavada A (2016) The meteorology of Gale Crater as deter‑
mined from Rover Environmental Monitoring Station observations andnumerical modeling. Part II: Interpretation. Icarus 280:114–138
Richardson MI et al (2007) PlanetWRF: a general purpose, local to global
numerical model for planetary atmospheric and climate dynamics. JGRPlanets 112(E9):E09001
Sankaranarayanan S, Shankar NJ, Cheong HF (1998) Three-dimensional finite
difference model for transport of conservative pollutants. Ocean Eng
25(6):425–442
Smagorinsky J (1963) General circulation experiments with the primitive equa‑tions: 1. The basic experiment. Mon Weather Rev 91:99–164
Smolarkiewicz PK, Margolin LG (1998) MPDATA: a finite-difference solver for
geophysical flows. J Comput Physics 140:459–480
Spiga A, Forget F (2009) A new model to simulate the Martian mesoscale andmicroscale atmospheric circulation: validation and first results. J Geophys
Res Planets 114:E02009
Steele LJ, Balme MR, Lewis SR, Spiga A (2017) The water cycle and regolithatmosphere interaction at Gale crater. Mars Icarus 280:56–79
Toigo AD et al (2012) The impact of resolution on the dynamics of the Martian
atmosphere: varying resolution studies with the MarsWRF GCM. Icarus
221(1):276–288

Page 16 of 16

Troen I, Mahrt LA (1986) A simple-model of the atmospheric boundary-layer
sensitivity to surface evaporation. Boun Layer Meteorol 37(1–2):129–148
Tyler Jr D, Barnes JR (2013) Mesoscale modeling of the circulation in the Gale
Crater region: an investigation into the complex forcing of convective
boundary layer depths. Int J Mars Sci Explor 8:58–77
Tyler D Jr, Barnes JR, Haberle RM (2002) Simulation of surface meteorology at
the Pathfinder and VL1 sites using a Mars mesoscale model. J Geophys
Res 107(E4):5018
Viúdez-Moreiras D, Gómez-Elvira J, Newman CE, Navarro S, Marin M, Torres J,
de la Torre M, the MSL team (2019) Gale surface wind characterization
based on the Mars Science Laboratory REMS Dataset. Part I: wind retrieval
and Gale’s Wind speeds and directions. Icarus 319:909–925
Viúdez-Moreiras D, Arvidson RE, Gómez-Elvira J, Webster C, Newman CE,
Mahaffy P, Vasavada AR (2020a) Advective fluxes in the Martian regolith
as a mechanism driving methane and other trace gas emissions to theatmosphere. Geophys Res Lett 47:e2019GL085694
Viúdez‐Moreiras D et al (2020b) Effects of a large dust storm in the near‐surface atmosphere as measured by InSight in Elysium Planitia, Mars.
Comparison with contemporaneous measurements by Mars Science
Laboratory. JGR-Planets 125:e2020JE006493. https://​doi.​org/​10.​1029/​
2020J​E0064​93
Webster CR et al (2015) Mars methane detection and variability at Gale crater.
Science 347:415–417
Webster CR et al (2018) Background levels of methane in Mars’ atmosphereshow strong seasonal variations. Science 360(6393):1093–1096
Yung Y et al (2018) Methane on mars and habitability: challenges andresponses. Astrobiology. https://​doi.​org/​10.​1089/​ast.​2018.​1917

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in pub‑lished maps and institutional affiliations.