MarsWRF Convective Vortex and Dust Devil Predictions for Gale Crater Over 3 Mars Years and Comparison With MSL‐REMS Observations

Convective vortices and dust devils have been inferred and observed in Gale Crater, Mars, using Mars Science Laboratory (MSL) meteorological data and camera images. Rennó et al. (1998) modeled convective vortices as convective heat engines and predicted a “dust devil activity” (DDA) that depends only on local meteorological variables, specifically the sensible heat flux and the vertical thermodynamic efficiency which increases with the pressure thickness of the planetary boundary layer. This work uses output from the MarsWRF General Circulation Model, run with high‐resolution nests over Gale Crater, to predict DDA as a function of location, time of day, and season, and compares these predictions to the record of vortices found in MSL's Rover Environmental Monitoring Station pressure data set. Much of the observed time‐of‐day and seasonal variation of vortex activity is captured, such as maximum (minimum) activity in southern summer (winter), peaking between 11:00 and 14:00. However, while two daily peaks are predicted around both equinoxes, only a late morning peak is observed. An increase in vortex activity is predicted as MSL climbs the northwest slopes of Aeolis Mons, as observed. This is attributed largely to increased sensible heat flux, due to (i) larger daytime surface‐to‐air temperature differences over higher terrain, enhanced by reduced thermal inertia, and (ii) the increase in drag velocity associated with faster daytime upslope winds. However, the observed increase in number of vortex pressure drops is much stronger than the predicted DDA increase, although a better match exists when a threshold DDA is used.


Introduction
Convective vortices occur during periods of strong convective heating of the surface, when the ground temperature exceeds the air temperature, warming the air above the surface. As this air rises, existing vorticity becomes more vertical and intensifies, with the air spiraling around a low-pressure region that develops in the vortex core (Balme & Greeley, 2006;Kanak, 2005;Spiga et al., 2016;Toigo et al., 2003). Rotation is randomly clockwise or counterclockwise, and vortices often occur in pairs or clusters (Balme & Greeley, 2006). They may extend to at least the height of the planetary boundary layer (PBL), which commonly extends to ~2-3 km in Earth's deserts and to~10-12 km over many parts of Mars and significantly higher in some regions; dust devils with heights up to 16.5 km have been identified in some regions (Fenton et al., 2016;Fenton & Lorenz, 2015). Dust devils are vortices that contain dust, making them visible; this dust is raised via strong tangential winds around the vortex core, assisted by the "suction effect" of the pressure drop and other thermophysical factors (see, e.g., Neakrase et al., 2016 for a review of all processes that have been proposed). Mars vortices, hence dust devils, can grow much larger than on Earth. Hence while the largest on Earth may reach a few tens of m in diameter, on Mars they may reach or order a km in diameter (Fenton et al., 2016). Detection of vortices is typically accomplished by detecting their pressure drop signature, with wind (speed and direction) and temperature also measurably affected when a vortex passes close enough to the sensors (e.g., Kahanpää et al., 2016;Murphy et al., 2016; see also section 2.1). In the case of dust devils, detection is also possible via imaging or by measuring changes in incoming solar radiation if the vortex is dusty enough (e.g., Lemmon et al., 2004; see also section 2.2).
Convective vortices are of particular interest on Mars because dust devils are thought to be one of the primary means of injecting dust into the thin Martian atmosphere, in which the thermal structure and hence circulation are strongly impacted by its presence (Gierasch & Goody, 1972;Smith, 2004). Note that Mars non-convective vortices have also been inferred from nighttime pressure data (e.g., Ordonez-Exteberria et al., 2018) and are likely a result of mechanically forced turbulence. However, they are less likely to raise dust high into the atmosphere or noticeably increase the atmospheric dust loading; the PBL depth is very small (hence vertical mixing is strongly inhibited) between approximately an hour after sunrise and an hour before sunset (e.g., Fonseca et al., 2017), hence any dust lifted is likely returned rapidly to the surface. Observations of active dust lifting centers during dust storm onset show little correlation with observations of where and when dust devils and their tracks occur (Cantor et al., 2006). Similarly, theory and modeling studies all suggest that fewer dust devils should occur as the atmosphere becomes dustier, which is in contrast to the rapid ramping up of dust injection as dust loading increases during storm onset (e.g., Kahre et al., 2006;Newman et al., 2002aNewman et al., , 2002b. However, it is estimated that dust raised by dust devils or other forms of small-scale convective dust lifting accounts for most of the dust raised outside of the "storm season" (Ls~0-180°) and hence that these processes are responsible for the observed temperatures at this time of year, which are tens of K warmer than a "clear" atmosphere would be (Basu et al., 2004;Kahre et al., 2006;Newman et al., 2002aNewman et al., , 2002b. Rennó et al. (1998) (henceforth R98) provided a theory that models convective vortices as convective heat engines and predict a "dust devil activity" (DDA) that depends only on local meteorological variables. Note that despite the name, this theory applies equally to vortices with or without dust. The theory of R98, described in section 3.2, has been used in several Mars General Circulation Models (GCMs) to predict the amount of dust lifting by dust devils; lifting set proportional to the DDA appears to successfully capture the global "background" dust amount in the atmosphere as a function of season when no major dust storms are present (Basu et al., 2004;Kahre et al., 2006). The theory has also previously been tested by comparing predictions of global models both with the spatiotemporal variation of dust devils and their tracks observed from orbit (see section 1.1.1) and with the variation of imaged dust devils and convective vortices (inferred from pressure data) observed from the surface (see section 1.1.2).

Theory of Convective Vortices and Dust Devils
The DDA is a measure of the energy that may be harnessed by vortices, but how this should be distributed into a vortex distribution is not clear at present. For example, should an increase in DDA result in more vortices, or an increase in the maximum vortex strength, or both? To date, comparisons with observations have largely focused on the number of vortices (or tracks) detected, but section 5.3 also investigates the relationship to pressure drop magnitude. Kahre et al. (2006) found that the spatial and seasonal patterns of DDA predicted using a GCM were broadly consistent with observations of dust devils and their tracks from orbit (Fisher et al., 2005): both show peak activity during local summer in both hemispheres, with Amazonis marked as a particularly active region. However, orbital data sets at the resolution required have limited spatial and seasonal coverage and do not typically provide information on time-of-day variation, as observations are often made at the same local time each orbit. In particular, afternoon observations are far more common than morning ones, although, e.g., the Colour and Stereo Surface Imaging System on the Trace Gas Orbiter is able to cover the full daytime portion of the diurnal cycle due to its non-Sun-synchronous orbit (Thomas et al., 2017).

Comparisons of Theory With Orbital Observations
A further issue is that smaller dust devils cannot be seen from orbit and dust devils with relatively low dust content may be missed. Also, dust devil tracks or even dust devils themselves may not be visible over certain surfaces, e.g., if the surface revealed when a layer of dust is removed has the same color and albedo, as might be true of a thick dust deposit. These arguments may explain why Kahre et al. (2006) found that small but nonzero DDA values are predicted for some seasons and locations at which no dust devils or tracks are observed , although this may also indicate a threshold DDA is needed for vortices to form or dust lifting to occur; see also section 5.3. In addition, if there is no loftable dust in a region then no dust devils or tracks will be seen, even if numerous strong vortices occur. There is no way to infer the presence of dust-free vortices from orbit, hence orbital observations are biased against regions with low dust availability. Orbital data sets of dust devils and their tracks thus contain numerous biases and are not ideally suited to testing theories of vortex or dust devil number, strength, timing of occurrence, etc.

Prior Comparisons of Theory With Landed (In Situ) Observations
Unlike orbital data sets, in situ measurements do not suffer from most of the above issues. On the surface, both dust-filled and empty vortices are most straightforwardly detected by measuring the variation of surface pressure, with rapid but short-lived pressure drops indicating the close passage of a vortex's low-pressure core (e.g., Murphy & Nelli, 2002;Schofield et al., 1997). Surface imaging of dust devils can then be used to cover a wider area; to provide estimates of the speed, diameter, and height of vortices; and to assess their dust-raising ability. The main disadvantages of in situ measurements are that spatial coverage is far more limited than from orbit: dust devils must be visible from the rover or lander, with increasing bias toward detecting larger and/or dustier dust devils as the distance increases; similarly, vortices can only be detected in pressure data if they affect pressure at the landed location. For pressure detections, it can be difficult to deconvolve the effects of strength, diameter, speed of passage, and distance, especially without reliable wind data or imaging, although several efforts have been made to do so (e.g., Jackson et al., 2018;Kurgansky, 2019;Lorenz, 2016). Yet overall, surface missions have the potential to capture the seasonal and diurnal variation of nearly the full strength spectrum of dust devils and vortices in a single location or region. For landers, imaging can also provide some limited information on any spatial variations in dust devils across the scene, while for rovers, spatial variations in both dust devils and dust-free vortices may be measured as the location changes over the mission. Kahre et al. (2006) found that the predicted diurnal variation of DDA at the Mars Pathfinder (MPF) landing site using their GCM's output was very consistent with the diurnal variation of daytime convective vortices inferred from MPF pressure data (Murphy & Nelli, 2002), with the peak occurring between noon and 1 pm LTST. Using a different GCM, Chapman et al. (2017) also noted a general agreement between predicted DDA and the observed diurnal cycle of number of convective vortices at the MPF and also at the Phoenix Lander sites (Ellehoj et al., 2010), although in both cases the model appeared to predict a slightly later peak time of day. However, Chapman et al. (2017) found a mismatch at the Viking Lander 2 site, for which their GCM predicted peak DDA in the late afternoon (around 17:00), whereas observations of wind fluctuations interpreted as being caused by vortices peaked between 10 am and noon (Ringrose et al., 2003). They also found a mismatch for Gale Crater, with their GCM predicting both a morning and afternoon peak in DDA over the first Mars year, whereas the average diurnal cycle found from pressure drop observations over the same period shows a single peak around noon . In fact, there appears to be a strong observed seasonal variation in the diurnal cycle in Gale Crater, as discussed in sections 5.1 and 5.2 below.
Prior to the current work, Kahanpää et al. (2016) was the only study to have compared the seasonal variation of predicted DDA with surface observations of seasonal vortex activity, as inferred from the first year of surface pressure data taken from the Mars Science Laboratory (MSL) Curiosity rover. Using output from a highresolution, nested mesoscale model (MarsWRF, which is also used in the current work, see section 3), they found an overall good match between the predicted DDA and observed seasonal cycle. The main exception was around southern summer solstice, during which time far fewer vortex pressure drops were inferred than predicted. The possible cause of this mismatch is discussed further in section 5.1.

Prior Modeling of the Circulation in the Gale Crater Region
Unlike prior landing sites where dust devils have been studied, the Gale Crater site of the MSL mission has significant mesoscale topographic relief. While direct observations of circulation are limited to the MSL instruments, the nature of the circulation has been modeled by several different groups (Tyler & Barnes, 2013Rafkin et al., 2016;Pla-Garcia et al., 2016;Newman et al., 2017;Fonseca et al., 2017;Richardson & Newman, 2018). The strong relief generates complex mesoscale circulations that interact both with the large-scale flows commonly observed at other landing sites and with the microscale (PBL) circulation and structure. Depending on the time of year and the strength of the large-scale flows, the circulation that the MSL rover observes can be dominated by flows induced within the crater itself. A significant true buoyancy slope flow is associated with near-surface upslope (anabatic) daytime and nighttime downslope (katabatic) winds (Tyler & Barnes, 2013), though this flow is significantly modified in practice by net hydrostatic rebalancing of the atmosphere between low and high topography as the near-surface scale height dramatically changes during the daily cycle Richardson & Newman, 2018).
The lateral advection of warm air "up slope" along with adiabatic cooling of this air will tend to reduce the thermal drive for convection. In addition, increased winds tend to generate additional near-surface turbulence through shear, which has the effect of cooling the surface through increased surface layer sensible heat fluxes; in more extreme form, this can be seen in the flows over the walls of Valles Marineris in numerical simulations, where the diurnal cycle of ground and near-surface air temperatures can be greatly reduced (Spiga et al., 2011). Both of these mechanisms will tend to decrease convection over the slopes of Aeolis Mons. By contrast, the mesoscale circulation set up within the crater during the daytime tends to converge over Aeolis Mons and over the crater rim, with downwelling over the crater trench. When it is dominant, this has the effect of suppressing the PBL depth over the trench (Tyler & Barnes, 2015) and hence damping PBL motions such as convective vortices there.
However, as described in some detail in section 4.3, the mesoscale flow is not always dominant, and in some cases the influence of large-scale flows on Gale Crater becomes significant (Newman et al., 2017;Pla-Garcia et al., 2016;Rafkin et al., 2016). It is thus important not to get too attached to a simplistic, idealized conception of the Gale Crater flow system and specifically the idea that the daytime boundary layer is always suppressed in the Gale Crater trench during the daytime.

Overview of the Work Presented Here
Dust devils and convective vortices have been observed in great numbers in Gale Crater by instruments on MSL, with their frequency varying strongly with location, time of day, and season. This data set thus provides an excellent test of theories and models that predict the intensity of convective vortex activity and the amount of dust raised by them. This paper focuses on comparing observations with thermodynamically derived predictions of DDA, obtained by combining the theory of R98 with meteorological variables output from the MarsWRF nested mesoscale model, run at a grid spacing of~1.4 km over Gale Crater.
Section 2 describes how MSL instruments are used to detect vortices and dust devils, discusses previously published results, and finally describes the methods used here to infer vortices from pressure data. Section 3 describes the MarsWRF model and how its output is used to predict the DDA, as defined in R98. Section 4 presents MarsWRF predictions across Gale Crater as a function of local time at local summer solstice, as well as seasonal variations at selected times of sol for the region in which MSL operates. Section 5 compares vortices measured by MSL's Rover Environmental Monitoring Station (REMS) with the DDA predicted at the time and location of the observations and examines what determines the DDA variation in the model. Section 6 concludes.

Vortex and Dust Devil Observations Made by MSL
MSL carries several instruments that have been used to detect convective vortices and/or dust devils: a meteorological sensor suite in the form of REMS and cameras in the form of the rover's navigation cameras (Navcams) and mast-mounted science camera (Mastcam).

Vortex Observations with REMS
The REMS suite includes sensors for measuring atmospheric pressure, temperature, relative humidity, UV radiation flux, and wind speed and direction. Measurements are performed at a rate of 1 Hz for typically up to ten hours per sol, where a sol is a Mars solar day. This includes 5 minutes of observations at the start of every hour in Local Mean Solar Time (LMST) plus typically eight 1-hour "extended blocks" per sol, where an hour is defined as 1/24 of a sol and a minute as 1/60 of an hour. The cadence of REMS extended blocks is described in detail in, e.g., Newman et al. (2017), but basically results in all hours of the sol being observed at least once at 1 Hz in every 6-sol period and almost daily coverage of the hour containing noon Local True Solar Time (LTST), i.e., when solar insolation peaks.
When some part of a convective vortex passes over the rover, its signature may be seen in (i) a rapid, shortlived pressure drop due to the vortex's low-pressure core. Its signature may possibly also be seen in (ii) a rapid, simultaneous change in wind direction as the rotating vortex moves across the sensor position, (iii) rapid, simultaneous air temperature fluctuations associated with the warm core of the vortex and the mixing of air around it, and/or (iv) an associated, relatively short-lived decrease in UV flux if the vortex contains sufficient dust that it blocks some solar radiation from reaching the surface. The latter can be used to estimate the amount of dust present in the vortex, although a significant amount of dust must sit in the atmospheric column above the rover to cause such effects, which is likeliest for larger dust devils that pass between the rover and the Sun direction. If a convective vortex does not pass directly over the rover, many of these observations are more muted and can be impossible to differentiate from general turbulence in the area, although they may be used to support the argument for vortex passage when a pressure signal is detected.
While the temporal (diurnal and seasonal) coverage of REMS measurements is very good, their spatial coverage is limited to the location of MSL. Fortunately, MSL is a rover. Thus by April 2018, it had not only monitored three complete Mars years but had also driven more than 18.5 km (approximately 9 km as the crow flies) from its original location in the trench of the crater. As shown in Figure 1, over this period, MSL crossed through the Bagnold Dunes and began to climb the slopes of Aeolis Mons. This means that the REMS data set also contains information on the spatial variation of vortex activity, albeit only along the rover's traverse.

Dust Devil Observations With Navcam/Mastcam
REMS measurements provide information about vortex activity at the rover location, but it is difficult to extract information about their size or motion (especially when simultaneous wind data are absent) or how much dust they contain. By contrast, MSL's cameras can be used to image dust devils directly, yet no clear sightings of dust devils were made in the first 2 years of the MSL mission. Until shortly after sol 1520, Navcam dust devil searches looked north across the crater trench. In retrospect, this appears to have been a poor choice, as the DDA is predicted to be generally low there compared to other directions visible from that location (see section 4.1), although dust availability may also have been an issue. Subsequent searches looking toward Aeolis Mons proved far more successful, although Navcam 360°surveys are performed regularly in addition to longer movies to avoid directional bias creeping in again. In fact, as the rover has climbed higher, it has become common to see dust devils when imaging to the north, as this now looks back down the slopes. However, the result is that the list of dust devils found in surveys only extends back as far as sol 1607,~2.4 years into the mission. The period since that time now exceeds one Mars year but includes a global and regional dust storm in MY34, both of which strongly affected atmospheric dust opacities but were not included in the present MarsWRF simulations. Given the strong impact of atmospheric dust opacity on both the predicted DDA (which decreases as opacity increases; Newman et al., 2002a) and observed vortex and dust devil occurrence (which strongly decreased during the MY34 global storm; Guzewich et al., 2018), a comparison with imaging is therefore deferred to a future paper in which MarsWRF is forced with realistic observed dust opacity maps for a given year, rather than using a dust scenario that represents a nonspecific year with no major storms. (2016), Kahanpää et al. (2016), and Ordonez-Exteberria et al. (2018) have all previously examined vortex activity in Gale Crater using REMS data. In Steakley and Murphy (2016), the signatures of 245 convective vortices with pressure drops exceeding 0.3 Pa were found in the first 707 sols of the mission, equivalent to an average of~1 vortex per sol after accounting for the incomplete temporal sampling of REMS measurements. They also found that pressure drops peaked between 11:00 and 13:00 LTST, with vortices more frequent during local spring and summer than in fall and winter.

Steakley and Murphy
Approximately the same period was investigated by Kahanpää et al. (2016), who found similar results in terms of the timing and seasonal behavior of pressure drops exceeding 0.5 Pa. In addition, they compared these results to DDA predictions made using output from the MarsWRF model (see section 3 for more details). They found that the DDA calculated from MarsWRF output predicted the seasonal variation of observed pressure drops (> 0.5 Pa) relatively well, except for an overprediction of the peak around local summer solstice (Ls~270°), which they attributed to the locally high thermal inertia of the specific area where the rover was located over much of this period, as discussed in section 5.1. Finally, they also interpreted a vortex burst on sol 664 as possibly being due to a front passing over Gale Crater.
Two full years of REMS data were studied by Ordonez-Exteberria et al. (2018), who found 635 daytime pressure drops exceeding 0.5 Pa over 1417 sols, corresponding to an average rate of 1.2 events per sol when corrected for sampling. These drops were most prevalent near local noon, but they also noted several weak events not correlated with strong surface-to-air temperature differences, as well as noting a large number of nighttime pressure drops, which were interpreted as being due to mechanically forced turbulence. Nighttime pressure drops were especially strong in regions of rough or steep topography and showed seasonal dependences which may link them to variations in local-and regional-scale winds. Finally, they noted roughly a doubling in the number of day and night pressure drops from the first to second Mars year of the mission.

Methodology Used in This Work
A method similar to that described by Kahanpää et al. (2016) is used to identify passing vortices. Moving through the REMS pressure data set in 20-s sliding windows, a vortex candidate is identified if either of the following criteria are met: • Minimum pressure > 0.5 Pa lower than the mean of the previous and next 20-s intervals • Minimum pressure > 0.3 Pa lower and mean pressure > 0.1 Pa lower than the mean of the previous and next 20-s intervals Note that the first of these criteria was not applied by Kahanpää et al. (2016). The criterion was added to identify relatively strong pressure drops with very short duration. Such events appeared to be common during the second and third Martian year of the MSL mission, unlike during the first year. The second criterion enables the detection of pressure drops with durations longer than 60s, which is an advantage compared to the methods applied by Steakley and Murphy (2016) and Ordonez-Exteberria et al. (2018). A visual inspection of every pressure drop is then performed to identify and remove any that do not look vortex-like or which are "double counting" of the same vortex by the algorithm. Figure S1 shows some examples of pressure drops retained and removed during this step. In common with Kahanpää et al. (2016), weak pressure drops appearing simultaneously with the REMS UV sensor passing into or coming out of the shadow of rover structures are removed as these events are apparently caused by the "shadow effect" of the REMS pressure sensor, first reported by Harri et al. (2014) and noted also by Ordonez-Exteberria et al. (2018). Nevertheless, by contrast with Kahanpää et al. (2016), pressure drops with magnitude exceeding 0.8 Pa are retained.
Magnitudes and half maximum durations of the vortex candidates are determined by fitting a linear combination of a Lorentzian function and a line to the REMS pressure data as described by Kahanpää et al. (2016). Vortices that arrive in pairs or larger "bursts" are for the most part handled individually; however a few events with multiple pressure minima are fitted with a linear combination of two Lorentzian functions. After this, events with a magnitude smaller than a threshold value are discarded. In contrast to Kahanpää et al. (2016) and Ordonez-Exteberria et al. (2018), this study uses a threshold of 0.6 Pa instead of 0.5 Pa, the rationale being that the pressure data of MSL years 2 and 3 are sometimes noisier than those of year 1, complicating the identification of pressure drops smaller than 0.6 Pa. For completeness, the distribution of pressure drops is also found using a method very similar to that of Ordonez-Exteberria et al. (2018) and using a threshold of 0.5 Pa. While this second method results infewer long duration, large pressure drops being discovered, the overall pattern of diurnal and seasonal variability is very similar using either method (see section 5.1).

The MarsWRF Model and Its Configuration for This Work
MarsWRF is a multiscale model of the Martian atmosphere, described in Richardson et al. (2007), Toigo et al. (2012), Newman and Richardson (2015), . This work uses a very similar model configuration to that described in Newman et al. (2017), in which MarsWRF was run as a global 2°model with five "nested" higher resolution regions roughly centered on MSL's landing site, with each nest smaller than its parent and with three times the horizontal resolution; see Figure 14 of that paper. In this work, however, only four nests are used, giving five domains in total and a resolution of~1.4 km in the innermost nest (domain 5), which covers all of Gale Crater. In addition, the present work primarily uses results from vertical grid B described in Newman et al. (2017), as it was found to produce the best match to winds and aeolian features within Gale Crater. However, some results are also shown that use vertical grid A, which more finely resolves the lowest portion of the boundary layer, to demonstrate the impact of model setup (see section 5.4).
The time-varying, three-dimensional atmospheric dust distribution seen by MarsWRF's radiative transfer scheme is prescribed using the Mars Climate Database Mars Global Surveyor dust scenario, which was developed to be representative of a year without any major dust storms Lewis et al., 2001). Twelve simulations, each lasting eight sols with the first then discarded as a "spinup" sol, were run at areocentric solar longitude, Ls = 0,30,60,90,120,150,180,210,240,270,300,and 330°, to fully sample the annual cycle of solar forcing.

Most of the spatially variable surface properties come from Mars Global Surveyor Mars Orbiter Laser
Altimeter and Thermal Emission Spectrometer observations (Christensen et al., 2001). These include MOLA topographic height (Smith & Zuber, 1996), the MOLA intrashot surface roughness map (Garvin et al., 1999), and maps of albedo, thermal inertia, and emissivity (Putzig & Mellon, 2007). Additional high-resolution topography data, from a Mars Express, High Resolution Stereo Camera instrument Digital Terrain Map, were used over Gale Crater for domain 5 (Neukum et al., 2004;Jaumann et al., 2007). Although some local surface properties have been measured by MSL (e.g., thermal inertia inferred from the diurnal cycle of surface temperature), these are not necessarily representative of the surface properties even a few meters away, whereas vortex pressure drops detected by REMS may be caused by vortices that passed several hundred meters from the rover and originally formed at an even greater distance. Due to drivability requirements and science interest, the surface properties measured along the rover traverse are sometimes not even representative of the typical local region around it, although a lack of vortices over a long period when MSL explored "Yellowknife Bay" is indeed most likely attributable to much higher thermal inertias over that entire region (see section 5.1). Overall, however, combining locally and remotely measured surface properties appropriately into an improved map of the region around MSL would be very challenging. In addition, the overall pattern of DDA predicted in Gale Crater is found to be more controlled by topography than by variations in albedo or thermal inertia (section 4.2.2). Hence, despite the somewhat low resolution, only maps of surface properties derived from orbital data are used in this work.

Dust Devil and Vortex Theory and Combination with MarsWRF Output
In R98, convective vortices (and hence dust devils) are modeled as convective heat engines. The DDA is defined as: where F s is the surface sensible heat flux (the heat input to the base of the vortex) and η is the vertical thermodynamic efficiency of the dust devil convective heat engine (the fraction of the input heat turned into mechanical work). After some manipulation, η is given as: and where p s is the ambient surface pressure, p top is the ambient pressure at the top of the convective boundary layer, and χ is the ratio of the universal gas constant and the heat capacity at constant pressure, equal to 0.25 for Mars. The PBL top is calculated inside MarsWRF's boundary layer scheme as the location where the bulk Richardson number falls to a critical value of 0.5 (Hong & Pan, 1996).
In MarsWRF, the sensible heat flux is given by: where C d is a drag coefficient that depends on the stability of the near-surface atmosphere, ρ is air density, u * is drag velocity, and T surf and T air are the surface and 1.5-m layer air temperatures, respectively. Both the sensible heat flux and u * are calculated inside MarsWRF's surface layer scheme (Zhang & Anthes, 1982), which uses Monin-Obukhov similarity and accounts for four stability categories: stable, mechanically induced turbulence, unstable forced convection, and unstable free convection.

Vortex and Dust Devil Predictions Across Gale Crater Using MarsWRF
This section presents DDA predictions across Gale Crater and explores how the DDA and its contributing terms are affected by surface properties and the circulation. Sections 4.1 and 4.2 focus on local (i.e., southern) summer solstice, Ls = 270°, the time of peak solar forcing and greatest vortex activity in Gale Crater, with seasonal variations explored in section 4.3. Section 5 then compares MSL observations of vortex pressure drops with predicted DDA along the rover traverse over 3 Mars years, using the context provided in section 4 to interpret those results. Figure 2 shows the predicted DDA across the crater at Ls~270°over eight time periods. Predicted DDA outside these periods is negligible (below 0.02). In terms of temporal variation, the predicted DDA peaks in the 12:00-14:00 period, with almost none predicted between 08:00 and 09:00 and only slightly more between 16:00 and 17:00. In terms of spatial variation, the lowest predicted DDA occurs in the crater trench to the N, NE, and NW of the MSL landing site and to the SE and S of Aeolis Mons, although there is also a band of low predicted DDA that runs eastward across the northern slopes of Aeolis Mons, starting roughly from the MSL landing site. The highest predicted DDA occurs largely on the NW, SW, and SE crater rim (as well as in some regions outside the crater) and in a north-south strip running through most of the western side of the 10.1029/2019JE006082 crater trench. There is also moderate DDA predicted over most of Aeolis Mons and in the SW quadrant of the trench.

Overview of Predicted DDA at Southern Summer Solstice
From an MSL perspective, based on these plots (see also the "zoomed-in" plots discussed in section 4.3), one might expect that looking from the landing site toward the slopes of Aeolis Mons or to the SW, between about local noon and 15:00, would have provided the best chance of observing dust devils in near-field MSL imaging. Conversely, looking anywhere between the NW and NE (i.e., toward the nearest portion of the crater rim), or toward the slopes of Aeolis Mons roughly due east of the landing site, would provide the worst chance, regardless of time of day. It is therefore unfortunate that MSL's early dust devil imaging focused on looking across the plains of the trench roughly to the north; this may have contributed to the lack of dust devils imaged in the first portion of the mission, as discussed in section 2.2. However, the changing position of MSL over the course of the mission (see Figure 1) also brought the rover higher up the NW slopes of Aeolis Mons, into the region with moderate predicted DDA, which would be expected to provide more near-field imaging of dust devils in several directions and more detections of in situ vortex pressure drops. Inside the crater, both are lowest over much of the trench (especially the regions picked out in section 4.1 as having low DDA) and are highest over much of Aeolis Mons and in a north-south strip at~137°E. However, the DDA minimum that stretches east across the northern slopes of Aeolis Mons is directly linked to the strong minimum in thermodynamic efficiency there. Note that the sensible heat flux peaks in the 12:00-13:00 plot, whereas the thermodynamic efficiency peaks in the 14:00-15:00 plot. This is because the former has a faster response to changes in solar heating, as it largely involves the surface heating up. However, the thermodynamic efficiency relates to the PBL depth, which takes some time to increase Journal of Geophysical Research: Planets following increased surface warming. Hence there is typically a delay of about 2 hours between the time of peak forcing (around local noon) and peak PBL depth. The different timings are discussed further in section 5.

Vertical Thermodynamic Efficiency and Its Contributors
The top two rows of Figure 4 show the vertical thermodynamic efficiency and its main contributor, the pressure thickness of the PBL P surf -P top (see equation (3)), and demonstrate how closely the spatiotemporal variation in the former follows that of the latter. Overall, in this season, the pressure thickness of the convective boundary layer within the crater is always greatest in the trench and grows significantly as the day progresses, with the thinnest boundary layer over Aeolis Mons and increasing only moderately over the times shown. Boundary layer depth is more often presented in height coordinates (bottom row of Figure 4), which displays a very similar pattern. Interestingly, this is almost the reverse of the spatial distribution noted in, e.g., Tyler and Barnes (2015), which focused on an earlier season: Ls~150°, the time of MSL's landing. While some observations have shown that PBL depth on Mars is larger where the pressure is lower, such as over high topography (e.g., Hinson et al., 2008Hinson et al., , 2019, the nature of the southern summer circulation at Gale Crater produces an inversion of this result. The cause of the spatial variation in PBL depth and its dependence on season is discussed further in section 4.3. Figure 5 shows the main contributors to the surface sensible heat fluxnear-surface air density ρ, drag velocity u * , and the surface-to-air temperature difference T surf -T air (see equation (4))in three daytime periods. Sensible heat flux is also influenced by differences in the drag coefficient C d , which varies with the stability of the boundary layer, but the additional variations introduced are less than a few percent. From a comparison of all four rows of Figure 5, the pattern of predicted surface sensible heat flux is most influenced by the variation of the surface-to-air temperature difference, followed by that of u * . Near-surface air density varies spatially by less than a factor of two over the whole domain in the 10:00-11:00 period and by even less as the day progresses and hence has little impact on the pattern of sensible heat flux via its presence as a multiplier in equation (4).

Surface Sensible Heat Flux and Its Contributors
On the contrary, there appears to be a strong negative correlation between the pattern of near-surface air density and sensible heat flux, but this is largely due to the crater topography which strongly influences both. Although surface temperature strongly controls near-surface air temperatures via radiative heating and sensible heat flux, daytime air temperatures measured at the same height, x, above the local surface will be cooler for an air parcel at height x above the slopes of Aeolis Mons than for one at height x above the crater trench. The air parcel above the higher ground is cooler due to horizontal (along constant geopotential surfaces) mixing with neighboring air parcels that are at higher altitudes above their own local surface (and are thus colder due to the generally negative daytime lapse rates in the real atmosphere). This results in larger T surf -T air over Aeolis Mons than over the crater trench and hence closely links surface topography to the DDA.
Other surface properties, such as thermal inertia and albedo, also have an impact on the pattern of T surf -T air shown in the bottom row of Figure 5. The top row of Figure 6 shows the albedo and thermal inertia maps used by MarsWRF in the domain 5 nest. Aeolis Mons has the lowest thermal inertias inside the crater, with the next lowest being in a thin strip running approximately north-south at~137°E. Intermediate values occur over most of the southern trench, with the highest values of thermal inertia over most of the northern trench and in an area just to the SE of Aeolis Mons. The surface albedo is also low in a north-south strip at 137°E and~138.5°E, as well as over the southern half of Gale Crater. Figure 7 compares results from the original simulation with one using uniform albedo and thermal inertia across domain 5. The variable surface properties in the original simulation appear responsible for (i) the north-south strip of high surface temperatures and high surface-to-air temperature differences at~137°E, due to lower albedo and thermal inertia causing the surface to warm more rapidly after sunrise, and for (ii) slightly reducing the surface-to-air temperature difference in the northern crater trench, due to high albedos and thermal inertias that result in lower daytime surface temperatures. Conversely, high albedos and low thermal inertias over Aeolis Mons appear to largely cancel out and to have little net effect on the overall pattern of predicted DDA. Overall, then, the bulk of the pattern of surface-to-air temperature difference inside Gale Crater appears largely insensitive to the albedo and thermal inertia map. In particular, Aeolis Mons has the greatest values of T surf -T air , with smaller values occurring in the crater trench, regardless of whether varying or uniform surface thermal properties are used. Thus the topographic control on air temperature, described above, appears to dominate over the influence of other surface properties on surface temperature.
The daytime pattern of u * is a result of constructive and destructive interference between (i) local daytime slope flows inside the crater, up the slopes of Aeolis Mons and the crater rim; (ii) regional daytime slope flows across the dichotomy boundary where Gale Crater sits, up the slope from approximately N to S (or NNE to SSW); and (iii) the global circulation at this time of year, when the Hadley circulation is strongest and consists of upwelling in the southern hemisphere, downwelling in the northern hemisphere, and a return flow from N to S across the equator. The resulting early afternoon circulation for this region and for Gale Crater in particular is shown in the bottom row of Figure 6, in the form of wind vectors at 1.5 m above the surface in domains 4 and 5 of the model. Regionally, northerly (i.e., from the north) winds dominate, due to the global and regional daytime flows being in the same direction in this season. There are some signs of the flow diverting to go around Gale Crater, but it also breaks into the crater over the NW rim, penetrating partway across the trench there. By contrast, local daytime upslope winds over much of the northern

10.1029/2019JE006082
Journal of Geophysical Research: Planets rim slopes appear to oppose the entry of wind across the rest of the northern rim and into the northern trench. Local daytime upslope flows on the steeper, north side of Aeolis Mons (which sticks up high above the altitude of the northern rim) reinforce the larger-scale northerly flow as it passes over the middle and higher slopes, with upslope winds on the gentler southern slopes too weak to oppose it except in one small region. Comparing to the third row of Figure 5, it is clear that the northwesterlies by the NW crater rim, and strong northerlies over much of Aeolis Mons, dominate the pattern of u * inside the crater.
Hence overall, it appears that the change in sensible heat flux along the rover traverse (from the landing site in the crater trench to Aeolis Mons' slopes) is predominantly a result of Gale Crater's topography, which produces both (a) air density variations driving surface-to-air temperature differences that peak over Aeolis Mons and (b) slope winds that also peak over much of Aeolis Mons, with the effect of albedo and thermal inertia maps being only secondary.

Cause of the Spatiotemporal Behavior of DDA in Southern Summer in Gale Crater
In summary, during southern summer in Gale Crater, the spatial pattern of sensible heat flux is controlled largely by thermal processes and wind patterns driven by the crater topography (and secondarily by albedo and thermal inertia variations), while the pattern of thermodynamic efficiency is caused by the suppression of the PBL over Aeolis Mons compared to the crater trench in this season. However, the latter is predicted to vary considerably with time of year, as explored further in section 4.3. Overall, the variation in sensible heat

10.1029/2019JE006082
Journal of Geophysical Research: Planets flux largely dictates the pattern of DDA, with the thermodynamic efficiency acting to diminish the peak over Aeolis Mons and weaken the minimum in the crater trench.
As described at the top of section 4.2, the sensible heat flux peaks closer to local noon, while thermodynamic efficiency peaks after~14:00. The net result is a single peak in predicted DDA between~noon and 14:00 in this season. In some other seasons, however, this later peak in thermodynamic efficiency, plus a late increase in sensible heat flux due to strong afternoon winds, results in a double-peak structure being predicted (see section 5.2.1).

Seasonal Variations in the Pattern of Predicted DDA in Gale Crater
This section examines the predicted seasonal and spatial variations across the NW quadrant of Gale Crater, which contains the MSL traverse. Figure 8 zooms in on that portion of domain 5 to examine in more detail the spatial variation at four times of year: Ls = 0°, 90°, 180°and 270°. Sensible heat flux and vertical thermodynamic efficiency peak at different times of sol (see discussion early in section 4.2) and hence are shown in Figure 8 at the times when this peak typically occurs (12:00-13:00 and 14:00-15:00, respectively), with DDA shown at both times of sol. Figures S2, S3, and S4 additionally show results at three times of sol for, respectively, DDA, sensible heat flux, and vertical thermodynamic efficiency.
As expected, the season with the weakest solar heating, local winter (Ls = 90°), has the lowest predicted DDA (left column), peaking near the top of Aeolis Mons at both times shown. The largest DDA over this region are predicted in local summer (Ls = 270°), with intermediate DDA at the two equinoxes (Ls = 0 and 180°). However, the peak values in summer occur close to the western edge of the region, at least 10 km from the rover traverse, out of the range of dust devil imaging and much too far away to affect the REMS pressure signal. These peaks are clearly tied to increases in the sensible heat flux there (third column), as discussed in section 4.2.
In the region closer to the rover traverse, sensible heat flux is generally stronger in all seasons over Aeolis Mons than in the crater trench. The pattern of surface-to-air temperature difference is very similar in all seasons, as it is largely determined by the crater topography and other surface properties. This means that seasonal variations in the pattern of sensible heat flux are largely dictated by seasonal variations in the pattern of drag velocity and hence winds. Around summer solstice, for example, the region where MSL has begun climbing the slopes of Aeolis Mons has particularly weak wind speeds (see third row of Figure 5, where a "tongue" of lower wind speeds extends from the trench up the slopes). Thus, at least in this simulation, u * is predicted to increase along MSL's route less than it would have had the rover climbed, e.g., a more western slope. Section 5 further discusses the relative importance of the wind field in capturing observed variations in vortex activity due to MSL's changing position.
As discussed in section 4.2.1, in summer the daytime PBL depth, hence vertical thermodynamic efficiency, is predicted to be smaller over the Aeolis Mons slopes than in the crater trench. However, the opposite is largely true over the three other seasons shown in Figure 8. In fact, the daytime PBL depth is predicted to peak over Aeolis Mons from Ls~330 to 180°in the 12:00-13:00 period (see Figure S4) and from Ls~30 to 150°in the 14:00-15:00 period. This is consistent with a number of studies carried out for MSL's landing season (Ls~150°), which predicted that the PBL depth would be larger over Aeolis Mons than in the trench (e.g., Tyler & Barnes, 2013). The change in pattern with season is due to whether the crater circulation is primarily driven by regional or local (crater interior) flows. Around local winter solstice (Ls = 90°), for example, the winds measured by MSL were strongly linked to local upslope (downslope) flows onto (off) Aeolis Mons (Newman et al., 2017;Richardson & Newman, 2018) with no clear regional influence. Strong, symmetric daytime flows up these slopes would be more likely to produce a higher daytime PBL depth over Aeolis Mons. Around local summer solstice (Ls = 270°), however, available MSL wind and aeolian observations (Baker et al., 2018;Viúdez-Moreiras et al., 2019) and modeling (e.g., Rafkin et al., 2016) suggest that strong regional and larger-scale flows penetrating into the crater from the north have more importance. As a result, the circulations inside and outside the craterespecially its northern halfare far more connected. This in turn results in an increased PBL depth on the northern rim and trench in summer compared to the usually "suppressed" crater interior boundary layer (e.g., Rafkin et al., 2016) but conversely disrupts the upslope winds that were responsible for increasing the PBL depth over Aeolis Mons in winter.

Journal of Geophysical Research: Planets
Overall, if MSL continues to scale the slopes of Aeolis Mons for another 3 years, Figure 8 suggests that increasing numbers of vortex pressure drops would be observed in most seasons. In local summer, however, this will be due to the increased sensible heat flux, as the vertical thermodynamic efficiency in that season is expected to decrease quite strongly as MSL climbs the slopes further.

Direct Comparison Between MarsWRF Predictions and MSL Observations
A quantitative prediction of the vortex activity observed by MSL can be made by sampling the predicted DDA at the rover location for each Ls and LTST hour. This can then be compared directly with observations of vortices and dust devils to assess the model's performance andif there is agreementto infer reasons for the observed seasonal and interannual changes. Section 5.1 provides the statistics of vortex pressure drops over the first 1980 sols (nearly 3 Mars years) of the mission, extracted using the method described in section 2.4. Because MSL landed at~4.6°S at Ls~150°, shortly before local (i.e., southern) spring equinox (Ls = 180°), each Mars year (YRs 1-3) is broken into four wide periods centered on local (i) spring equinox (Ls = 135-225°), (ii) summer solstice (Ls = 225-315°), (iii) fall equinox (Ls = 315-45°), and (iv) winter solstice (Ls = 45-135°). Section 5.2 compares these statistics with those obtained from MarsWRF output and demonstrates that the model predicts many aspects of the observed vortex activity. Section 5.3 examines the relationship between DDA and vortex numbers as a function of vortex strength, while section 5.4 examines the impact of simulation setup on results.

The Statistics of Observed Pressure Drops in 12 Periods Spanning Sols 1 to 1980
The stacked histograms in Figure 9 show the inferred number of pressure drops per hour for each hour between 07:00 and 17:00 LTST, for four seasonal time periods in the first 3 years of MSL's mission. These were obtained by normalizing the number of observed pressure drops by the number of measurements in each period. The pressure drops are also binned according to the magnitude of the detected pressure drop of each event, using four logarithmic magnitude bins that cover the range from the detection threshold of 0.6 Pa to the maximum vortex pressure drop magnitude detected over this period of just under 5.6 Pa. The results for the first 2 years are consistent with those presented in previously published work (see section 2.3) and are also consistent with results for all 3 years obtained using a method very similar to that of Ordonez-Exteberria et al. (2018) (see Text S1 and Figures S5 and S6).
In most Mars years, the largest number of vortex pressure drops occurs in the period surrounding summer solstice and the smallest in the period surrounding winter solstice, with one main exception. In YR1, the period around summer solstice was relatively low in both number and magnitude of pressure drops compared to the period surrounding spring equinox. As suggested to explain similar results found in Kahanpää et al. (2016), this may have been due to the rover sitting at one location from MSL sol 166 to 272 (Ls~250-317°), which is very similar to the YR1 summer solstice period (defined as Ls = 225-315°) shown in the second plot in the top row of Figure 9. This location ("Yellowknife Bay") was measured to have a much higher thermal inertia than nearby regions (Martínez et al., 2014), which would have produced a slower rise in surface temperature during the day, reduced the daytime sensible heat flux, and very likely resulted in the reduced vortex activity observed by REMS. By contrast, the thermal inertia maps used in MarsWRF are derived from orbital data and do not capture such local variations; thus MarsWRF did not predict such a decrease in the DDA (see discussion in section 3.1).
Two other periods also come close to being exceptions. In YR3, a comparable number of pressure drops were measured around winter solstice as around the previous spring equinox. However, this is most probably a result of the strong increase in vortex activity as the rover climbed rapidly higher on the slopes; i.e., there were likely far more pressure drops recorded in the period around the following spring (not considered in this study). Second, the period around fall equinox in YR2 has a smaller number of pressure drops overall than the period around that summer solstice, but there is a greater peak in number of drops between 11:00 and 13:00. This is due to several sols when the atmosphere was anomalously convective around local noon, as also noted by Ordonez-Exteberria et al. (2018). Typically, the number of 1-hour measurement periods containing N vortices (shown in Figure S7) decreases logarithmically with N up to N~10 and then declines more slowly, but a "hump" around N~18 corresponds to such especially restless periods.
In most seasons and Mars years, the peak number of observed vortex pressure drops occurs between 11:00 and 14:00 LTST. It appears to skew earlier in this period around the equinoxes (first and third columns) and later around the solstices (second and fourth columns), although again the season spanning local summer in YR1 is anomalous here. Finally, there is a strong overall year-to-year increase in the number and magnitude of observed vortex pressure drops for most periods and times of day, with the greatest activity occurring when MSL sits highest on the slopes of Aeolis Mons.

10.1029/2019JE006082
Journal of Geophysical Research: Planets 5.2. MarsWRF-Predicted DDA at the Rover Location Figure 9 also shows MarsWRF predictions of DDA for the same time-of-day range and seasonal periods for which the REMS pressure drop data are shown. Each panel is produced by stepping through that seasonal period sol by sol, sampling the minute-by-minute DDA predictions in the model at MSL's location for that sol (interpolating between results for the two MarsWRF simulations that are closest to it in Ls), and then averaging the results at each time of day over all sols contained in the period. Comparing the observations and MarsWRF predictions shows that the general seasonal and diurnal variation is largely captured in the variation of DDA predicted using MarsWRF output. The general trend of increased vortices from year to year for most seasons and times of sol is also captured by the model; however, the observed number of vortex pressure drops increases faster than the predicted DDA as MSL climbs up the slope of Aeolis Mons. This is discussed further in section 5.3.
For both the observed number of pressure drops and predicted DDA, the peak generally occurs in the period containing summer solstice, is reduced in the periods around both equinoxes, and is smallest in the winter solstice periods. The exceptions are around YR1 summer solstice, when observed activity is low (see discussion in section 5.1), and around YR3 winter solstice, when the observed activity is comparable to that around the preceding spring equinox. The latter is also predicted by MarsWRF, however, and (as discussed in section 5.1) is likely due to the year-to-year increase as MSL climbs Aeolis Mons. Both peak DDA and number of vortices occur around the solstices in the 13:00 to 15:00 period (again, except for the summer solstice of YR1, when observations dip) and occur earlier around the equinoxes. However, the clear predicted "double-peak" structure in DDA around spring equinox, which is also visible around fall equinox, is generally not seen in the observed number of vortices or pressure drop magnitudes. By contrast, there is a hint of a double-peak structure in the observations around summer solstice (particularly in YR3, whereas in the model it is most noticeable in YR1).
In summary, there is general agreement but also some significant differences between the predicted DDA and observed vortex activity. However, given the several areas of agreement, it is possible to dig into the model results to examine what may be responsible for the observed vortex variability with season, time of day, and location on Mars.

Controls on the Seasonal and Year-to-Year Variation of DDA
As discussed above and shown in Figure 9, MarsWRF generally does a good job of capturing the observed seasonal variation of pressure drops. Figure 10 shows the diurnal variation of each relevant variable in MarsWRF for each seasonal period at the rover location and demonstrates the relative contribution of sensible heat flux and vertical thermodynamic efficiency. The former typically peaks shortly before or around noon, whereas the latter peaks between~14:30 and 16:00, shifting later as the season shifts from spring to summer to fall. These results reveal that the peak DDA occurs around summer solstice due largely to a stronger afternoon sensible heat flux than around the equinoxes, which in turn is due largely to a stronger surface-to-air temperature difference at this time of year. Finally, around winter solstice, the surface-to-air temperature difference and PBL depth are both far lower than in other seasons at most times of sol, with no afternoon peak in predicted PBL thickness, which results in the lowest DDA being predicted for this season relative to other times of year.
Figure 10 also demonstratesas discussed in section 4.3that an increase in sensible heat flux from year to year, not thermodynamic efficiency, is responsible in most seasons for the predicted increase in DDA as MSL climbs up the slope of Aeolis Mons. Only around fall (and, to a lesser extent, spring) equinox is an appreciable increase predicted in peak PBL depth (hence thermodynamic efficiency) in the afternoon from year to year, whereas peak sensible heat flux is predicted to increase significantlyand monotonicallyfrom year to year in every season. Differences between years are discussed further in section 5.3.  (1) and (4), energy is available for vortex formation (hence nonzero DDA values occur) when the sensible heat flux is positive, i.e., when the surface is warmer than the atmosphere above it. The model in fact predicts no DDA before 07:50 or after 17:00 in any season. These predicted time periods of nonzero DDA contain all of the observed daytime pressure drops exceeding 0.6 Pa (see Figure 9).
Figure 10 reveals that the predicted double-peak structure in some seasons is somewhat due to the different timings of the peak in surface-to-air temperature difference and thermodynamic efficiency but primarily due to the double-peak structure in drag velocity in some seasons. Surface-to-air temperature difference peaks either in the hour before noon (around both equinoxes in all years and around summer solstice in YR3) or in the hour after noon, whereas thermodynamic efficiency peaks at~15:00 around spring equinox and summer solstice, up to an hour later around fall equinox, and is flat from~10:30-15:30 around winter solstice. In addition, drag velocity peaks strongly in the late morning and mid-afternoon, with a dip in the early afternoon, around both equinoxes, particularly spring. The net effect is both a late morning and early to midafternoon peak in the predicted DDA in many years and seasons, especially around spring and to a lesser extent fall equinox. Figure 9 allows the daily timings of observed peaks in vortex activity to be compared with those predicted by MarsWRF as a function of season. The best agreement occurs around winter solstice, when both the model and REMS data show a nearly symmetric distribution about the noon-13:00 period, especially in YR3. Around summer solstice, by contrast, the model predicts a less symmetric distribution, with DDA peaking between~12:00 and 15:00 and a hint of a secondary peak between~10:00 and 11:00 that fades from year to year and is gone by YR3. This is visible in the observations too, although in YR3, the morning peak is still present, and in YR1, the number of pressure drops >0.6 Pa peaks in the late morning rather than afternoon. However, as noted in section 5.1, YR 1 suffered from relatively few REMS measurements and the rover sitting in a location with an anomalously high thermal inertia for over a hundred sols, so the data set is hampered by low statistics and possibly biased by local surface properties.
Around spring equinox, the morning and afternoon peaks are predicted to be rather similar in magnitude, with a distinct dip between them shortly after noon. However, the observations show only a hint of a twopeak structure, with a secondary maximum in the afternoon in YR2 but no sign of this at all in YRs 1 or 3. Around fall equinox, two peaks are again predicted, but with a smaller dip between them which shrinks further with position up the slope (i.e., with year). As in spring, there is barely a hint of two-peak structure in the observations, and this disappears by YR3. The net effect around both equinoxes is an overprediction of the afternoon DDA, primarily due to an apparently erroneous peak in u * in the afternoon (see yellow lines in Figure 10). Interestingly, despite using output from a different model, run at far lower resolution (5°x5°grid spacing), Chapman et al. (2017) predict a similar double-peak structure in DDA when averaging over the whole of YR1. However, given the strong sensitivity of winds to the crater topography (which is not resolved in their model) and the fact that the double-peak behavior changes when MarsWRF is set up differently (see section 5.4), this agreement is likely coincidental.

Seasonal and Year-to-Year Changes as a Function of Vortex Strength
This section now ignores diurnal variability and focuses on the total number of pressure drops per sol. This enables easier comparison with MarsWRF predictions of seasonal and interannual variations but also provides improved statistics for looking at the match as a function of pressure drop magnitude. According to R98, the maximum tangential wind velocity around a convective vortex is proportional to the square root of its central pressure drop. The central pressure drop can hence be used as a measure of vortex strength.
As the vortices do not pass directly over MSL, the detected pressure drop magnitudes are smaller than the actual central pressure drops of the vortices. Nevertheless, the temporal distribution of the detected pressure drop magnitudes reflects the distribution of the vortex strengths. Figure 11a shows the total inferred number of vortex pressure drops >0.6 Pa per sol and MarsWRF predictions of mean DDA as a function of Ls over the 3 years used in this work, again binned into four broad seasons in each year. This clearly demonstrates the degree to which the DDA underpredicts the year-to-year increase in the number of vortices observed by MSL. However, if the observations and the DDA are plotted using a minimum DDA value of 0.17, the match is much better (Figure 11b). This is especially true if the summer solstice period in YR1 is discounted as anomalous, as discussed in section 5.1. This suggests that a threshold DDA may apply below which no vortices occur (or at least none strong enough to be detected by the methods used here) and that DDA may be better able to predict the total number of vortices if the existence of such a threshold is taken into consideration. To examine this in a quantitative manner, however, such a DDA threshold should be applied to each short time period (e.g. every hour of every sol), rather than being applied to the mean DDA value over a period lasting >100 sols (as is effectively shown here). Figure 11 also shows the comparison between DDA and the inferred number of vortex pressure drops >1.05 Pa (second row), >1.83 Pa (third row), and > 3.2 Pa (fourth row). The match is increasingly poor as the pressure drop magnitudehence strength of the vortexincreases. This suggests that a higher threshold DDA may apply to stronger vortices and/or that the relationship between vortex strength and DDA may not be linear. Comparing Figures 11d, f, and h with Figure 11b (which is dominated by the most abundant, i.e., weakest, vortices) suggests that increasingly high thresholds (somewhat equivalent to starting the y-axis at a higher value in Figure 11d, higher again in Figures 11f, and still higher in Figure 11h) might improve the match to observations of increasingly strong vortices. If one assumes that a stronger vortex (with a deeper central pressure drop and hence higher tangential wind speeds) is required to raise dust from the surface, Figure 11. REMS total daily inferred pressure drops, separated according to the drop magnitude (histograms), compared to the mean DDA over the 9 hours during which convective vortices are observed (08:00-17:00) at MSL's location predicted by MarsWRF (dark green lines), as a function of Ls over 3 Mars years, averaging over broad (90°-wide) Ls bands. In the left column, the DDA axis starts at 0.0, whereas in the right, it starts at 0.17.

10.1029/2019JE006082
Journal of Geophysical Research: Planets one might expect the temporal variation of imaged dust devils to match that of larger pressure drops (e.g., Figures 11e or 11g, rather than 11a) and thus to have a threshold-dependent and/or nonlinear relationship with DDA also. This would be consistent with the findings of Kahre et al. (2006) that MER Spirit observed zero dust devils in some seasons, whereas the predicted DDA never dropped to zero at that location. The relationship between predicted DDA and dust devils imaged by surface cameras will be explored in future work.

Impact of Simulation Setup on Results
Given some of the mismatches between the observed and modeled vortex activity, one may ask whether the MarsWRF predictions or the theory of R98 are responsible. As discussed in section 5.3, the present study does suggest that a threshold DDA may be needed to connect R98 to pressure drop statistics, but this study is not well suited to determining the accuracy of R98. A better approach would be to perform terrestrial experiments in which the local and regional atmospheric state can be measured and predicted with much higher accuracy than is available for Mars, and compared to more complete and extensive measurements of vortices and dust devils.
This section therefore focuses on potential MarsWRF errors. As noted in section 3.1, all results shown above used output from a simulation that has previously been shown to represent winds and aeolian features inside the Gale Crater reasonably well. However, it does so by using a vertical grid with rather poor vertical resolution in the lower boundary layer (Newman et al., 2017). It thus makes sense to begin investigating the source of the mismatchesespecially the timing of peak vortex activityby examining the results using a different vertical grid in the model, specifically, vertical grid A of Newman et al. (2017) which has three additional model layers centered below~150 m. Figure 12 shows the match to observed pressure drops as a function of time of sol and season using vertical grid A and may be compared directly to Figure 9 showing grid B results. The main differences are in the times of peak DDA, with Figure 12 not predicting the erroneous Figure 12. As in Figure 9 but now using output from a MarsWRF simulation using vertical grid A (Newman et al., 2017). See text for more details.

10.1029/2019JE006082
Journal of Geophysical Research: Planets mid-afternoon peak around the equinoxes seen in Figure 9. In the original simulation, as shown by Figure 10 , the late morning peak is largely due to a peak in u * (hence sensible heat flux) at this time, while the midafternoon peak is a combination of a second peak in u * and the growing PBL depth (hence vertical thermodynamic efficiency). The new simulation does not have this strong double-peak structure in u * and hence does not predict the erroneous mid-afternoon peak (see Figure S8). However, the new simulation predicts peak DDA at equinox between 13:00 and 14:00, which is later than the observed peak in number of vortices.
Overall, there is no clear improvement in the match to diurnal variability when using the new simulation and no improvement to seasonal or year-to-year variability at all.
The next step beyond what has been done here would be to use dust opacity maps based on regional (e.g., Montabone et al., 2015) and local (e.g., Guzewich et al., 2018) observations of atmospheric dust content, in order to vary MarsWRF's atmospheric dust distribution more realistically. Work on this is underway, including modeling the impact of the global and regional dust storms in MY 34, the latter of which was also experienced by the InSight lander and recorded in its more sensitive measurements of pressure. The results of these MarsWRF simulations will be used to study the impact of dust opacity changes on the predicted DDA, which will be compared with the observed vortex pressure drops and dust devil statistics at both the MSL and InSight landing sites during the storm(s).
Ideally, future work would also compare observed vortex pressure drops with the explicit simulation of vortices in a Large Eddy Simulation (LES) of Gale Crater. While MarsWRF may also be run in this way, this is a difficult experiment to perform due to the need for the LES to include the whole of Gale Crater to properly capture the local topographic influences, and yet be run at~50 m grid spacing in order to capture the majority of turbulent structures and resolve larger convective vortices. In practice, an area four times the size of the crater would be needed to avoid boundary condition issues and edge effects; hence this would require~6,000 grid points in longitude and in latitude, a very expensive simulation. In addition, the LES would need to be driven with time-varying regional wind profiles taken from a lower resolution MarsWRF nested mesoscale simulation. Finally, even at such a high resolution, vortices smaller than~200 m diameter would not be properly resolved. However, the benefits of being able to model the PBL more explicitly than in a mesoscale simulation (in which a PBL subgrid-scale parameterization is required) and to compare MSL observations directly with model predictions (e.g., the distribution of vortex pressure drop magnitudes, timings, etc.) make this idea very appealing.

Conclusions
The thermodynamic theory of R98 is combined with output from the MarsWRF atmospheric model, run at a grid spacing of~1.4 km, to predict the variation of dust devil (or equivalently, convective vortex) activity, DDA, for Gale Crater, Mars. This is then compared with the statistics of daytime surface pressure drops >0.6 Pa measured by REMS on the MSL rover over 3 Mars years, which are associated with the passage of convective vortices. In the R98 theory, DDA is proportional to the sensible heat flux, which increases with the surface-to-air temperature difference and drag velocity, and to the vertical thermodynamic efficiency of the PBL, which increases with the pressure thickness of the PBL.
MSL landed at~4.5°S at Ls~150°, so each Mars year is broken into four wide periods centered on southern (i) spring equinox (Ls = 135-225°), (ii) summer solstice (Ls = 225-315°), (iii) fall equinox (Ls = 315-45°), and (iv) winter solstice (Ls = 45-135°). Except for YR1, when there were large data gaps and the rover remained in an anomalously high thermal inertia region for over a hundred sols, the greatest number of inferred vortex occurrences and largest pressure drop magnitudes were observed in the period around summer solstice. MarsWRF output also predicts peak DDA around summer solstice, which is attributed to both high surface-to-air temperature differences from mid-morning through mid-afternoon and large boundary layer thicknesses in the mid-afternoon. In terms of diurnal variations for this season, MarsWRF output predicts peak DDA between noon and~15:00, with a small secondary peak between 10:00 and 11:00 in YR1 only, whereas REMS pressure drops suggest more of a two-peak structure in all years.
Significant numbers of vortices are also observed and predicted around both equinoxes, but while observations show a peak in vortex activity between 11:00 and 13:00, MarsWRF predicts both peak DDA between 11:00 and 12:00 and a second peak in the midafternoon after~14:00, especially around spring. In the model,

10.1029/2019JE006082
Journal of Geophysical Research: Planets the late morning peak is largely due to a peak in u * (hence sensible heat flux) at this time, while the midafternoon peak is a combination of a second peak in u * and the growing PBL depth (hence vertical thermodynamic efficiency). Another MarsWRF simulation without the strong double-peak structure in u * does not predict the erroneous midafternoon peak; however, the predicted DDA still peaks around equinox between 13:00 and 14:00, which is later than the observed peak in number of vortices.
Finally, the period around winter solstice has the fewest dust devils observed and predicted in YRs 1 and 2, with MarsWRF output showing both the smallest sensible heat flux (largely due to a smaller surface-to-air temperature difference) and PBL depth, with no midafternoon peak in the latter as occurs in all other seasons. The winter solstice period has the most symmetric diurnal distribution of pressure drops around noon-13:00, as predicted by the model which shows this being due to the relative flatness of the PBL depth with local time.
Observations show an increase in number and also magnitude of vortex pressure drops from year to year in all seasons, as the rover made its way out of the trench and up the lower slopes of Aeolis Mons, with a yearto-year increase also found in the predicted DDA. The rate of the year-to-year increase in number of vortex pressure drops >0.6 Pa is underpredicted by the model if vortex numbers are assumed to be simply proportional to DDA. However, the match may be far better if a threshold DDA is chosen below which no pressure drops are expected (or are too small in magnitude to be detected by the methods used here). This result, as well as the greater mismatch between DDA variation and variation in number of larger magnitude pressure drops, suggests that the relationship between DDA and vortex pressure drop statistics is likely thresholddependent and/or nonlinear.
The model does not generally attribute the year-on-year increase in DDA to an increase in PBL depth over the mound (as in some previous work). Indeed, a detailed study of how the PBL depth is predicted to change with rover position shows it only increasing slightly around the equinoxes as MSL climbs the slopes, and basically static from year to year in the other seasons. Instead, the greater DDA and vortex occurrence as the rover climbed Aeolis Mons are primarily attributed to an increase in surface-to-air temperature and u * as MSL moved higher on the NW slope. Both of these are found to be primarily a function of the crater topography. The surface-to-air temperature difference increases higher up the slope due to cooler nearsurface daytime air temperatures caused by mixing with cooler adjacent high-altitude air, with an additional small effect due to lower regional thermal inertias over the slopes. In the case of u * , the rover experiences stronger peak daytime upslope winds as it climbs the slopes.
The spatial variation in predicted DDA may also explain why no dust devils were imaged for about 2 years at the start of the mission, when imaging campaigns were looking away from the mound toward a region of lower predicted DDA. Dust availability may also have been a factor, but the results of this study suggest that seasonal DDA predictions would be very helpful for directing imaging campaigns for future landed missions, particularly those landing in regions with strong topography or other major surface property contrasts.

Erratum
In the originally published version of this article, several instances of text were incorrectly typeset in the key points and in the abstract. The errors have since been corrected and this version may be considered the authoritative version of record.