Sub‐Diurnal Methane Variations on Mars Driven by Barometric Pumping and Planetary Boundary Layer Evolution

In recent years, the Tunable Laser Spectrometer within the Sample Analysis at Mars (TLS‐SAM) instrument on board the Mars Science Laboratory (MSL) Curiosity rover has detected methane variations in the atmosphere at Gale crater. Methane concentrations appear to fluctuate seasonally as well as sub‐diurnally, which is difficult to reconcile with an as‐yet‐unknown transport mechanism delivering the gas from underground to the atmosphere. To potentially explain the fluctuations, we consider barometrically induced transport of methane from an underground source to the surface, modulated by temperature‐dependent adsorption. The subsurface fractured‐rock seepage model is coupled to a simplified 1‐D atmospheric mixing model to provide insights on the pattern of atmospheric methane concentrations in response to transient surface methane emissions, as well as to predict sub‐diurnal variation in methane abundance for the northern summer period, which is a candidate time frame for a MSL Curiosity sampling campaign. Our analysis suggests that there is a lower limit to the subsurface fracture density that can produce the observed methane patterns, below which the atmospheric methane variations would be out of phase with the observations. The best‐performing model scenarios indicate a significant, short‐lived methane pulse just prior to sunrise, the detection of which by TLS‐SAM would be a potential indicator of the contribution of barometric pumping to Mars' atmospheric methane variations.

• Barometrically driven atmospheric methane abundance timing controlled by fracture topology and atmospheric dynamics • There is a lower limit to fracture density (0.01%) that can produce observed methane patterns • A late morning or early evening Tunable Laser Spectrometer within the Sample Analysis at Mars sample could constrain diurnal methane pattern and transport processes Supporting Information: Supporting Information may be found in the online version of this article.
In addition to apparent seasonal fluctuations in methane, some evidence suggests that atmospheric methane varies on a sub-diurnal time scale as well (Webster et al., 2021).TLS-SAM primarily conducts experiments at night due to mission operational constraints (primarily the thermal requirements of TLS-SAM sample handling system; Pla-García et al. ( 2019), but also the large power requirement of the instrument, which limits the use of other instruments), and in fact all TLS-SAM detections of methane thus far have been from nighttime measurements.Two lone non-detections in 2019 were reported from daytime measurements (Webster et al., 2021) during northern summer at Gale crater on sols 2615 and 2644 (approximate location: −4.73°S 137.38°E).These sporadic non-detections occurred on either side of a normal background methane value collected at night, implying a diurnal to sub-diurnal variability in atmospheric methane.Confirming and characterizing this apparent diurnal variability of methane has been highlighted by the TLS-SAM team as the next key step to understanding methane abundance and circulation at Gale crater (Moores, Gough, et al., 2019;Webster et al., 2021).
The primary goal of this work is to facilitate the science goals of ongoing and future atmospheric sample collection missions on Mars by determining an optimal intra-sol timing for collection.It is possible that reduced electrical power over time in conjunction with TLS-SAM pump life will likely place limits on scientific operations.It is therefore important to maximize the scientific return of whatever remaining TLS-SAM measurements there may be, especially with regard to characterizing the apparent diurnal variability in methane.Recent models (Giuranna et al., 2019;Luo et al., 2021;Pla-García et al., 2019;Viúdez-Moreiras, 2021;Viúdez-Moreiras et al., 2021;Webster et al., 2015Webster et al., , 2018a;;Yung et al., 2018) suggest a local source of methane within Gale crater, with circulation trapping methane at night and dissipating it during the day.Characterizing the diurnal variability of methane provides insight into the underlying mechanisms driving the methane fluctuations.Based on the apparent seasonal trend with infrequent measurements, the logical time of year to make relevant measurements is in the northern Summer period between solar longitude (L s ) 120-140°, coincident with the time of year of the previous measurements indicating diurnal variations.At the time of writing, this period is approaching in the months of September-October 2023, which presents an ideal opportunity for collecting in situ atmospheric methane data at Gale crater.
Running TLS-SAM experiments at strategically optimal times will improve the probability of gathering useful atmospheric data to answer key questions about methane at Gale crater.Numerical models of methane emissions and mixing within the atmosphere have the potential to inform this goal of determining ideal times to collect samples.The general consensus in the planetary science community is that if methane is present in Mars' atmosphere, its source is most likely located underground (Carrier et al., 2020; The National Academies of Sciences, Engineering & Medicine, 2022).This presents the question of how methane from deep underground can reach the surface rapidly enough to generate the observed short-term atmospheric variations.Some of the possibilities that have been proposed include: a relatively fast methane-destruction mechanism, modulation mechanisms that change the amount of free methane in the atmosphere and near-surface (e.g., regolith adsorption), and rapid transport mechanisms capable of delivering gases from depth (e.g., barometric pumping).This paper focuses on the latter two of these, and uses simulations driven by high resolution pressure and temperature data as forcing in order to provide insight on the timing of sub-diurnal methane fluxes driven by barometric pumping. 10.1029/2023JE008043 3 of 26 Barometric pumping is an advective transport mechanism wherein atmospheric pressure fluctuations greatly enhance vertical gas transport in the subsurface (Nilson et al., 1991).Low atmospheric pressure draws gases upwards from the subsurface, with air and tracer movement taking place primarily in the higher-permeability fractures rather than the surrounding, relatively low-permeability rock matrix (Figure 1).High atmospheric pressure pushes gases deeper into the subsurface, with some molecules diffusing into the rock matrix, in which the barometric pressure variations do not propagate efficiently.Over multiple cycles of pressure variations, this fracture-matrix exchange produces a ratcheting mechanism (Figure 1) that can greatly enhance upward gas transport relative to advection and diffusion alone (Harp et al., 2018;Massmann & Farrier, 1992;Neeper & Stauffer, 2012a;Nilson et al., 1991;Takle et al., 2004).Barometric pumping has been studied in a variety of terrestrial contexts on Earth, such as: CO 2 leakage from carbon sequestration sites (Carroll et al., 2014;Dempsey et al., 2014;Pan et al., 2011;Viswanathan et al., 2008) and deep geological stores (Etiope & Martinelli, 2002;Rey et al., 2014), methane leakage from hydraulic fracturing operations (Myers, 2012), radon gas entry into buildings (Tsang & Narasimhan, 1992), contaminant monitoring (Stauffer et al., 2018(Stauffer et al., , 2019)), and radionuclide gas seepage from underground nuclear explosions and waste storage facilities (Bourret et al., 2019(Bourret et al., , 2020;;Carrigan et al., 1996Carrigan et al., , 1997;;Harp et al., 2020;Jordan et al., 2014Jordan et al., , 2015;;Sun & Carrigan, 2014).In the context of Mars, barometric pumping in fractures was first hypothesized as a potentially effective transport mechanism for underground methane by Etiope and Oehler (2019).Although a few modeling papers (Klusman et al., 2022(Klusman et al., , 2024;;Viúdez-Moreiras et al., 2020) have investigated barometric pumping in the context of methane transport on Mars, our recent paper (Ortiz et al., 2022) is, to our knowledge, the first to consider the explicit role of subsurface fractures and the ratcheting mechanism.In that paper, we demonstrated that barometric pumping in fractured rock is capable of producing significant surface fluxes of methane from depths of 200 m, and the timing and magnitude of those fluxes were reasonably consistent with the timing of high-methane periods measured by MSL Curiosity.The emphasis on timing in that paper was on reproducing the observed seasonality of surface fluxes.We highlighted in our discussion that the timing of surface fluxes could be further modulated by processes that retard gas transport and therefore included adsorption in shallow regolith to produce a more complete transport model in this work.
Adsorption is a reversible phenomenon in which gas or liquid molecules (the "adsorbate") adhere to the surface of another material (the "adsorbent").Particle transport (e.g., methane) through porous media (e.g., martian

10.1029/2023JE008043
4 of 26 regolith), is retarded by adsorption onto the pore walls.Adsorption is aided by adsorbents with high specific surface area, which have more sites onto which the particles can adsorb.It is believed that much of the martian regolith consists of fine mineral dust particles (Ballou et al., 1978;Chevrier & Mathé, 2007;Rampe et al., 2020), which have a large specific surface area (Meslin et al., 2011), making the regolith relatively amenable to adsorption.Furthermore, adsorption reactions are generally temperature-dependent, with lower temperatures favoring adsorption and higher temperatures favoring desorption.Specifically, both the rate of adsorption and the equilibrium surface coverage are higher at lower temperatures for many systems (Adamson, 1979;Gough et al., 2010;Himeno et al., 2007;Pick, 1981).
Several previous papers have investigated whether the temperature dependence of regolith adsorption could explain the seasonal variations in methane in the martian atmosphere because of this temperature dependence.
Work by Gough et al. (2010) used laboratory-derived constants to determine the seasonal variation of methane across Mars due to adsorptive transfer to and from the regolith.Extrapolating to martian ground temperatures, the adsorption coefficient measured for methane gas was relatively low, though the authors concluded that the mechanism could still be capable of contributing to rapid methane loss.Meslin et al. (2011) used a global circulation model to determine the seasonal variation of methane due to adsorptive transfer into and out of the regolith, finding that at Gale's latitude, this seasonal variation in methane was less than a few percent, and therefore not likely the cause of the methane fluctuations.Another paper (Moores, Gough, et al., 2019) investigated regolith adsorption, but with methane provided by a shallow (30 m) microseepage source, and found that their one-dimensional adsorptive-diffusive numerical model was able to produce the observed seasonal variation.
More recently, research by Klusman et al. (2022) followed the analysis of Moores, Gough, et al. (2019) pertaining to adsorption, while also considering the role of barometric pumping as the primary transport mechanism for the shallow subsurface, and were able to produce the seasonal variation of methane when invoking high regolith permeabilities (10 −10 m 2 ).
In this manuscript, we consider the barometrically induced transport of a subsurface methane source to the planet's surface that is modulated by temperature-dependent adsorption/desorption.Our two-dimensional simulations consider the explicit role of discrete, interconnected fractures in promoting advective transport, with additional seasonal modulation provided by temperature-dependent regolith adsorption.To elucidate the effects of subsurface architecture (i.e., the degree of fracturing in the rock, quantitatively represented in terms of fracture density, and defined as the ratio of fracture volume to total bulk rock volume), we simulate gas flow and transport through rocks with fracture density ranging from 0% (unfractured), to 0.035% (highly fractured), with the chosen range derived from rover photographs at Gale crater (Text S6.2 in Supporting Information S1).The subsurface seepage model is coupled to a 1-D atmospheric diffusion and mixing model to provide generalized insights on the pattern of atmospheric concentrations of methane in response to transient surface methane emissions, as well as to predict sub-diurnal variation in methane abundance for the northern summer season.In reality, the atmospheric circulation of methane involves complex three-dimensional flow and transport.Nonetheless, our 1-D atmospheric mixing model representing the bulk behavior of gases in the planetary boundary layer (PBL) column can provide a reasonable approximation of the temporal evolution of methane abundance in the near-surface atmos phere in response to a fluctuating ground surface emission.

Methods: Fractured-Rock Heat and Mass Transport Simulations With Coupled Atmospheric Mixing
We used fractured-rock heat and mass transport simulations to determine the approximate timing of transient methane surface fluxes driven by barometric fluctuations throughout the Mars year.Calculations are performed within the Finite-Element Heat and Mass (FEHM) simulator, a well-tested multiphase code (Zyvoloski et al., 1999(Zyvoloski et al., , 2017(Zyvoloski et al., , 2021)).FEHM has been used extensively in terrestrial barometric pumping studies (Bourret et al., 2019(Bourret et al., , 2020;;Jordan et al., 2014Jordan et al., , 2015;;Neeper & Stauffer, 2012a, 2012b;Stauffer et al., 2019), and was previously modified by the author to adapt to conditions at Mars in a related paper examining barometric pumping of methane (Ortiz et al., 2022).We have made a simplifying assumption that there is no liquid water in the domain, which is a reasonable assumption over the range of subsurface depths we consider (to 500 m).The base of the cryosphere at Gale crater-below which temperatures exceed 0°C and are capable of supporting liquid water-is estimated to be at ∼3,850 m depth (Klusman et al., 2022), though liquid brines could exist at somewhat shallower depths due to a lower freezing point.If liquid water were present, it would reduce the available air-filled 10.1029/2023JE008043 5 of 26 porosity (as ice) and cause temporary immobile storage due to phase partitioning of gas-phase constituents into the liquid.The effect of our assumption is therefore an upper bound on advective gas transport, excepting retardation due to adsorption.Gravity and atmospheric gas properties are modified for this study to replicate Mars conditions.
Our simulations require several steps: (a) heat flow simulations to generate the subsurface temperature profiles, (b) subsurface mass flow and transport simulations of Mars air and methane driven by barometric fluctuations, with regolith adsorption terms dictated by the subsurface temperature changes from step 1, and (c) 1-D atmospheric mixing of methane emitted from the subsurface into a transient PBL column in order to calculate methane mixing ratios.Parameter value ranges for each model component are provided in Table 1.
Initial testing of a coupled energy and mass transport model indicated that due to conduction dominance (the fracture volume fraction is very small), the temperature field can be adequately described using a decoupled 1-D conductive heat transfer model.We therefore run the heat transport simulations to generate time-dependent temperature profiles with depth.We then run the 2-D, fractured-rock mass flow and transport simulations to calculate the fluxes of martian air and methane driven by barometric fluctuations.The flow model assumes isothermal conditions, while the transport model considers temperature variations in its calculation of adsorption coefficients.The assumption of isothermal conditions in the flow model is justified based on verification tests, which indicated that the martian air flow properties were not significantly modified by ignoring temperature effects (Text S2.4 in Supporting Information S1).Mass flow and transport equations in the fractures are coupled to those in the rock matrix to simulate the overall behavior of gases in fractured rock.These approaches are standard in subsurface hydrogeology-the governing equations and computational approach are described in detail below in Section 2.2.Finally, we simulate the 1-D atmospheric mixing of methane by coupling the surface methane emissions to a diffusive transport model within a PBL column of time-varying height (Section 2.4).This step allows us to infer atmospheric methane concentrations generated in response to the time history of surface fluxes emitted in the subsurface seepage model.

Heat Flow Model
Although the mass flow and transport simulations use a 2-D domain, we found that simple matrix conduction dominated over fracture convection, which had a negligible influence over subsurface temperatures (Text S2.3 in Supporting Information S1), justifying the simulation of transient subsurface heat transport using a 1-D model.The 1-D approach also facilitates computational efficiency due to the high degree of mesh refinement required to accurately simulate subsurface temperatures (Text S2.1 in Supporting Information S1).The single-phase heat conduction equation (Carslaw & Jaeger, 1959) is as follows: where , where κ is the thermal conductivity of the material [W m −1 K −1 ], c is the specific heat capacity [J K −1 kg −1 ], and ρ r is the density of the material [kg m −3 ]).

Boundary and Initial Conditions: Heat Flow Model
We prescribe an initial surface temperature of −46.93°C (226.22K), which is the mean surface temperature at Gale crater (Klusman et al., 2022) as obtained from the Rover Environmental Monitoring Station (Gómez-Elvira, 2019).Ground surface temperatures fluctuate about this mean value, so this temperature is also used as the reference temperature for CO 2 properties (Mars atmosphere is 95% CO 2 ) in the equation of state for the mass flow model.At ground surface, we prescribe temperature as a time-varying Dirichlet boundary condition.We generated a synthetic temperature record representative of the surface temperatures collected by REMS.We extended the time series of generated temperatures so that the simulations can spin up with a sufficiently long record.At the bottom of the domain, we prescribe temperature as a constant Dirichlet boundary condition assigned based on the geothermal gradient and depth of the domain being considered.

Subsurface Mass Flow and Methane Transport Model
The flow and transport simulations are set up similarly to those presented in Ortiz et al. (2022), with some exceptions listed in the subsequent paragraph.Transient barometric pressures are prescribed at the ground surface and serve as the primary forcing condition.Methane is produced at a constant rate within a 5-m-thick zone at variable depths within the domain depending on the scenario, and is allowed to escape the subsurface domain only at the ground surface boundary.
In contrast to the simulations previously published (Ortiz et al., 2022), these simulations include the effects of temperature-dependent regolith adsorption.We model regolith adsorption as a Langmuir adsorption process, following Gough et al. (2010) and Moores, Gough, et al. (2019), described in greater detail in the following subsection (Section 2.2.1).The martian air, and the tracer gas (methane, CH 4 ) have properties consistent with the mean ambient pressure and temperature conditions at Gale crater.
As in the heat flow model, we extracted the dominant frequency and amplitude components of the barometric pressure record collected by the Curiosity MSL-REMS (Gómez-Elvira, 2019) using Fourier analysis.We then generated a synthetic barometric pressure record using these components (Text S1 in Supporting Information S1), which allows us to treat the problem in a more general way while extending the time series of the pressure forcing to achieve cyclical steady-state in the surface fluxes.

Governing Equations and Boundary Conditions
The governing flow equations for single-phase flow of martian air in the fracture network are given by:  is flux into the fracture.Note that Equation 2 is an aperture-integrated two-dimensional equation for fracture flow and Equation 3 is the local cubic law for laminar fracture flow (Zimmerman & Bodvarsson, 1996).
Governing equations for flow in the matrix are given by: where ∇ is the 3-D gradient operator, ϕ is the porosity [dimensionless; m 3 /m 3 ], k m is matrix permeability [m 2 ], and P m is the air pressure in the rock matrix [Pa].Note that P f = P m on the fracture-matrix interface (I), and the pressure gradients ∇P m at the fracture-matrix interface control the interfacial flux on the right-hand side of Equation 5. We make the assumption that the bulk movement of air through the rock matrix behaves according to Darcy's law (Equation 5).In the case of a low-permeability rock matrix, the pressure gradients and fluxes induced in the matrix by barometric pressure variations are typically small.
The governing equations for transport of a tracer gas (e.g., methane) in a fracture are given by: where Governing equations for transport in the rock matrix with adsorption are given by: where ρ r is the rock density [kg m −3 ], s is the adsorbed concentration [mol kg r ], k eq is the Langmuir equilibrium distribution coefficient [kg v mol −1 ], and    is the tracer source in the matrix [mol m −3 s −1 ], and C f = C m on the fracture-matrix interface.The distribution coefficient k eq is temperature-dependent, and its formulation in the model is described in more detail in Section 2.2.1.
The flow and transport simulations use martian air and methane properties consistent with the mean surface temperature at Gale crater (−46.93°C).The bottom of the domain is a no-flux boundary.The left and right lateral boundaries are no-flux boundaries.The top-surface boundary is forced by the synthetic barometric pressure record we generated using frequency and amplitude components representative of the pressure record collected by MLS-REMS (see Text S1 in Supporting Information S1).Vapor-phase methane and martian air are allowed to escape the domain from the top boundary.We prescribe a continuous methane production rate (9.6 × 10 −7 mg CH 4 m −3 sol −1 ) within a 5-m-thick zone at the bottom spanning the lateral extent of the domain (Figure 2a).This rate is consistent with liberal estimates of the maximum methane production rate by serpentinization reactions on Mars (Stevens et al., 2015) in addition to measurements of methanogenic microbes at depth in Mars-analog 10.1029/2023JE008043 8 of 26 terrestrial settings (Colwell et al., 2008;Onstott et al., 2006).Methane-clathrate-hydrate (MCH) reservoirs have also been hypothesized as a potential source of methane (Gloesener, 2019;Gloesener et al., 2021;Pla-García et al., 2019), though the associated release of methane would likely be episodic and therefore not well-represented by a continuous subsurface source.Our model assumes direct source rock-to-seepage pathway similar to that described in Etiope et al. (2013), rather than a source-reservoir-seepage system.We considered a range of methane source depths (labeled as "methane production zone" in Figure 2a) from 5 to 500 m below ground surface.For source depths ≤200 m, a standard 200 m depth model domain was used.For the cases with source depth 500 m, we used a model domain of depth 500 m.
The flow and transport simulations are performed in three steps: (a) initialization, (b) "spin-up," and (c) the main flow and transport runs.We initialize the flow model using a constant surface pressure for 10 8 years to create a martian air-static equilibrium gradient throughout the subsurface.This duration is chosen because it is sufficiently long; after 10 8 years, we can confidently assert that there are negligible changes to the modeled subsurface pressures.The initialization simulation is run without methane in the domain.We used this martian air-static pressure equilibrium as the initial state for the flow and transport simulations.
We then run a spin-up simulation lasting 50100 sols, equivalent to 75 Mars Years (MY).The purpose of the spin-up simulation is to establish the memory of surface pressure and temperature fluctuation periodicity in the subsurface.Additionally, it allows for the methane generated in the source zone to sufficiently populate the subsurface and reach a cyclical steady-state in terms of surface flux.We verify in each case that the system has reached a cyclical steady-state equilibrium by identifying a linear trend in cumulative surface mass outflow that does not change year-over-year.In practice, 75 MY is more than enough time for this equilibrium to develop-linear regression on annual mass outflow for the last 15 MY yields a slope of 1.2 × 10 −13 kg/MY, representing a mean annual change of <0.00002% in the worst case examined, where a slope of 0.0 represents exactly no change year-over-year.The domain is initially populated with a uniform concentration of methane gas (C 0 = 9.6 × 10 −5 mol  kg −1  ) to allow the subsurface to more efficiently reach a quasi-equilibrium by pumping out excess methane from the system in the early stages of the simulation.Adsorbed methane concentration is initially zero everywhere.Finally, we run the flow and transport simulations starting from the conditions established in the initialization and spin-up runs.To be conservative, the final simulations are run for 75 MY, and implement the same mechanisms as the spin-up simulations.

Temperature-Dependent Langmuir Adsorption Model Implementation
The Langmuir adsorption isotherm can be used to adequately describe the adsorption-desorption process on Mars analogs (Moores, Gough, et al., 2019).This is partly due to the fact that for methane at the low average temperatures on Mars, the surface coverage θ (i.e., the fraction of the adsorption sites occupied at equilibrium), is estimated to be quite low (of order 10 −10 ), so that the Brunauer-Emmett-Teller (BET) formulation is unnecessary.The equilibrium rate constant k eq (ratio of sorbed phase to gas phase concentration) for the adsorption isotherm is defined as: where k eq is the equilibrium rate constant [kg v mol −1 ], C i is the vapor-phase concentration of the tracer gas i (which in this case can be assumed to be methane) [mol ], P i is the partial pressure of the tracer gas, R a and R d are the absolute rates of adsorption and desorption, and θ is the surface coverage [m 2 /m 2 ].The equilibrium surface coverage θ eq can be approximated using the k eq at a given partial pressure of methane   4 (or concentration C CH4 ) and temperature T: The equilibrium constant can be adapted to a partial-pressure basis: where γ is the uptake coefficient (determined experimentally), η is the evaporation coefficient, ν is the mean molecular speed,  ML 4 is the number of methane molecules per m 2 of adsorptive surface required to form a monolayer, h is Planck's constant, and k B is Boltzmann's constant.For a reversible adsorption-desorption process, the evaporation and uptake coefficients are approximately equal (η ≈ γ) and drop out of Equation 10 (Gough et al., 2010).The monolayer coverage variable  ML 4 is calculated as 5.21 × 10 18 molecules m −2 based on the size of an adsorbed methane molecule (19.18 Å) (Chaix & Dominé, 1997).
Implementation of temperature-dependent adsorption in FEHM is relatively straightforward.Because the simulation time is quite long, it is more computationally efficient to sequentially couple the temperature field to the mass flow and transport simulations.We performed several verification tests to ensure that the martian air flow properties were not significantly modified by ignoring temperature effects (Text S2.4 in Supporting Information S1).
Using the subsurface temperatures acquired from the heat flow simulation, at each node we assign a distribution coefficient for the adsorption reaction that varies with depth and time.In this way, the flow and transport simulations are non-isothermal insofar as they account for temperature-dependent adsorption.Gough et al. (2010) reported on the results of laboratory studies of methane adsorption onto JSC-Mars-1, a martian soil simulant, and determined the ΔH for methane adsorption using experimentally determined values of the uptake coefficient (γ), which is the ratio between the adsorption rate and gas molecule collision rate.They found that the observed energy change, ΔH obs , for methane adsorption onto JSC-Mars-1 is 18 ± 1.7 kJ mol −1 .
Although not identical to the overall adsorption enthalpy, ΔH tot , it is a lower limit for this process that is similar to the overall adsorption enthalpies reported by others for similar systems (Gough et al., 2010).From this, we have calculated the values of k eq as it varies with temperature and tabulated them into a format usable by FEHM.
Because the surface temperature perturbations do not propagate very far into the subsurface (Figure S7 in Supporting Information S1), we actively calculate the time-dependent Langmuir distribution coefficient k eq only 10.1029/2023JE008043 10 of 26 for the upper 5 m of regolith, and we assign a temporally and spatially constant average k eq value for the remainder of the subsurface.This has the added benefit of reducing the computational costs of the simulation.
While we considered only the adsorption and tracer transport of methane gas in our model, the presence of other gas species could certainly affect the overall amount of methane adsorbed via competition for adsorption sites on the regolith.It is true that the CO 2 comprising most of Mars' air could adsorb to the regolith as well; however, the difference in adsorption enthalpies of methane and CO 2 leads to a three order-of-magnitude difference in their adsorption timescales, favoring methane, so it is unlikely to have a noticeable effect.Some researchers Moores, Gough, et al. (2019) have considered competitive adsorption between methane and water vapor in the uppermost 1 m of the regolith, though this was outside the scope of the current work.The absence of competing gas species adsorbates in our model means in essence that we are representing the conservative upper bound of adsorption effects within a barometrically driven transport system.

Geologic Framework and Numerical Mesh
We assigned the background rock matrix a porosity (ϕ m ) of 35%, which is in the range estimated by Lewis et al. (2019) based on consideration of the low bedrock density at Gale crater.We set the background rock permeability (k m ) to 1 × 10 −14 m 2 (0.01 Darcies).This is slightly more permeable than the conservative 3 × 10 −15 m 2 prescribed by previous research modeling hydrothermal circulation on Mars (Lyons et al., 2005), which is reasonable, as permeability tends to decrease with depth (Manning & Ingebritsen, 1999) and our domain (200-500 m) is much shallower than the domain considered there (∼10 km).We assumed a fracture porosity (ϕ f ) of 100% (i.e., open fractures); we calculated fracture permeability (k f ) as k f = b 2 /12 = 8.3 × 10 −8 m 2 assuming a fracture aperture (b) of 1 mm for all fractures in the domain.Rover photographs of bedrock fractures often show fracture apertures in the range of 1-2 cm (Figures S12 and S13 in Supporting Information S1).However, these photographs are nearly always of fractures expressed at the planet's surface, where they are potentially exposed to freeze-thaw cycles and dehydration of the surrounding rocks, which will cause the fracture apertures to expand.These processes are not as active below the surface, so fracture apertures at depth will be comparatively narrower.Furthermore, at least in the shallow subsurface, fractures tend to be somewhat infilled by dust and/or unconsolidated material (Figure S12 in Supporting Information S1) such that the effective permeability of the fracture is less than that predicted by the cubic law ( = 2 12 , where k f is fracture permeability [m 2 ]).These factors combined with the fact that lithostatic pressure, a force that tends to close fractures, increases with depth, lead us to prescribe uniform 1 mm fracture apertures as an approximate value for Mars' subsurface.

Numerical Mesh and Fracture Generation Algorithm
We generated the fracture networks in our scenarios to be somewhat representative of Mars' subsurface.Because the subsurface on Mars is so poorly characterized, we estimate the fracture density (i.e., the ratio of fracture volume to bulk rock volume) based on rover photographs depicting surface expression of fracture networks at Gale crater (Figure S13 in Supporting Information S1) and extrapolated their distribution into the subsurface.To address the likelihood of variable subsurface architecture, we consider the following range of fracture densities: 0% (unfractured), 0.001%, 0.005%, 0.01%, 0.02%, and 0.035%, shown in Figure 3.
The model is set up in FEHM as a two-dimensional planar domain 50 m wide and with variable domain depth.For scenarios with methane source depth ≤200 m, we use a mesh with domain depth 200 m.For the scenario with source depth 500 m, we use a mesh of depth 500 m.The computational mesh was generated using the Los Alamos National Laboratory (LANL) developed software GRIDDER (https://github.com/lanl/gridder,2018).Mesh discretization is uniform in the x and y directions such that Δx = Δy = 1 m.We randomly generated orthogonal discrete fractures using the 2-D Lévy-Lee algorithm (Clemo & Smith, 1997), a fractal-based fracture model (Geier et al., 1988) produced by random walk.An orthogonal fracture network is a general case, though it can be a reasonable assumption since in mildly deformed (i.e., less tectonically active) bedded rocks, fractures are commonly oriented nearly vertically, with either two orthogonal azimuths or a single preferred azimuth (National Research Council, 1997).The Lévy-Lee model generates a fracture network with a continuum of scales for both fracture length and spacing between fractures.A more detailed description of the algorithm can be found in Text S6.1 of Supporting Information S1.This mesh was then mapped onto a 3-D grid and extended across the width of the domain in the y direction-a single cell across-since FEHM does not solve true 2-D problems.This mapping essentially embeds the fractures in the rock matrix via upscaling of properties (see Section 2.3.2),allowing transfer of fluids and tracers to occur at the fracture-matrix interface.This mesh was then mapped onto a uniform grid.

Upscaling of Fracture Properties
Fractures in our model domain are embedded in the rock matrix via upscaling of permeability and porosity.Fracture permeability k f is upscaled using: where b is the assumed fracture aperture (m) and Δx is the grid/cell block size (m).Upscaled to the grid dimensions of the numerical mesh, the modeled (effective) fracture permeability was 8.3 × 10 −11 m 2 .We upscale fracture porosity using a flow-weighted scheme (Birdsell et al., 2015): giving a model (effective) fracture porosity of 0.001 (0.1%) at the scale of the computational grid (Δx = Δy = Δz = 1 m).The upscaled relationships (Equations 11 and 12) consistently allow the simulation of the governing equations (Equations 2-7) for fractures and matrix using a porous media simulator such as FEHM.This approach is widely used for simulation of flow and transport in fractured rock (Chaudhuri et al., 2013;Fu et al., 2016;Haagenson & Rajaram, 2021;Pandey & Rajaram, 2016).

Atmospheric Column Mixing Model
Methane vented from the subsurface of Mars mixes within the lower atmosphere, where it can be collected as an atmospheric sample by the TLS-SAM instrument.We simulate atmospheric mixing of methane using a one-dimensional, vertical column diffusive transport finite-difference model in order to make general observations about how the instantaneous surface flux translates to atmospheric abundance of methane (Figure 2b).The atmospheric mixing model is sequentially coupled to the subsurface model as a post-processing step.We then use an optimization routine to determine the range of atmospheric transport parameters that minimize the error of calculated methane abundance compared to the TLS-SAM background measurements.This routine is performed for each fracture density case.
We represent the atmospheric mixing using a one-dimensional vertical (z-axis) diffusive transport model (Equation 13).Surface flux from the subsurface transport model is specified as a time varying flux boundary condition in the atmospheric transport model at the ground surface (z = 0 m).The methane diffuses within the atmospheric column, the height of which is equal to the height of the PBL, which varies in thickness hourly and seasonally in 30° increments of solar longitude L s (Guzewich et al., 2017;Newman et al., 2017).
It is important to acknowledge the three-dimensional nature of atmospheric circulation within Gale crater.Methane mixes vertically within the atmospheric column, but it is also advected horizontally and indeed can be transported out of the crater during the day (Pla-García et al., 2019;Viúdez-Moreiras et al., 2021).The complex topography contributes to some variability in PBL depth, and horizontal winds (including slope winds) can sweep through the crater to great effect on the air composition near the surface (Pla-García et al., 2019;Viúdez-Moreiras, 2021).Representing such effects is outside the scope of the present analysis, but the general mixing conditions controlled by PBL variations and bulk atmospheric diffusivity can be reasonably captured using our 1-D atmospheric mixing model.
At night, the PBL height is largely suppressed (<300 m), is approximately constant in height (Guzewich et al., 2017;Newman et al., 2017), and experiences relatively quiescent conditions when compared to the daytime conditions.Despite the suppressed PBL height, the nighttime atmosphere at Gale crater is still turbulent, though to a much lesser degree than during daytime.As the ground surface and atmosphere heats up during the day, the PBL rapidly expands to heights of several kilometers and undergoes a much greater amount of vertical mixing.
In our atmospheric mixing model, we therefore conceptualize the PBL at Gale crater as belonging in either one of two states: "collapsed" or "expanded," each having its own set of atmospheric mixing parameters (Figure S10a in Supporting Information S1).In this way, our approach is conceptually similar to the non-local mixing scheme formulated by Holtslag and Boville (1993), which is implemented in the GEOS-Chem model (GEOS-Chem, 2023; Lin & McElroy, 2010).We note that Holtslag and Boville (1993) used a scheme for vertical diffusion in the PBL to be compared with a GCM based on local vertical gradients, also incorporating non-local (vertical) transport effects for heat and moisture.The governing equations are as follows: where C is the atmospheric methane concentration [kg m −3 ], t is time [s], D c,e is the turbulent/eddy diffusion coefficient [m 2 s −1 ] with the subscript representing a PBL state of either c (collapsed) or e (expanded), z is the vertical coordinate [m], k c,e is a first-order loss rate coefficient [s −1 ].The eddy diffusion coefficients are constant with respect to height, and are therefore independent from changes in dynamic stability conditions.The PBL state is defined as collapsed when h PBL < h thresh , and expanded when h PBL ≥ h thresh , where h PBL is the height of the PBL, and h thresh is the threshold PBL height [m] marking the transition between collapsed and expanded states (chosen to be 300 m).Inspection of the data (Guzewich et al., 2017;Newman et al., 2017) revealed that the maximum PBL height at night was near, but did not exceed, 300 m and exhibited extremely small slopes (i.e., change in height over time).Furthermore, every occurrence of the PBL crossing the 300 m threshold represented the approximate time of occurrence of a change of atmospheric state consistent with our definition of "expanded" or "collapsed" (i.e., the beginning of rapid growth or collapse of the PBL).Because the use of "expanded" or "collapsed" are already simplified designations representing the bulk conditions of a given state, we found the 300 m PBL threshold to be sufficient for our analysis.The loss rate parameter k c,e in this case implicitly combines the effects of photochemical loss as a baseline (assuming a lifetime of methane in Mars' atmosphere of ∼300 years (Atreya et al., 2007;Krasnopolsky et al., 2004)) and any currently unknown destruction mechanisms that may act more rapidly.This loss rate parameter is conceptually identical to the reciprocal of the effective atmospheric dissipation timescale term used in the atmospheric mixing model described by Moores, Gough, et al. (2019).
The diffusive transport equation is solved numerically using a backward Euler finite-difference method scheme, which is implicit in time.The domain is discretized spatially such that Δz = 1 m, and discretized temporally such that each time step Δt = 0.04 sols.For comparison with TLS-SAM methane abundance measurements, modeled abundances are calculated everywhere and recorded at a height of z = 1 m above ground surface to represent the concentration at the height of the TLS-SAM inlet (Mahaffy et al., 2012).

10.1029/2023JE008043
13 of 26 Computation of the transient concentration profiles is complicated slightly by the fact that the model dimensions vary in time via PBL expansion/contraction.At each time step, we modify the number of nodes based on h PBL (t).The methane concentration profile C(z) at the previous time step is translated to the current time step as an initial condition by compressing/extending the profile in proportion to the change in column height such that mass is conserved.For example, when the model domain expands, the vertical concentration profile likewise expands, causing the maximum concentration to be reduced since the profile is spread over a larger area with mass conserved (Figure S10b in Supporting Information S1).This expansion and contraction of C(z) during PBL state transitions can be conceptualized as vertical advection of the tracer within the atmospheric column induced by PBL extension and collapse.
Independent of the state of the PBL (collapsed/expanded), the specified flux boundary conditions are as follow: where j(t) is the time-varying surface mass flux emitted [kg m −2 s −1 ] from the subsurface transport model, and the subscripts represent either indicate collapsed (c) or expanded (e) PBL states.
Atmospheric mixing simulations were run with a spin-up period of 3 MY in order to reach a cyclical steady-state with regard to atmospheric methane abundance.Atmospheric mixing was then simulated for 1 MY, with concentrations recorded at the height of the TLS-SAM inlet (z = 1 m) in order to compare to background methane abundances observed by MSL Curiosity (Webster et al., 2021).Simulations were set up within a differential evolution optimization routine to determine the range of atmospheric transport parameter combinations that best match the observed abundances.Error was quantified in terms of the reduced chi-square statistic,   2  (Press et al., 2007).The parameters optimized were the diffusion coefficients for the collapsed and expanded states (D c and D e , respectively), as well as the methane loss terms for the collapsed and expanded states (k c and k e , respectively).Intuitively, we expect that D e ≥ D c since the expanded state of the PBL is characterized by increased heating and turbulent eddies, which will tend to mix atmospheric tracers more rapidly than would conditions in the more stable collapsed state (Lin et al., 2008).Similarly, we also would expect k e ≥ k c , which accounts for the fact that horizontal advection out of the atmospheric column should be greater in the expanded state than in the collapsed state.We therefore constrained the optimization routine such that: where D c and D e have units of [m 2 s −1 ], k c and k e have units of [s −1 ], and k photochemical is the assumed photochemical loss rate of 1/300 years (∼10 −10 s −1 ).The collapsed-state diffusion coefficient D c has a lower bound on the order of magnitude of free-air methane diffusion in Mars' atmosphere.This lower bound is, in fact, rather conservative, as the binary diffusivity of CH 4 -CO 2 at overnight pressures (800 Pa) and temperatures (180 K) at Gale crater (G.M. Martínez et al., 2017) is approximately 9.4 × 10 −4 m 2 s −1 (Moores, King, et al., 2019).The upper bound is chosen conservatively as double the diffusion coefficient required for methane to fully mix across the depth of the PBL (h PBL ≈ 300 m when in a collapsed state) in 1 hr, which we presume to be the shortest reasonable length of time this condition could be reached.Diffusivity in the expanded state (D e ) is assumed to always be greater than or equal to D c , with an implied maximum value of 10 4 m 2 s −1 .This is a conservative upper bounds considering the estimated eddy diffusivity at higher altitudes in Mars' atmosphere (30-100 km), which are of order 2 × 10 3 m 2 s −1 (Rodrigo et al., 1990) and likely greater than the average diffusivity in the lower atmosphere.

Non-Uniqueness of the Solution
The lack of high-frequency methane abundance data means that this problem is rather poorly constrained.In the analysis described above, we arrive at an optimal solution that minimizes error of the simulated abundances 10.1029/2023JE008043 14 of 26 compared to the sparsely collected observations by modifying four atmospheric transport variables: D c , D e , k c , and k e .The magnitude of the eddy diffusion coefficient (D c,e ) controls how rapidly methane released from the ground surface will mix upwards across the atmospheric column, thereby diluting itself.One can intuit that for the fluxes produced in each subsurface fracture density case, there might be a range of combinations of parameter values that would produce similar annual/seasonal atmospheric abundance patterns, but that would look quite different at the diurnal time scale.We attempt to address this non-uniqueness below in order to provide a more holistic view of the potential diurnal methane abundance patterns dependent on atmospheric mixing rates.
For the fractured subsurface cases that produce the best overall fit to the observed methane abundances in the differential evolution algorithm, we analyze the surrounding parameter spaces that produce similar results with regard to overall reduced   2  value.The reduced   2  statistic is used extensively in goodness of fit testing, and has been applied previously by Moores, Gough, et al. (2019) and Webster et al. (2018b) for comparing modeled methane abundance to TLS-SAM measurements (see Press et al. (2007) (Bevington, 1969).
The "best" fit in each fracture density case is characterized by the minimum   2  .For a given fracture density case, we subset the simulation outcomes to the parameter combinations with error in the range:   2  ≤ ( min  2  ) + 0.5 .The 0.5 was arbitrarily chosen to provide a reasonable sample size of candidate solutions, and corresponds to an approximately 8% change in goodness-of-fit probability as calculated by the   2  statistic.Candidate solutions in this range therefore have similar levels of fit to the "best" scenario, and generally sample a wide range of parameter values and combinations.We then divide this parameter space into 4 scenarios: (a) lowest D c , (b) highest D c , (c) smallest k e /k c ratio, and (d) largest k e /k c ratio.A summary of the effective parameter value ranges are given in Table 1, and the actual parameters estimated from the optimization algorithm and used in the various scenarios are detailed in Table 2.The end-member scenarios for diffusivity are conceptually similar to the transport end-members investigated by Moores, King, et al. (2019), in which they considered both a completely static, stably stratified near-surface air layer, in addition to a well-mixed near-surface air layer.

Results and Discussion
We present numerical simulations of transient surface methane flux caused by barometric pressure-pumping into Mars' atmosphere from a constant underground source.We simulated this transport mechanism acting in a range of subsurface fracture network topologies by varying the fracture density in our domain (Figure 3).We then translate methane flux (i.e., surface emissions) into atmospheric abundance (i.e., mixing ratio, in ppbv) by supplying the computed methane fluxes to the atmospheric diffusion model described in Section 2.4.We assess our simulations by comparing their fit to TLS-SAM's observed background methane abundance fluctuations (Webster et al., 2021), which included two non-detections at mid-sol measurements in northern summer.We identify the best-fitting simulations by computing the reduced chi-square statistic for the modeled methane abundance variation over one Mars year (L s 0-360°).Note that the TLS-SAM measurements were taken over multiple Mars years (MY).The parameter optimization approach proceeds based on the overall   2  value (Table 2), which is calculated using all background TLS-SAM measurements.The optimization approach therefore inherently selects scenarios that best match both the seasonal and sub-diurnal variations.However, due to the paucity of measurements taken at different times of day (i.e., those that would be indicative of sub-diurnal methane variations), the optimization approach is more likely to select parameter combinations that more closely match the seasonal variations observed rather than the sub-diurnal variations.To address this, we pick out the fracture density cases that match the seasonality well (Overall   2  in Table 2), and examine the surrounding parameter space to observe changes in sub-diurnal methane variations that were measured in northern summer (Summer   2  in Table 2).We do not explicitly optimize the parameter space to reduce error of sub-diurnal variations in the northern summer period.
Though we investigated a range of methane source depths, because our simulations reach a cyclical steady-state, there was negligible variance in the timing of surface fluxes caused by varying source depth since the subsurface becomes equivalently populated with methane gas.Therefore, the primary source of variance in the timing of surface flux pulses was the fracture density.The best-fitting cases had a fracture density of 0.01% (Figures 4 and 5), followed 10.1029/2023JE008043 15 of 26 closely by cases with fracture density 0.035% (Figures S18 and S20 in Supporting Information S1) and 0.02% (Figures S17 and S19 in Supporting Information S1).The main focus of this paper is on characterizing the timing of methane variations, so the source depth does not matter for the rest of the analysis presented here.The effect of source depth would be more pronounced in the case of a source term that produces methane episodically instead of continuously, such that subsurface concentrations were not at cyclical steady-state.For example, one can assume that the barometric pumping mechanism is essentially always active, and the year-over-year intensity is relatively constant.If the methane source produces methane for a very short time duration (e.g., on the order of a sol), the source depth would have to be quite near to surface to reach the atmosphere in significant quantities.If such a source produces methane very infrequently, then one would also not expect the signal to match the observations, regardless of depth.On the other hand, a longer-duration subsurface source that produces methane more frequently would more closely approximate a continuous source and could be located deeper and could feasibly produce the observed methane variations.
For each fracture density case, the optimization algorithm arrives at a "best" solution in terms of error using some combination of atmospheric transport parameters.However, due to the non-uniqueness of potential solutions  .Therefore, we investigate several atmospheric transport end-members in the candidate parameter space for each of the fracture density cases, the three best of which (fracture density 0.01%, 0.02%, and 0.035%) are presented in Table 2.These scenarios are described in Section 2.4.1, with parameter values detailed in Table 2.It is worth noting that the subsurface cases we investigate with low fracture density (0%, 0.001%, and 0.005%) produce methane abundance patterns that are almost completely out of phase with the observed abundance pattern, regardless of the choice of atmospheric transport parameters.These results are included in Supporting Information S1.
As a general discussion related to evaluating the appropriateness of the modeled diffusivities, atmospheric mixing time is one metric by which we can estimate whether a given set of parameters is realistic.The approximate time  2 and Section 2.4.1.Panel (e) is the "best" fitting scenario (corresponds to top row in Table 2), and panel f is the corresponding surface methane flux.

10.1029/2023JE008043
18 of 26 required for a system to reach a fully mixed state in response to an instantaneous point source located on a boundary (Fischer et al., 1979) is described by: where t ss is the time [s] of full mixing (i.e., when maximum deviation from the steady-state concentration profile is <1%), L is the length of the domain

Seasonal Methane Variation
The best overall fit to TLS-SAM measurements arose in the cases where fracture density was 0.01%.Several features are apparent in the abundance plots (Figures 4a-4e) showing seasonal atmospheric abundance changes on Mars.Note that the gray band apparent in the plot is the result of large diurnal variations in the simulated abundance.The black line represents the night-time average abundance (calculated between 0:00 and 2:00 LMST) for the sake of visualization, since a significant majority of measurements were performed in this window.It should be noted that the error is calculated based on the simulated instantaneous methane abundance values rather than this night-time average.Our analysis assumes that the seasonal variation in atmospheric methane reported in Webster et al. (2021) is real, though we again note that the infrequent measurements pose a problem to identifying a statistically significant seasonal variation, as reported by Gillen et al. (2020).
Generally, the "best" fit scenario (Figure 4e) represents the seasonal methane variations well throughout the Mars year, especially the elevated abundances in northern summer (L s 90-180°) and gradual decline in northern autumn (L s 180-270°).However, exceptions occur in several time periods.The first occasion is from L s 32-70°, marking the approximate middle of northern spring.Over this interval, the simulated values generally overestimate atmospheric abundance.Second, the simulation underpredicts abundance at L s ∼ 216°, in northern autumn.The difference between simulated and observed abundances at this point is less pronounced, as the simulated diurnal abundance (shown in gray) falls very nearly within one standard error of the mean (SEM) for this measurement, as indicated by the error bars on the plot.Third, the simulations also underpredict atmospheric abundance at L s = 331°, the middle of northern winter.
The results composite in Figures 4a-4d shows the effect of the atmospheric transport end-members investigated for fracture density 0.01%.The general character of the seasonal methane abundance variation remains in each scenario, though the details vary somewhat.Scenarios with smaller D c (such as scenarios a and d) have a greater range of diurnal abundance (gray band).Smaller D c in general means that the mixing of methane across the depth of the atmospheric column takes longer.This allows methane concentrations near the emission surface (e.g., at z = 1 m, where the TLS-SAM inlet is located) to build to higher values before subsequent mixing.Scenarios with smaller D c also seem to produce a more pronounced increase in atmospheric methane abundance during northern winter.Scenarios with higher diffusivity (e.g., scenario b) begin to approach an instantaneous mixing condition.Instantaneous mixing may be a reasonable approximation under conditions where the PBL is extremely unstable.On Earth, for example, this can be the case during hot days with significant convection creating thermal instability, or stormy days, during which increased wind speeds near the ground surface generate increased drag and promote the formation of turbulent eddies with vertical mixing.However, under most conditions an instantaneous mixing assumption will tend to overestimate vertical mixing (Lin & McElroy, 2010).We initially used a more simplified instantaneous mixing approach similar to what done in Moores, Gough, et al. (2019), but opted for a diffusive partial mixing model as being more realistic of general atmospheric conditions (discussed in more detail in Text S4 of Supporting Information S1). 10.1029/2023JE008043 19 of 26

Sub-Diurnal Methane Variation
With the goal of determining useful timing of TLS-SAM measurements, we also examined our simulations over shorter time scales, looking at the diurnal variations in methane abundance in northern summer (Figure 5).Northern summer is the only season in which TLS-SAM has performed daytime enrichment method measurements, generally collected around noon (Webster et al., 2021).All other measurements have been collected close to midnight, so this is therefore the only season in which we have clues as to the possible sub-diurnal shape of methane variations.Direct observation of a sub-diurnal shape has not been possible due to instrument operational constraints of TLS-SAM, which cannot make multiple measurements on the same sol.The defining characteristic of these results (Figure 5) is the sharp drop-off in atmospheric abundance that occurs between approximately 8:00 and 16:00 local time (LMST), which coincides with the elevated PBL height seen in the bottom panel of the same figure.Note that we use a 24-hr time convention for the remainder of the discussion, where 0:00 -11:59 LMST represent the morning from midnight to just before noon.In our model, the drop-off in abundance is controlled largely by the mid-day extension of PBL height, and also the generally 2-3 order of magnitude difference between D e and D c (Table 2).When the PBL collapses in the early evening (∼17:00 LMST), it remains relatively shallow (i.e., atmospherically quiescent compared to mid-day) through the night until early the next morning.The atmospheric mixing ratio responds accordingly by rebounding somewhat after the PBL collapse, after which point it holds relatively steady into the following morning.
The "best" scenario in terms of error shown in Figure 5e generally reproduces the observed summer methane abundances.The model slightly underpredicts methane abundance relative to that observed at L s = 158.6°(yellow circle), though the modeled concentration is within one SEM of the measured value.The mid-day non-detections (L s 120.7 and 134°) are generally captured by the model, as well as the positive TLS-SAM detection that was collected between them (L s 126.3° at 23:56 LMST).The latter point distinguishes this case from the higher-fracture-density cases (0.035% and 0.02%), which where not able to match this intermediate observation regardless of the scenario considered (Figures S19 and S20 in Supporting Information S1).An accurate match to the observed abundances is thus controlled by both the assumed subsurface fracture density and the parameters in the atmospheric transport model.
For the case shown in Figure 5f (i.e., the "best" case in terms of error), elevated daytime fluxes have a somewhat bimodal pattern (i.e., two primary methane flux pulses).The first occurs between 4:00 and 6:00 LMST, and has substantially greater magnitude (by a factor of 5-11) for the dates with non-detections (L s = 120.7,134°) and at L s 158.6° than it does on the dates of the other measurements.The second primary methane pulse occurs between 15:30 and 17:00 for L s = 103.4,126.3, and 142.4°, and less strongly (by a factor of 1.4-5) between 16:00 and 18:00 for the L s = 120.7,134° (non-detects) and L s = 158.6°.The timing of the surface flux pulses is dictated entirely by the subsurface fracture density; that is, the fracture topology.The surface flux pulses are produced in response to the small morning barometric pressure drop occurring at approximately 3:00, and the large mid-day pressure drop occurring between 7:40 and 16:00.If the subsurface were a homogeneous medium, we would expect a surface flux pulse roughly coincident with the pressure drop, having a Gaussian shape in time.This is actually observed in our model as fracture density increases: for example, in the case where fracture density = 0.035%, the surface flux has fewer individual spikes, and is characterized by a more "diffuse" flux pattern with center-of-mass near the middle of the large mid-day pressure drop (Figure S20f in Supporting Information S1).The sparse fracture network in the present case (fracture density 0.01%) does not release methane at the surface in sync with the pressure drops-trace gases must work their way tortuously through individual fractures.The surface pressure wave propagates through the fractures and is attenuated by the rock matrix, leading to varying degrees of phase lag in the subsurface signal.Over multiple barometric pressure cycles, methane gas is brought closer to the surface through different fracture pathways-the variety of travel pathways leads to different surface breakthrough times depending on the pressure propagation and gas transport history within each fracture.This helps explain why the individual flux pulses shown in this case vary so much in magnitude despite being forced by relatively similar atmospheric pressures.
Examination of the end-member scenarios reveals some key differences imbued by the choice of atmospheric transport variables (Figures 5a-5d).In terms of   2  , there is little to distinguish the end-member scenarios examined, although scenario c clearly performed worse than the rest over this time frame.Scenarios a and d used small values of D c (of order ≤0.01 m 2 s −1 , which is on the order of magnitude implied by a 1-sol crater mixing time, and two orders of magnitude greater than binary CH 4 -CO 2 diffusion), the effect of which is apparent in the rapid ORTIZ AL.
10.1029/2023JE008043 20 of 26 spike in methane abundance between 4:00 and 7:00 LMST.This spike is a direct result of the methane surface flux pulses occurring between 4:00 and 6:00 LMST; the smaller values of D c cause the sensor at z = 1 m to more readily feel the effects of these pulses before they eventually mix by diffusion into the rest of the atmospheric column.The effect of these early morning methane pulses is greatly muted in scenarios b and c, which had much greater values for these mixing coefficients (of order ≥6 m 2 s −1 ).We further interrogated the candidate solution parameter space generated by the differential optimization algorithm in order to understand the interaction between atmospheric mixing parameters, with results in Text S7.4 of Supporting Information S1.Diffusion coefficients D c and D e , unsurprisingly, are positively correlated such that smaller D c corresponds to a smaller D e .The candidate solution space contains diffusion coefficient values such that range of the ratio D e /D c is between 59 and 678 (Figure S22 in Supporting Information S1), with a mean value of 351.We provided bounds to the algorithm for this ratio in 1 ≤ D e /D c ≤ 1,000, so the atmospheric mixing model evidently favors comparatively large daytime eddy diffusivities compared to those during the collapsed state, although the absolute magnitudes of these diffusivities do not overly affect the results in terms of error.A linear regression on D e = f(D c ) yields a slope of 10.8, with an adjusted R 2 value of 0.85.Also unsurprisingly, first-order methane loss rate parameters k c and k e are inversely correlated in order to preserve mass balance in time.The range of the ratio k e /k c is 1.01-3.21(Table 2) having mean value 1.46, with the overall best scenarios in terms of error coming out of ratios close to unity.A linear regression on k e = f(k c ) yields a slope of −1.1, with an adjusted R 2 value of 0.67.

Effects of Dust Devil Pressure Drops on Flux Timing
As part of making predictions about timing of atmospheric methane measurements, we also considered the effects of dust devil vortices on surface flux of methane in the vicinity of the rover.We considered this because MSL Curiosity is currently climbing Aeolis Mons (a.k.a.Mt.Sharp), and will be doing so for the remainder of the mission.Observational data and Mars Weather Research and Forecasting (MarsWRF) General Circulation Model (Richardson et al., 2007) simulations of Gale crater indicate a gradual increase in vortex detections during most seasons as the MSL Curiosity rover ascends the slopes of Aeolis Mons (Newman et al., 2019;Ordóñez-Etxeberria et al., 2020).The primary reason for this is related to the increase in topographic elevation, which encourages vortex formation because of the greater difference in daytime surface-to-air temperatures (Newman et al., 2019).More discussion on this is provided in Text S5 of Supporting Information S1.
We describe these dust devil simulations in Text S5 of Supporting Information S1.We considered pressure drops associated with dust devils over a range of duration and intensity.As expected, the greatest surface flux is caused by dust devils with the longest duration (25 s) and largest pressure drop (5 Pa; Figure S11 in Supporting Information S1).However, the total mass of methane emitted in this scenario was 9.4 × 10 −10 g, which has a negligible effect on atmospheric methane abundance in our model.Overall, dust devils likely do not make much of a difference in surface methane emissions.This makes sense, as the diurnal pressure variations by comparison have magnitude of order several tens of Pa, with the primary pressure drop occurring over an interval of several hours.We can therefore likely ignore the effects of dust devils on overall timing of methane variations, which is encouraging since we are unable to predict the occurrence of individual vortices.

Implications for Future Measurements
Confirming and characterizing the apparent diurnal variability of methane has been highlighted by the TLS-SAM team as the next key step to understanding methane abundance and circulation at Gale crater.At the time of writing, Mars' summer period approaches, the timing of which is coincident with prior measurements that suggested subdiurnal methane variations (L s 120-140°).This makes northern summer a prime candidate for potential corroboration of the hypothesized subdiurnal methane variations.The TLS-SAM wide range pumps have performed exceptionally well, and have already exceeded their flight lifetime requirements, but we need to be prudent in planning their use in future measurements.This compels the need to choose strategic sampling times in order to learn as much as possible about methane seepage and circulation patterns at Gale crater.Strategic atmospheric sampling using TLS-SAM during this northern summer time frame has the potential to validate and contextualize the results of our coupled subsurface-atmospheric mixing model as well as the previous measurements suggesting diurnal methane variations.
With the goal of more robustly characterizing diurnal methane variability, we would propose a set of enrichment runs in the period L s 120-140°, which occurs September-October 2023.In the interest of conserving TLS-SAM pump life, we propose initially performing a minimum of two measurements.The first proposed measurement would establish a baseline for the second in addition to providing comparison to measurements conducted in previous MYs, while the second measurement would aim to extend the current characterization of diurnal methane variability.The measurements we propose would correspond to the approximate time of year of the previous two mid-sol samples, as well as the apparent generally elevated methane abundance occurring in northern summer.Ideally, the samples would also be coordinated such that they coincide with TGO solar occultations on any of either 25 September, 27 September, 9 October, or 11 October 2023 for potential cross-comparison of measurements.Both enrichment runs should be performed identically to each other with the exception of local time conducted.A version of the dual-enrichment run modified slightly from the procedure of previous measurements (Webster et al., 2018a) would provide better quantification of background methane and better conserve pump life without deviating significantly from previous run procedures (see Text S3 in Supporting Information S1 for a more complete description of the modified procedure).
The first sample we propose should ideally be performed around L s 126° to coincide with time-of-year of the previous MY positive detection on sol 2626, which was conducted between the two daytime non-detections in 2019 (Webster et al., 2021).This would serve as a baseline observation, both for the sake of comparison to the following measurement, as well as to the previously established baseline abundance for this period.Performing the measurement within the 23:00-3:00 LMST time range would make this measurement immediately comparable to most measurements from previous MYs, and additionally would refresh the baseline for the current MY and second run.
The second measurement would ideally be collected at a previously unmeasured time, and would be chosen to provide new insight into the methane emission and mixing mechanisms at play, in addition to extending the characterization of the apparent diurnal variability.We envision two primary candidate timing windows for this proposed measurement, which we hereafter refer to as I and II.Window I would take place between 6:30 and 10:00 LMST with the goal of further constraining the drop in observed methane abundance that seems to occur between midnight (0:00 LMST) and 11:20 LMST.Prior work using atmospheric transport models (Figure 8 in Moores, King, et al. (2019), Pla-García et al. (2019), and Viúdez-Moreiras (2021)), in addition to the present work, predict that this drop occurs some time mid-/late-morning due to the upward extension of the PBL column and reversal of horizontal flows from downslope flows-converging during nighttime toward the bottom of the crater-to upslope flows-diverging during daytime through the crater rims and Aeolis Mons.A measurement in Window I would further constrain the timing of the apparent drop in methane abundance; for instance, elevated methane levels late in this window would aid the argument that PBL extension and the accompanying transition to divergent flows are strongly linked to the daytime drop in abundance.Methane abundance noticeably higher than the baseline measurement near midnight would imply additional flux in the intervening morning hours based on our model.However, if the magnitude of the difference is not overly large, it could be difficult to parse out the effects of a morning flux pulse (e.g., Figures 5a and 5d), gradual overnight methane accumulation, or simply sol-to-sol abundance variation.
Window II encompasses the time between 18:00-21:00 LMST, and a sample therein would serve to characterize the hypothesized rise in methane levels at sunset, post-PBL collapse (∼17:00).A measurement early in this window (18:00-19:00) could provide useful information regarding potential surface release mechanisms.If methane builds up rapidly to concentrations consistent with or above nighttime values, it could be indicative of daytime methane emissions, such as those caused by barometric pumping, though not exclusively due to this mechanism.Along that methane abundance noticeably greater than nighttime values (e.g., Figures S19a and  S19d in Supporting Information S1) would suggest either the occurrence of mid-/late-afternoon flux pulses, or that the magnitude of nighttime emissions is less than that estimated in other studies (or is nonexistent), both of which would also be consistent with barometric pumping.Abundances lower than observed nighttime values, on the other hand, could suggest gradual evening/overnight methane accumulation, which may point to an emission mechanism other than barometric pumping, which produces primarily daytime fluxes.

Conclusions
This study investigates the transport of subsurface methane in fractured rock into Mars' atmosphere driven by barometric pressure fluctuations at Gale crater.The subsurface seepage model is coupled with a 1-D atmospheric mixing model in order to simulate atmospheric concentrations within an evolving PBL column in response to transient surface emissions and compares them to MSL abundance measurements.Though our 1-D atmospheric mixing model does not fully represent the complex 3-D atmospheric circulation conditions (e.g., horizontal exchange, winds, and slope flow) that exist at Gale crater, the bulk mixing times estimated by Pla-García et al. (2019) are reasonably accounted for.Atmospheric transport variables are chosen by an optimization routine such that they minimize the error compared to MSL TLS-SAM measurements, which include seasonal and sub-diurnal abundance variations.The simulations are evaluated based on how well they represented seasonal and diurnal variations in atmospheric methane concentrations, including daytime non-detections observed by MSL TLS-SAM.Part of the investigation involves simulating subsurface transport in rocks covering a range of fracture densities.To that end, a lower bound on subsurface fracture density of 0.01% is established, below which the seasonal atmospheric variations driven by barometric pumping are out-of-phase with observations.We examine the sub-diurnal atmospheric methane variations produced by our simulations in Mars' northern summer, a time period chosen due to its coincidence with previous measurements suggesting the presence of large diurnal abundance fluctuations.Several key features were identified in the best-performing simulations.Simulations indicated a pre-dawn methane surface flux pulse (4:00-6:00 LMST) that may be detectable before PBL thickness increases and upslope (divergent) circulation develops.Detection of a large methane spike would be suggestive of barometric pumping, and would add to the evidence supporting a localized emission source in the interior of Gale crater, such as the highly fractured Murray outcrops as mentioned in Viúdez-Moreiras et al. (2021) or the crater floor mentioned in Pla-García et al. (2019).Another feature identified was a large abundance depression during mid-sol between 11:00 and 17:00 coincident with PBL extension and divergent slope flows, followed by a rapid rebound in methane abundance following PBL collapse in the early evening and development of downslope flows converging toward the bottom of the crater.As a way to test our proposed transport mechanism and extend the current characterization of diurnal methane variation, we propose a set of two TLS-SAM enrichment measurements for the middle of Mars' northern summer (L s = 120-140°), with the option of either a mid-/late-morning or an early evening measurement.Each measurement has high potential to better-constrain the current understanding of the timing of either the apparent morning drop in methane or evolution of nighttime methane increase, respectively, and both measurements have modest potential to incrementally support or refute the influence of a barometric pumping mechanism on diurnal methane variations at Gale crater.
The modeled methane abundances presented in this work are controlled by two factors: the subsurface transport pattern driven by barometric pumping and the PBL dynamics.Though driven by the same barometric signal, surface methane flux patterns in our model varied significantly with subsurface fracture network topology (i.e., fracture density).Fracture density controls the degree to which the atmospheric pressure signal propagates into the subsurface, both in terms of overall depth and phase response.So important is the communication of the atmospheric pressures with the subsurface that cases we considered with very low fracture density (≤0.005%) produced surface flux and abundance patterns that were almost completely out of phase with TLS-SAM observations.In our coupled 1-D atmospheric mixing model, we chose a handful of atmospheric transport parameters to approximately describe the PBL mixing dynamics, which essentially controlled the rate at which mixing from the surface methane emission would occur in the atmospheric column at different times of day.The atmospheric methane abundance was highly sensitive to these parameters, which exerted a great influence on both the seasonal and sub-diurnal abundance patterns.Despite this, our sensitivity analysis showed that no combination of atmospheric transport parameters in our model could generate abundances that were in-phase with the observed patterns for the low fracture density cases (≤0.005%).This implies an important interplay between the influence

Figure 1 .
Figure 1.Schematic of the barometric pumping mechanism, which has ratcheting enhanced gas transport due to temporary immobile storage.The upward advance of the gas during barometric lows is not completely reversed during subsequent barometric highs due to temporary storage of gas tracer into rock matrix via diffusion.Adapted from Figure 1 in Harp et al. (2018).
the 2-D gradient operator (operating in the fracture plane), ρ is the air density [kg m −3 ], t is time [s],  ⃗  is the in-plane aperture-integrated fracture flux [m 2 s −1 ],    is the volumetric flux [m 3 /(m 2 s)] of air in the rock matrix,    denotes the normal at the fracture-matrix interfaces pointing out of the fracture (I), b is the fracture aperture [m], μ is the dynamic viscosity of air [Pa s], P f is air pressure within the fracture [Pa], k f is fracture permeability [m 2 ], g is gravitational acceleration [m s −2 ], and z is elevation [m].The right-hand side of Equation 2 represents the fluxes across the fracture-matrix interface, where positive    ⋅ f and C m are tracer concentrations [mol  kg −1  ] in the fracture and matrix, respectively; D is the molecular diffusion coefficient of the tracer [m 2 s −1 ];    is the normal at the fracture-matrix interfaces pointing out of the fracture (I); and    is the tracer source in the fracture plane [mol m −2 s −1 ].The first term on the right-hand side of Equation 6 represents the tracer mass fluxes across the fracture-matrix interfaces.Note that the mass fluxes across fracture-matrix interfaces include advective and diffusive fluxes.Even in the absence of significant air flow in the matrix, diffusive flux exchanges between the fracture and matrix persist and are included in our formulation.

Figure 2 .
Figure 2. Schematics of model domains used in flow and transport simulations.(a) The subsurface fracture-rock flow and transport model.Fracture network generated using the Lévy-Lee algorithm.Fractures are shown in red, with rock matrix in blue.A methane source located in the methane production zone produces methane at a constant rate.(b) Schematic of the coupled subsurface-atmospheric mixing model.Methane is emitted into the atmosphere from the subsurface fractured-rock transport model.Mixing of methane occurs via 1-D vertical diffusion within the atmospheric column (light blue region), the volume of which varies seasonally and hourly based on the evolution of the planetary boundary layer height, h PBL (t).The atmospheric mixing model is described in detail in Section 2.4.

Figure 3 .
Figure 3. Schematic of the subsurface model domain showing subsurface fracture densities used in this study.

Figure 4 .
Figure 4. Composite of atmospheric mixing end-member scenarios simulating atmospheric methane abundance for the case with fracture density 0.010% showing seasonal methane variation.Panels (a-e) compare simulated (stars, lines) to measured (circles) atmospheric methane abundance values plotted against solar longitude, L s [°].Night-time averages of the simulated abundance (thick black line) are plotted to aid visualization because of the large diurnal variations present (gray band).Measured abundances are from Webster et al. (2021).Note that some measurements were collected in different Mars years.Panel letters a-d correspond to lettering of atmospheric transport parameter end-member scenarios described in Table2and Section 2.4.1.Panel (e) is the "best" fitting scenario (corresponds to top row in Table2), and panel f is the corresponding surface methane flux.

Figure 5 .
Figure 5. Composite of atmospheric mixing end-member scenarios simulating atmospheric methane abundance for the case with fracture density 0.010%.Panels (a-e) compare simulated (stars, lines) to measured (circles) atmospheric abundance values in local time, LMST, for northern summer, which highlights the day-night difference in abundance largely caused by the elevated planetary boundary layer (PBL) height h PBL .Simulated abundances of the sols with non-detections are indicated by dashed lines.Measured abundances are from Webster et al. (2021).Note that all measurements were taken on different sols and, in some cases, different Mars years, with the solar longitude, L s [°] of the measurement indicated on the plot by its color.Panel letters (a-d) correspond to lettering of end-member scenarios described in Table 2 and Section 2.4.1.Panel (e) is the "best" fitting scenario in terms of model error (corresponds to the top row of Table 2), and panel (f) is the corresponding surface methane flux.Surface flux in local time (solid and dashed lines as above) is plotted along with PBL height (dotted line).Atmospheric pressure (blue line) is plotted without visible scale, but the minimum and maximum values shown are approximately 703 and 781 Pa, respectively.The pressure time series shown is from L s = 120.7°;pressures on the dates of the other measurements are different but similar in shape.Taking into account the derived crater mixing times (t ss ) calculated from D c and D e and comparing to that estimated by Pla-García et al. (2019) (t ss = 1 sol), suggests that scenarios a and d are likely to be more closely representative of actual conditions with regard to mixing timescales.
[m], and D is the diffusion coefficient [m 2 s −1 ].Three-dimensional atmospheric modeling performed by Pla-García et al. (2019) determined that the mixing time scale for martian air within Gale crater is approximately 1 sol.Applied to the present model, this implies a collapsed-state diffusion coefficient D c ≈ 0.4 m 2 s −1 (where L ≈ 250 m), a minimum expanded-state value of D e = 25.2 m 2 s −1 occurring at L s = 130° (where max L = 2,045 m), and a maximum expanded-state value of D e = 219 m 2 s −1 (where max L = 6,017 m).The implied value of D c calculated above additionally is of the same order of magnitude as the eddy diffusion coefficient at z = 1.3 m estimated by G.Martínez et al. (2009).We therefore give preference in the discussion to parameter-space solutions in our mixing model that have diffusivities of similar orders of magnitude (0.1 ≤ D c ≤ 1.0 m 2 s −1 and 25 ≤ D e ≤ 500 m 2 s −1 ).
Considering these simulations in terms of crater mixing time (t ss ) of ∼1 sol estimated by Pla-García et al. (2019) also favors the scenarios with smaller D c .For an approximate collapsed-state PBL height of 250 m, mixing times for Table 2 scenarios are as follows: (best) 0.05 sols, (a) 4.3 sols, (b) 0.04 sols, (c) 0.07 sols, and (d) 0.75 sols.However, the collapsed state only accounts for part of each sol.The maximum diurnal PBL height during the expanded state varies from 2,045 to 6,017 m throughout the Mars year.For max h PBL = 2,045 m-which occurs in northern summer-the inferred mixing time t ss is (best) 0.01 sols, (a) 0.8 sols, (b) 0.004 sols, (c) 0.14 sols, and (d) 0.28 sols.For max h PBL = 6,017 m-which occurs during northern winter-the inferred mixing time t ss is (best) 0.07 sols, (a) 6.56 sols, (b) 0.04 sols, (c) 1.18 sols, and (d) 2.4 sols.Scenarios a and d most closely approximate the presumed crater mixing time, though it should be noted that there can be significant variation in mixing times throughout the Mars year (Pla-García et al., 2019; Yoshida et al., 2022), and our atmospheric mixing model is not set up to account for these variations due to representing D e with a single value.
Note.Effective minimum and maximum values for the range of values for D e and k e are derived from the ratios prescribed to the parameter estimation algorithm, described in Equations 17 and 19, respectively.Values and Considered Ranges of the Model Parameters Used in the Different Aspects of the Numerical Modeling Webster et al., 2021) TLS-SAM abundance values, modeled abundance values, and the standard error of mean (SEM) uncertainties of the TLS-SAM data (Table2inWebster et al., 2021).A value of Note.D c and D e are in units of [m 2 s −1 ], and k c and k e are in units of [s −1 ].Scenarios are described as follows according to the parameter space discussed in Section 2.4.1:(best) parameters with overall best fit to TLS-SAM data, (a) lowest D c , (b) highest D c , (c) smallest k e /k c ratio, and (d) largest k e /k c ratio.

Table 2
Description of Parameters Used in Various Atmospheric Mixing Scenarios for the Three Best-Performing Fracture Densities