Introduction
The Martian atmosphere was established four billion years ago and most of the carbon in the Martian atmosphere has subsequently escaped over time (Mahaffy et al. Reference Mahaffy, Webster, Atreya, Franz, Wong, Conrad, Harpold, Jones, Leshin, Manning, Owen, Pepin, Squyres and Trainer2013; Webster et al. Reference Webster, Mahaffy, Fkesch, Niles, Jones, Leshin, Atreya, Stern, Christensen, Owen, Franz, Pepin and Steele2013). The existence of methane emissions into the atmosphere of Mars is an enticing possibility worthy of exploration in situ. The detection of atmospheric methane on Mars may indicate biogenic origin (Chastain & Chevrier Reference Chastain and Chevrier2007) or not (Parnell et al. Reference Parnell, Boyce and Blamey2010). There were three independent measurements of Martian atmospheric methane – one from the Planetary Fourier Spectrometer onboard Mars Express (Formisano et al. Reference Formisano, Atreya, Encrenaz, Ignatiev and Giuranna2004) (including at Elysium Planitia), another from the Fourier Transform Spectrometer at the Canada–France–Hawaii Telescope, Mauna Kea, Hawaii (Krasnopolsky et al. Reference Krasnopolsky, Maillard and Owen2004), and finally, simultaneously from three ground-based telescopes with high-dispersion infrared spectrometers (Cryogenic Near-IR Facility Spectrograph on the NASA Infrared Telescope Facility, Mauna Kea, Hawaii; Near-IR Spectrograph on the Keck-2 Observatory, Mauna Kea, Hawaii; and the Phoenix Near-IR spectrometer on the Gemini South telescope, Cerro Pachon, Chile) (Mumma et al. Reference Mumma, Villanueva, Novak, Hewagama, Bonev, DiSanti, Mandell and Smith2009). The methane measurements imply that an active methane source exists on Mars and that some sort of destructive process is taking place. The methane has been tentatively detected with an approximate concentration of 10 ppbv (but with large uncertainties of ±5 ppb) requiring a production rate of ~1 × 105 molecules cm−2 s−1 or 126 tonne y−1 averaged over the Martian surface. The methane appears to be localized, however, as extended plumes of ~30 ppbv average concentrations with a maximum plume size of ~170 km from discrete sources releasing around 0.6 kg s−1, perhaps resembling hydrocarbon seeps (though again, there is much uncertainty). Similar measurements were incorporated into a global circulation model (GCM) to predict from where Martian methane plumes may have originated (Mischna et al. Reference Mischna, Allen, Richardson, Newman and Toigo2011). More recently, measured peaks of methane are lower at 8 ppb. Indeed, the empirical observations have been questioned on the basis of Earth's atmospheric distortion (Zahnle et al. Reference Zahnle, Freedman and Catling2011). The measured seasonal variability of 6 months does not agree with the 300-year long chemical lifetime of methane expected in the Martian atmosphere (Lef`evre & Forget Reference Lef`evre and Forget2009). Martian dust storms could act as possible methane sinks reducing their atmospheric lifetime (Farrell et al. Reference Farrell, Delory and Atreya2006). Methane clathrates are a possible source of the methane on Mars (Chastain & Chevrier Reference Chastain and Chevrier2007). Methane emissions correlate with water vapour measurements suggestive of metastable clathrates while explaining the observed 6 monthly timescale of methane production/losses (Chassefi`ere Reference Chassefi`ere2009). Methane, of course, is a potential biomarker since on Earth, methane is 90–95% biogenic in origin. The existence of deeply buried microbial life 2–3 km below the Earth's surface suggests a potential analogous biotic origin to methane on Mars from extant or extinct microbial biota buried in its subsurface. However, the most favoured hypothesis is that Martian methane has a geological origin (serpentinization) given the lack of organics on the surface (Atreya et al. Reference Atreya, Mahaffy and Wong2007). Serpentine is formed from the aqueous alteration of olivine in the presence of CO2 at modest temperature and is expected to be widespread on Mars (Parnell et al. Reference Parnell, Boyce and Blamey2010). Serpentinization by water releases hydrogen as a potential energy source for life catalysed by Fe-Ni (Munterer Reference Munterer2010). Other possibilities than serpentinization appear implausible. Volcanic emissions of methane would require recent Martian volcanism of which there is currently no evidence and asteroid/comet/interplanetary dust particle delivery would require recent delivery for which again there is no evidence, while any subsurface methanogenic population must be marginal due to the paucity of subsurface energy sources today. Nevertheless, biogenic origin remains a possibility. Whether of geological or biological origin, the methane may be trapped as methane clathrates in water ice to be released under reduced pressure or higher temperature. In situ measurements of methane by the Mars rover Curiosity have yielded inconclusive results (Webster & Mahaffy Reference Webster and Mahaffy2011).
Our goal in this paper is to investigate how Mars rovers might be deployed to make in situ measurements of methane to confirm its presence and, further, to exploit that information to determine the nature of methane plumes and their origin. The 30 kg Kapvik micro-rover prototype was deployed in an analogue mission at asbestos mine sites to simulate teleoperations for searching and sensing methane on Mars (Qadi et al. Reference Qai, Cloutis, Samson, Whyte, Ellery, Bell, Berard, Bovin, Haddad, Lavoie, Jamroz, Kruzelecky, Mack, Olsen, Perrot, Popa, Rhind, Sharma, Stromberg, Strong, Tremblay, Wilhelm, Wing and Wong2015). A Mars rover may autonomously search for methane sources, perhaps with a scout rover, to increase the science return of the mission (Gallant et al. Reference Gallant, Ellery and Marshall2013). One or more scout rovers might be devoted to methane source localization while the main rover performs more generalized exploration. Any Mars aircraft will be challenging to implement due to the thin atmosphere, but there do exist marginal solutions to this approach, offering the prospect of vertical plume tracking (Fielding & Underwood Reference Fielding and Underwood2004). We do not consider this here. There is a diverse set of robotic plume-tracing algorithms that have been developed for terrestrial purposes, many of which assume laminar flow plumes. However, any local methane plumes on Mars are likely to be turbulent in nature.
The earliest robotic techniques for locating odour sources consisted of only instruments for measuring the concentration of the substance of interest to implement reactive responses. The Braitenberg architecture uses two chemical sensors directly connected to each wheel that proportionally drive their wheel motors, effectively gradient-tracing a plume (Braitenberg Reference Braitenberg1984). In turbulent plumes, the instantaneous gradient and the direction to the plume source are often not aligned, so gradient-descent algorithms will not per se yield the source. Further information is required to find the plume source. Robotic scent tracking may be aided by a wind sensor to measure the forcing function acting in the fluid. Knowledge of the forcing function allows information about the plume's structure to be inferred. Plume measurements can then be compared with models of the plume's structure and used to estimate the model parameters using a Kalman filter (Ishida et al. Reference Ishida, Nakamoto and Moriizumi1998). The goal is to predict the most likely location of the plume source given inferred wind measurements from a moving vehicle. Average plume concentration 
$\bar C$ flows over time in the x direction (centreline) and spreads laterally in the y direction given by the static Sutton Gaussian plume model (Sutton Reference Sutton1953):
$$\bar C(x,y) = C_{\max} {\rm e}^{ - y^2 /2{\rm \sigma} _y^2}, $$
where 
$C_{\max} = (Q/\sqrt {2\pi} \sigma _y u)$ = plume centreline concentration, 
$\sigma _y = \displaystyle{{1} \over {2}}Ax^{(2 - n)/2} $ = dispersion in the y direction (plume width) (m2 s−1), Q = source strength (kg s−1), u = u x = wind speed in the x (centreline) direction (m/s), A = a constant, n = unitless parameter that governs the shape of the plume, most commonly, n = 1 providing the best match to natural plumes. This model predicts time-average plume behaviour, but instantaneous odour distribution differs from the time-averaged plume. Hence, a chemical filament may be modelled in terms of fluctuating location:
$$\dot p = u(\,p) + n(t)$$where p = (x,y) = filament location, u(p) = (u x,u y) = average wind velocity, n(t) = (n x,n y) = time-variable Gaussian noise. Integration yields instantaneous filament location at time t f after emission from source location p 0 = (x 0,y 0) at time t i = t 0:
$$p_{t_{\rm f}} = \int_{t_i} ^{t_{\rm f}} {u(\,p(t)){\rm d}t} + \int_{t_i} ^{t_{\rm f}} {n(t){\rm d}t} + p_0. $$From f-1 discrete flow velocity measurements rather a continuous model, the source location at t 0 is given by:
$$p_0 = p_{t_{\rm f}} - \sum\limits_{i = 1}^{\,f - 1} {u(\,p(t_i ))\delta t} - \int_{t_i} ^{t_{\rm f}} {n(t){\rm d}t}. $$Unfortunately, the time of emission t i = t 0 is unknown but it is assumed that emission is continuous. A local ensemble transform Kalman filter has demonstrated that combining discrete temperature measurements from the Mars Global Surveyor (MGS) thermal emission spectrometer with a Mars GCM at different scales gives enhanced forecasts of both temperature and winds (even though the latter were not measured within the Kalman filter) (Hoffman et al. Reference Hoffman, Greybush, Wilson, Gyarmati, Hiffman, Kalnay, Ide, Kostelich, Miyoshi and Szunyogh2010). This demonstrates the potential for advanced fluid forecasting. However, the Kalman filter approach works most effectively with laminar plumes or over long-duration averaged measurements because of its reliance on linear models (though unscented Kalman filters relax this constraint (Wan & van der Merwe Reference Wan and van der Merwe2000). A reactive plume source localization system based on naive physics with a bounded map of air flow was demonstrated to control instantaneous heading and velocity of vehicles (Kowaldo & Russell Reference Kowaldo and Russell2006). This Kowaldo–Russell class of algorithms implemented an infotaxis strategy (Vergassola et al. Reference Vergassola, Villermaux and Shraiman2007).
The Farrell–Pang source likelihood map (SLIM) was based on hidden Markov models (HMM) and was subsequently improved using Bayesian methods. HMM predict the likelihood of an odour plume according to position and its source based on measurements of plume concentration history and flow velocity (Farrell et al. Reference Farrell, Pang and Li2003). The plume model predicts the probability of observed concentration measurements using Bayes rule. SLIM constructs a probabilistic estimate of source location. Plume detection and non-detection events are used to update the SLIM using Bayes theorem. SLIM updates the source likelihood vector based on the probability that odour migrates between specified cells at a given time (Pang & Farrell Reference Pang and Farrell2006). These methods were found to be reliant on a path planner to ensure that the vehicle continued to move in the direction of the source. With SLIM, the robot is implicitly constrained in its search for the source to a restricted area in its map (Ellery et al. Reference Ellery, Nicol and Cloutis2012). Boundary constraints are used for terrestrial applications when it is assumed that some previous information is available about the source being searched for. In planetary exploration, such constraints may not be available due to a lack of prior information. Any information about Martian methane plumes from remote sensing data will likely have resolutions of ~km (assuming that IR sensing will have inferior resolution to the current high visible image resolution of 20–100 m available from orbit) which is not well suited to local rover navigation scales. Probabilistic algorithms such as those involved in SLIM are much more computationally expensive than reactionary algorithms. Planetary rovers are heavily constrained by limited computational resources onboard. Hence, many current robotic state-of-the-art SLIM algorithms are not suitable for searching Martian methane plumes. Reactive algorithms on the other hand generate responses based on current measurements alone and have low computational overheads and so are more readily implemented on planetary rovers. Reactive plume tracing is common in nature. Lobsters forage by scent and their behaviours have been replicated robotically (Grasso et al. Reference Grasso, Consi, Mountain and Atema2000). Significantly, long range olfactory-based navigation has been observed in sea birds (Nevitt Reference Nevitt1999). Many insects forage for mates using olfactory methods in response to pheromones including moths (Mafra-Neto & Carde Reference Mafra-Neto and Carde1994), wasps (Kerguelen & Card Reference Kerguelen and Card1996) and mosquitoes (Zwiebel & Takken Reference Zwiebel and Takken2004). Moth tracking of odour plumes in search of potential mates is currently the best-documented biological scent-tracing behaviour such as (Farkas & Shorey Reference Farkas and Shorey1972) and Mafra-Neto & Carde (Reference Mafra-Neto and Carde1994).
Here, we describe an exploration of reactive control laws for a scout rover to attempt to autonomously locate a plume source of Martian methane. We had previously simulated the plume environment to capture the odour source localization problem and found a suitable algorithm for autonomous source localization via evolutionary robotics using genetic algorithms (de Croon et al. Reference de Croon, O'Connor, Nicol and Izzo2013). In this paper, we examine the rationale for the use of such genetic algorithms by exploring reactive control laws in a plume model. We begin with a detailed exposition on the methane plume model followed by a statistical analysis comparing our model with the ideal Sutton model. We then tested our reactive Braitenberg-type strategy for plume tracking to search for the plume source and found the technique lacking.
Plume model overview
Most models of long-term averages of atmospheric plumes are characterized as Gaussian, smooth and constant. The methane plumes that we expect to encounter on Mars are patchy, stochastic and time-variant. Within the plume, there are many empty spaces of clean air whose concentrations are that of the background. Unlike Gaussian plumes, the concentration gradient does not, in general, point in the direction of the source. Gradient-descent techniques are robust for static Gaussian plumes but fail to locate the source if plumes are turbulent. Jones (Reference Jones1983) adapted the statistical analysis of models of puffs of gas (Chatwin Reference Chatwin1982) and applied it to modelling turbulent air plumes validated by sensor measurements along the plume centreline. Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) used the same statistical scheme to validate their filament-based model of a moth's pheromone plume. Animals and insects must use anemotaxis rather than chemotaxis patterns to locate the source of a chemical under turbulent conditions. We must do similarly.
The concentration at any point in the plume over time will vary due to: (i) reactions that either release or destroy the substance, (ii) advection of the plume as a whole downwind and (iii) molecular and turbulent diffusion that is responsible for changing the shape of the individual puffs from which the plume is made. Reactions refer to the locations and rates where the plume chemical is created or destroyed. Advection is the process of the chemical being transported through a fluid flow (in this case, wind) and is dependent on the spatial chemical gradient. The wind will vary with both position and time generating turbulent diffusion. Molecular diffusion is the natural dissipation of concentrations due to molecular collisions. It is dependent on the second derivative of the concentration with respect to position. A concentration distribution over time with no reaction or advection would result in a uniform distribution as time tends to infinity. The classic diffusion-advection-reaction (DAR) equation (equation 5) models the temporal dynamics of the plume. It describes the change of methane concentration, C = C(p,t) as a function of position and time in terms of constant molecular diffusion constants in the x and y directions, D x and D y (it is assumed that D = D x = D y), instantaneous wind speeds in the x and y directions, u = (u x,u y), and the rate of methane production/destruction due to the reaction process of interest Q = Q(p,t):
$$\displaystyle{{\partial C} \over {\partial t}} = \displaystyle{{\partial} \over {\partial x}}\left({D\displaystyle{{\partial C} \over {\partial x}}}\right) + \displaystyle{{\partial} \over {\partial y}}\left({D\displaystyle{{\partial C} \over {\partial y}}}\right) - u_x \displaystyle{{\partial C} \over {\partial x}} - u_y \displaystyle{{\partial C} \over {\partial y}} + Q.$$We adopted a value of diffusion coefficient D of 0.01 m2 s−1 derived as a rough approximation of its computation from the Chapman–Enskog equation for a Martian near-equatorial summer:
$$\eqalign{D & =\displaystyle{{1.858x10^{ - 3} T^{3/2} \sqrt {1/M_{{\rm CH}4} + 1/M_{{\rm CO}2}}} \over {p\sigma _{12}^2 \Omega}} \cr & = 3.54 \times 10^{ - 3} {\rm m}^{2} / {\rm s}\sim 0.01\,{\rm m}^{2} /{\rm s}\;,}$$
where T = absolute temperature = 290 K, M = molar mass (g/mol) = 44.01 g mol−1 for CO2 and 16.04 g mol−1 for CH4, p = pressure = 0.0069 atm, 
$\sigma _{12} = \displaystyle{{1} \over {2}}\left( {\sigma _{{\rm CO}2} + \sigma _{{\rm CH}4}} \right)$ = average collision diameter = (3.3 + 3.8)/2 = 3.55 A, Ω = collision integral ≈ 1. Methane is released periodically at a frequency of 1/T where T = time period between pulse crests and d q = pulse duration such that d q < T. The amount q released at every time step in the solver is such that the entire impulse releases a total quantity Q of methane:
$$Q = \int_x {\int_y {\int_t {q\partial x\partial y\partial t}}}. $$We adopted a value of 1 kg s−1, consistent with the estimated value of 0.6 kg s−1 mentioned in the Introduction section. The DAR partial differential equation is numerically integrated over time and over a spatial mesh to determine the concentration at any given location within the mesh. We imposed a number of simplifications in solving the DAR equation that reflects certain characteristics of Martian methane plumes:
(i) the Martian rover is a surface traversing vehicle with no ability to alter its vertical position – we use a two-dimensional (2D) model only,
(ii) the simulation area is flat ground with no perturbations in the wind due to uneven ground or obstacles,
(iii) one of the anomalies in the Martian methane measurements is the rapidity at which methane is destroyed,
(iv) the wind model does not necessarily follow conservation laws because the model is designed to focus on replicating methane distributions.
The timescales for the simulation are much smaller than the methane lifetimes so no destruction processes are modelled. Rather, methane is conveyed out of the simulation environment by wind but no influx of methane was brought into the simulation environment from upwind sources. This imbalance acts as a sink and stabilizes the solver. Although recent work suggests that non-uniform vertical distribution of methane driven by large-scale Hadley circulation could provide data useful for methane location (though methane is uniformly mixed within a few weeks) (Viscardy et al. Reference Viscardy, Daerden and Neary2016), a surface rover will not have high resolution vertical measurements available for tracking.
The wind vector u is constructed from three components, u = a + v + τ where a = large-scale advection (base wind) responsible for moving the entire plume, v = medium-scale vortex function responsible for mixing individual puffs of methane and τ = small-scale turbulence responsible for individual motion of the methane puffs. Each of these represents a specific scale of transportation: advection affects the plume as a whole; vortex affects individual puffs and turbulence affects methane molecules. Advection varies in magnitude and direction over time causing the plume to form a sinuous shape. The advection function consists of a nominal (mean) wind direction and speed which varies ±30° over time. Vortices are local circular wind patterns that alter the shape of the puffs and move them cross-wind. Turbulence is a discrete, stochastic process added to the system. The wind is varied by superimposing effects from a continuously variable noisy stream function forming a 3D surface. Wind direction follows lines of constant altitude on the surface and the strength corresponds to the gradient. Equation (8) shows how the two components of the advection vector are determined from the stream function Z(x,y) (White Reference White1998):
$$a_x (x,y) = - \displaystyle{{\partial Z(x,y)} \over {\partial y}}\;\hbox{and}\;a_y (x,y) = \displaystyle{{\partial Z(x,y)} \over {\partial x}}.$$The advection stream function Z a(x,y) consists of the superposition of a linear plane representing the nominal wind and a noisy surface that adds texture. The noisy surface is computed using the Fourier series:
$$\eqalign{Z_a (x,y) & = \sum\limits_{n = 1}^{N} \left(A_n \sin \left(\displaystyle{{n\pi x} \over {L_x}}\right) + B_n \cos \left(\displaystyle{{n\pi x} \over {L_x}}\right)\right.\cr & \quad \left. + C_n \sin \left({\displaystyle{{n\pi y} \over {L_y}}}\right) + D_n \cos \left({\displaystyle{{n\pi y} \over {L_y}}}\right) \right),}$$where L x and L y are the x and y sides binding the simulation area, respectively. The Fourier coefficients (A n, B n, C n, D n) are varied by filtering white Gaussian noise through the closed loop transfer function H(s) = (G(s)/(s 2 + as + b)) where G(s) = open loop transfer function. Medium-scale mixing of methane puffs is modelled as vortices by adding a Gaussian pulse to the advection potential Z centred at the vortex location. The amplitude A, diameter d and position (x,y) of the vortices are specified by the vortex stream function Z v(x,y):
$$\eqalign{Z_v (x,y) & = A\left({\displaystyle{{\exp \left( { - 8(x - x_0 )/{\rm d}} \right)} \over {1 + \exp \left( { - 8(x - x_0 )/{\rm d}} \right)}}}\right) \cr & \times \left({\displaystyle{{\exp ( - 8(y - y_0 )/{\rm d})} \over {1 + \exp ( - 8(y - y_0 )/{\rm d})}}}\right).}$$Turbulence τ is added to the base wind as a stochastic change in x and y wind directions. To prevent the model from diverging, feedback is used to drive the wind model towards the base wind. The second derivative of the turbulence is calculated by multiplying a gain by a matrix of random numbers generated from a normal distribution R which is corrected by feedback from the error between the previous full wind model u = (u x,u y) and the current base wind model u m:
$$\displaystyle{{{\rm d}\tau} \over {{\rm d}t}} = K_{\rm 1} R_x + K_{\rm 2} (u - u_{\rm m} ).$$Wind speed on Mars varies with an average in the range 0–10 m s−1 – we used a low value of 1 m s−1 to ensure that the plume did not disperse too rapidly thereby overly taxing gradient-descent methods. The plume model parameters are summarized in Table 1:
Plume model parameters
Our plume model clearly indicates the wispy nature of the methane plume with significant gaps and puffs (Fig. 1) – simple gradient following would not be an effective search strategy for locating the source.
Snapshot concentration C (p) at arbitrary time t of our methane model in x and y directions (in metres).
Statistical methods of describing plumes
Most models of Martian methane plumes have been based on GCM with large plumes on planetary scales (Mischna et al. Reference Mischna, Allen, Richardson, Newman and Toigo2011). Rover missions do not cover such global scales but are localized to small regions. The Mars Exploration Rovers were initially designed for mission traverse distances of 600 m though they have both surpassed such distances to several kilometres (7.7 m for Spirit and to date 43.8 km for Opportunity over decadal timescales) – regional coverage over decades but still local over timescales of weeks. Models of methane plume distributions estimate that measurements which are not directly in the plume path would need to be made within a few metres of the plume axis to be discernible from background levels of methane (Olsen et al. Reference Olsen, Cloutis and Strong2012). We are, therefore, interested in local plume models rather than global models. In local methane plumes, the general planetary forcing function (advection) is no longer applicable. Instead, a turbulent environment is more representative of the situation in which a robotic rover would find itself while seeking a methane source. Many investigations on the structure of turbulent plumes exist. Well-known characteristics of chemical plumes include the tendency to meander with instances where local concentrations are much higher than the time mean (Murlis et al. Reference Murlis, Elkinton and Carde1992). Turbulence tends to create rotational, diffusive and stochastic effects on which a single particle would wander through an area of influence whose length is constantly increasing (Lumley & Panofsky Reference Lumley and Panofsky1964). This turbulence causes patches of clean air in the plume where concentrations of methane are indistinguishable from background levels. Intermittency is defined as the percentage of time where concentrations are below a pre-defined threshold – it can be caused by both meandering and patchiness of the plume, and in tests of several minutes, care must be taken to ensure that high values are descriptive of the plume and not the local wind pattern during that time (Murlis et al. Reference Murlis, Elkinton and Carde1992). We use similar statistical methods as Jones (Reference Jones1983) and Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) to compare our model to past experimental and simulated results.
Time-averaged plumes are generally modelled as Gaussian. Since both plume dynamics and Mars rover movements have short timescales, the Gaussian model does not capture enough detail for use as a suitable simulation environment. The Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) model represents the plume as a sequence of puffs. The realism of the presented model has a long-term average that closely replicates the Sutton model with its characteristic tear-shaped contour map (Sutton Reference Sutton1953). To correlate the average concentration 
$\bar C$(p) of the plume with the Sutton model, the mean concentration (once the methane is well mixed) is defined between time t i and the final time of the simulation t f as:
$$\bar C(\,p) = \displaystyle{{1} \over {t_{\rm f} - t_i}} \int_{t_i} ^{t_{\rm f}} {C(\,p,t){\rm d}t}.$$ The Jones Reference Jones1983 experiment emplaced four sensors at discrete locations inside the atmospheric plume located in the nominal wind direction: 2, 5, 10 and 15 m downwind from the source. The concentration measurements, Γr(p,t) were extracted over r runs to yield the mean measured concentration, 
$\bar \Gamma $(p) at a specific location:
$$\bar \Gamma (\,p) = \displaystyle{{1} \over {t_{\rm f} - t_i}} \int_{t_i} ^{t_{\rm f}} {\Gamma ^r (\,p,t){\rm d}t}. $$The mean concentration at point p = (x,y) in a Gaussian plume from the source at time t is given by:
$$C(\,p,t) = \mathop {\lim} \limits_{n \to \infty} \displaystyle{{1} \over {n}}\sum\limits_{r = 1}^n {\Gamma ^r (\,p,t)} = \bar \Gamma (\,p,t),$$
where 
$\bar \Gamma (p,t)$ = summation of measured intermittent fluctuations. Both C(p,t) and 
$\bar \Gamma (p,t)$ are, in fact, equal under certain conditions (Jones Reference Jones1983). We also define 
$\tilde \Gamma $(p) = the conditional average concentration above a sensor measurement threshold where Γ(p,t)≥γ where γ = pre-defined threshold (determined by the sensor instruments). The average concentration 
$\bar \Gamma (p)$ of the plume matched closely with the statistics when the integral of equation (14) was calculated only during the times when the sensor was in the plume (i.e. when the concentration 
$\tilde \Gamma (p)$ was above a defined threshold). Table 2 summarizes the mean data from both Jones’ (Reference Jones1983) and Farrell et al.’s (Reference Farrell, Murlis, Long, Li and Ring2002) results:
Jones (Reference Jones1983) and Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) plume data
Jones’ (Reference Jones1983) experiments involved measurement of ionized air to trace coronal discharges from thin wires as plumes. Coronal density was measured as charge concentrations (C/m3). For this reason, we are interested in the relative patterns of the experimental results rather than the absolute values. Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) modelled the emission of pheromones into the air from the female gypsy moth (in molecules cm−3). This was dependent on the puff size. Although this makes direct comparison difficult, general trends can be ascertained. In both cases, the conditional means were greater than the general means for each location due to the removal of zero-concentration events from the conditional mean. Furthermore, the mean decreases as the sensor position was moved further from the plume unsurprisingly. Table 3 presents the mean concentration and the mean conditional probability of being inside the plume, comparable to the same data given by Jones (Reference Jones1983) and Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) summarized in Table 2.
Model simulated general mean and the conditional in the plume mean at four positions
Five statistical parameters may be employed to describe the methane plume. The nth central moment M n(p) may be recast as relative intensity of fluctuations to allow them to be compared across studies (equation 15):
$$M_n(p) = \displaystyle{{\left({\displaystyle{{1} \over {t_{\rm f} - t_m}} \int_{t_m} ^{t_{\rm f}} \left({{\mathop{\Gamma}\limits^{\leftharpoonup}}(p) - \Gamma (\,p,t)}\right)^n \partial t}\right)^{1/n}} \over {{\mathop{\Gamma}\limits^{\leftharpoonup}}(p)}}.$$ Relative intensity, skewness, kurtosis, etc. used to quantify peaking behaviour are derived from these moments. For n = 2,3,4, M n(p) yields three parameters: standard deviation σΓ, skewness S Γ and kurtosis K Γ, respectively. The final two of the five parameters are: (Γˆ/ 
$\tilde \Gamma $) = peak-to-mean ratio where Γˆ(p) = max(Γ(p)), and I = intermittency defined as the percentage of time for which Γ(p,t)≥γ. Without loss of generality, we assume that γ = 0. The Jones (Reference Jones1983) and Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) experiments and models determined the following values for these five parameters (Table 4).
Plume amplitude data for five statistical parameters at three positions
The trends with distance for each parameter in Table 3 are similar between our model and the experiment. However, the Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) model does not match well the intermittency in the Jones (Reference Jones1983) experiment. Re-creating the 85% intermittency from the Jones experiment was, however, a challenge during the design of our model.
Model analysis
We performed multiple simulations to ensure that they performed similarly under different initial conditions and with different random number seeds. At each instant during the simulation, the plume is characterized as sinuous and patchy. The general mean concentration map of the entire simulation, however, closely resembles a Gaussian plume model. Figure 2 compares the concentration maps between the instantaneous and the time-averaged cases.
Long-term average results from the simulation with the plume source at (10,0) in 2D space (in metres). The solid lines represent isopleths for the average mass concentration of the dynamic model where the dash line represents a best-fit Sutton model. The isopleths have concentration values of 0.25, 0.5, 1, 1.5 and 2 kg m−2 progressing closer to the source.
Figure 2 shows isopleths from the time-averaged map and the Sutton plume as a Gaussian plume for comparison (Sutton Reference Sutton1953). The values used were Q = 2.11, u = 0.90 and C y = 0.54. The Sutton model approximates well the time-averaged plume in the simulation. Whereas Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) show only one isoline of concentration, Fig. 1 shows multiple isopleths which better illustrates the comparison between the averaged dynamic model and the Sutton model.
As would be expected, the conditional means are greater than the general mean concentrations at each location. The means from Jones (Reference Jones1983), Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) and our model are also plotted in Figs. 3 and 4.
Relative mean concentration over four positions for three data sets (y axis is unitless).
Conditional mean concentration over four positions for three data sets (y axis is unitless).
The data are normalized to allow comparison between the two models and the experiment based on the results at the 2 m position. The current model follows the downward trend of the other two plots, but does not match the experimental data as well as the model of Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002), especially for the conditional mean. This discrepancy is because in both Jones’ experiment and Farrell's model, the intermittency for all positions was similar, whereas in our model the intermittency at 2 m is significantly lower than the other positions. Since this trend greatly affects both means, it is the likely source of the discrepancy.
The amplitude statistics are presented in Table 5 for comparison to the unfiltered values summarized in Table 3 in addition to the filtered values published in Jones (Reference Jones1983) and Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002). For consistency, the results are presented for the four positions and six filtering conditions. The filtering conditions were implemented to eliminate high-frequency effects to different degrees and explore their effects on the results. The purpose was to emulate the effects of real sensors which will in general implement low-pass filtering to reduce the effects of high-frequency noise, which prevents useful integration of signals. There was an effect on reducing the values of most of the statistical parameters, but they were smooth and gradual. Furthermore, although most of the statistical trends were retained, the intermittency decreased rapidly as the high-frequency threshold was reduced in frequency – this is to be expected since the threshold is affected by the cut-off frequency of the filter. Trends for many of the parameters fit Jones’ experiment. These trends include a decrease of all parameters as the filter bandwidth decreases and the increase of the skewness, kurtosis and peak-to-mean ratio as the distance from the source increases.
Five statistical parameters (standard deviation, skewness, kurtosis, peak-to-mean ratio and intermittency) for six models at four locations along the nominal wind position for a 10 min long simulation
Our model did not replicate the consistent levels of intermittency throughout the experiment exhibited by Jones (Reference Jones1983) and Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002). The 2 m downwind location had a much lower percentage of intermittency than the other three positions. In Jones’ (Reference Jones1983) experiment, the levels of intermittency were fairly consistent. For Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002), intermittency increased 20% for the three data points, a range that was twice as wide compared with Jones’ data for the same distance covered. Once the plume reaches the 5 m position, intermittency stays very consistent (a range of 3% over 10 m). Even so, the intermittency was still much lower than that of both the experimental data of Jones (Reference Jones1983) and the model of Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002). This is not altogether surprising as we were using Martian parameters with lower atmospheric density, a high diffusion coefficient and a higher dispersion rate. Indeed, this may be considered as one of the useful results of these Martian simulations in how Martian plume behaviour will be different than that on Earth, especially in relation to methane emissions. Another difference in results is the odd results that Jones (Reference Jones1983) gets for the 15 m downwind position. In his experiment, the skewness, kurtosis and peak-to-mean ratio all tend to increase with distance downwind except for the final position at 15 m where it decreases. The standard deviation is the only parameter that does not exhibit this reversal. Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) have not analysed their results at the 15 m downwind position though their results match the trends for five data points until the 15 m position. Our model shows that they do not reverse trend at the 15 m point for any of the parameters.
Reactive control behaviours
We consider two reactive control methods – a chemotaxis gradient-based search method and an anemotaxis method which relies on knowing the local wind vector – to find the plume source. Gradient methods become less effective at locating chemical sources as the Reynolds number of the transporting fluid increases. Reynolds number is defined as 
${\mathop{\rm Re}\nolimits} = (\rho vL/\mu )$ where ρ = fluid density = 0.01308 kg m−3, v = flow velocity = 1 m s−1, L = characteristic length = 1 m, μ = fluid dynamic viscosity = 1.422 × 10−5 kg ms−1 assuming Mars surface conditions, relative flow velocity of 1 m s−1 and characteristic rover dimension of 1 m, giving a low Reynolds number of 9200. Gradient-based plume tracing is effective only in low Re fluids where molecular diffusion is the dominant source of dispersal. High Re number implies that plume dispersion is dominated by turbulence which stretches and twists plume filaments causing increased intermittency in the plume. There is no single Reynolds number for Earth or Mars atmospheres as it depends on fluid velocity, atmospheric density, characteristic length and atmospheric viscosity – indeed, the Mars surface atmosphere would have a much lower Reynolds number than Earth sea level due to its much lower density making gradient methods more tolerant on Mars. The gradient method used in this work is implemented through the classical Braitenberg vehicle architecture. The two sets of wheels on each side of the mobile robot are controlled via a chemical sensor on the opposing side of the vehicle. The speed of each side, v l and v r, is determined by a proportional gain, K, on the chemical concentration measurements on the left side, C l, and the right side, C r, of the vehicle. The forward and angular velocities (v,
$\dot \theta $) are calculated from the wheel speeds and the rover width b:
$$v_{\rm r} = KC_{\rm l}, $$
$$v_{\rm l} = KC_{\rm r}, $$
$$v = \displaystyle{{1} \over {2}}(v_{\rm r} + v_{\rm l} ),$$
$$\dot \theta = \displaystyle{{1} \over {b}}(v_{\rm r} - v_{\rm l} ).$$We initially tested the Braitenburg vehicle on a static Gaussian plume of the Sutton (Reference Sutton1953) model using the parameters in Table 6.
Braitenberg parameters for testing on the Sutton plume model
The vehicle navigated the Sutton plume for various initial conditions: the initial position was varied between 40 and 60 m downwind and up to 15 m lateral to the centreline in either direction. The initial heading of the rover was always in the upwind direction at an angle in the range ± 2π/5(72o). The vehicle was able to navigate to the source in all 1000 trials as expected since the plume was static and did not change between the trials. A representative path showing the navigation of the Sutton plume is shown in Fig. 5.
Results of Braitenberg vehicle: representative path of the Braitenberg vehicle navigating the static Sutton plume in x and y coordinates (in metres). Isolines of plume concentration are shown in solid lines with the concentration value showed. The dash line is the path of the vehicle with the starting point marked as a square.
The Braitenberg vehicle was also tested on a dynamic plume model with representative paths through the plume simulation shown in Fig. 6 with the time-averaged plume in the background. The simplified gradient method used in the Braitenberg vehicle cannot find the methane source due to the added complexity of the turbulent plume. Rather, it tends to follow the puffs of methane as they are blown downwind causing the vehicle to diverge from the source.
Three representative paths taken by the Braitenberg vehicle through the dynamic plume model in x and y coordinates (in metres). The contour of the time-averaged plume is in the background. The rover paths are the solid lines. The termination conditions shown were either plume exit and/or migration away from the plume source for clarity – however, simulations were extended beyond these conditions indicating that these trajectories were pathological and never recovered ‘the scent’.
Wind speed was measured relative to the rover rather than wind ground speed. Mars rovers drive very slowly – 0.03 m s−1 top speed nominally but they usually drive slower depending on rock distribution, slope and regolith properties. This is two orders of magnitude less than ground wind speed, so it can be ignored relative to wind speed. We adopted a 1 m s−1 wind speed but 10 m s−1 is common on Mars even not accounting for the ubiquitous dust devils, which are likely to be highly dispersive. In running the rover trajectories, they were in the plume for extended periods of time – the reactive gradient-based approach was unable to take advantage of the accumulated data. This graphically illustrates that any future Mars rover mission that might seek to find the source locations of methane emissions must employ sophisticated plume tracking, mapping, modelling and predictive capabilities. Candidate approaches include SLIM and variations thereof, but these techniques are computationally intensive and are unlikely to be employable with near term rover-based computational capabilities.
Conclusions
This paper presents a model of chemical dispersion which captures the difficulty in Mars methane source localization. The model is consistent with Jones (Reference Jones1983) data unlike the Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) model. Jones (Reference Jones1983) and Farrell et al. (Reference Farrell, Murlis, Long, Li and Ring2002) dictate that Mars methane plume dynamics are likely to be complex. For this reason, traditional gradient-search methods valid for time-averaged plumes fail to find the source of such dynamic wispy plumes. More sophisticated robotic strategies are required, such as SLIM, but severely restricted computational resources onboard current rovers make this unlikely. However, one promising approach is that high computational capacity field programmable gate arrays are becoming space qualified. A variation on SLIM would be through the incorporation of large-scale Hadley-driven movement from orbit (Viscardy et al. Reference Viscardy, Daerden and Neary2016) with surface rover local measurements at the surface using traditional Kalman filter-based SLAM (self-localization and mapping). Of course, this assumes that methane detection sensors and rover autonomous navigation are sufficiently sensitive and responsive, respectively, to support such a methane tracking capability. Multiple rovers might be deployed to provide highly parallel search and tracking facilities – a sensor network effectively to implement network science. It may be arguable that every Mars mission should include methane sensors to complement temperature, pressure and wind vector sensors that are currently on almost every Mars lander/rover mission payload manifest. The prospect for network science with rover fleets on Mars is expensive and therefore bleak. Furthermore, this does not solve the basic problem that computational resources will always be in demand by all services onboard a rover and plume tracking will require SLIM-type capabilities.
In future work, we would recommend the adoption of deep learning techniques which require large data sets from which to train (it is unclear if sufficient data exists and the construction of such data sets would be challenging in terms of effort if they must be generated through experimental means) (Schmidhuber Reference Schmidhuber2015). However, we have generated simulated training data quite readily and successfully through models to compensate for sparse data sets in a different context (Cross et al. Reference Cross, Ellery and Qadi2013). The construction of plume models to output position-concentration data maps should be feasible and ideally suited to convolutional neural networks. The advantage of neural networks is that once trained, they consume minimal computational resources, processing and memory – it is the training process that would be implemented on the ground that is cumbersome. The use of a priori trained neural nets offers the possibility of implementing highly sophisticated algorithms without the traditional computational penalties. We have used Kalman filter-trained neural networks for SLAM-type tasks in rover navigation (Hewitt et al. Reference Hewitt, Ellery and de Ruiter2017) which could be applied in unscented form to plume tracking. Neural fields are spatial neural networks that can implement complex spatiotemporal patterns (Coombes Reference Coombes2005) that may be applied to plume tracking. We would recommend using more sophisticated leaky integrate-and-fire neuron models than the traditional simple switching neuron in order to explore whether stochastic resonance could be exploited in ‘noisy’ plumes (Castro & Sauer Reference Castro and Sauer1997; Mitaim & Kosko Reference Mitaim and Kosko1998). Scent-based foraging behaviour in higher animals beyond insects may offer insights into such neural learning of maps in the mammalian hippocampus. Furthermore, there may be similarities with mammalian eye movements which are characterized by the search for high information parts of a scene (saliency) and a random component (saccades) (Itti et al. Reference Itti, Koch and Niebur1996). There are a host of possibilities for solving the plume source tracking problem particularly bio-inspired approaches.
Acknowledgements
This work was supported through the Scent of Science project for the European Space Agency's Advanced Concepts Team and the National Science & Engineering Research Council of Canada Canadian Astrobiology Training Programme.
