Fully Interactive and Refined Resolution Simulations of the Martian Dust Cycle by the MarsWRF Model

The MarsWRF model is set up with fully interactive dust at 5° × 5° and 2° × 2 resolution. The latter allows for a better representation of topography and other surface properties. An infinite reservoir of surface dust is assumed for both resolutions. For 5° × 5°, surface dust lifting by wind stress takes place over broad areas, occurring in about 20% of the model's grid cells. For 2° × 2°, it is more spatially restricted, occurring in less than 5% of the grid cells, and somewhat reminiscent of the corridors Acidalia‐Chryse, Utopia‐Isidis, and Arcadia‐West of Tharsis. The onset times of major dust storms—large regional storms or global dust storm events (GDEs)—do not exhibit much interannual variability, typically occurring at around Ls 260°. However, their magnitude does show significant interannual variability—with only small regional storms in some years, large regional storms in others, and some years with GDEs—owing to the interaction between major dust lifting regions at low latitudes. The latter is consistent with observed GDEs having several active dust lifting centers. The agreement between themodel's surface dust distribution and observation‐based dust cover index maps is potentially better for 2° × 2°. For the latter, there is also significant surface dust lifting by wind stress in the aphelion season that is largely confined to the Hellas basin. It has a recurring time pattern of 2–7 sols, possibly resulting from the interaction between midlatitude baroclinic systems and local downslope flows. Plain Language Summary Mars General Circulation Models (MGCMs) simulate the Mars climate and atmosphere for many Martian Years (MYs). A challenge is to configure such models to produce Mars global dust storm events (GDEs) in few but not all MYs. That is because GDEs are known to occur once every few MYs, on average. We set up the MarsWRF MGCM in this way using “interactive dust,” meaning the model freely lifts, transports, and deposits surface dust (assuming an inexhaustible amount of available surface dust). We use different horizontal model grid point resolutions and compare their results in terms of dust storm source regions and changes in surface dust loading. For the high‐resolution experiment, we find that GDEs are likely to develop if regional dust storm activity around two equatorial source regions on the planet, namely, south of Chryse Planitia and in the northern Hellas basin, combine with one another. The latter is consistent with knowing from observations that GDEs have several active dust lifting centers. Also, we find that the model's surface dust distribution in the high‐resolution experiment agrees potentially better with observation‐based dust cover maps.


Introduction
The dust cycle is a crucial component of the Martian climate (e.g., Zurek et al., 1992). The reason is that dust greatly affects the thermal and dynamical state of the Martian atmosphere. Dust particles are radiatively active, absorbing solar radiation and absorbing/emitting infrared radiation. Thus, they heat up/cool down the atmosphere and drive dynamical processes (e.g., Gierasch & Goody, 1972;Pollack et al., 1979). The effects of the dust on climate may be enhanced by feedback mechanisms like thermal gradients at the edge of regional dust storms, which trigger stronger winds, lifting even more dust. Therefore, a profound knowledge of the dust cycle is of key importance for understanding the behavior of the Martian atmosphere.
The dust cycle on Mars exhibits strong variability on annual and interannual time scales (e.g., Newman & Richardson, 2015). The annual variability is largely a result of the Martian hemispheric dichotomy in topography (Richardson & Wilson, 2002). Around perihelion, when it is spring/summer in the southern hemisphere, the atmosphere is relatively warm and dusty, with the higher elevations of the southern hemisphere enhancing the global overturning circulation, thus leading to increased dust lifting (Richardson & Wilson, 2002). Around aphelion, on the other hand, lower amounts of dust are present in the atmosphere, and while regional and planet-encircling dust storms are unlikely, local dust storms may still occur (Zurek et al., 1992). While any Martian year (MY; the convention of Clancy et al., 2000, is used in this work) has local and regional dust storms, global dust storm events (GDEs) occur, on average, once every few MYs (e.g., Haberle et al., 1982;. As an example, the last GDE took place in the middle of 2018. It was monitored from orbit by the Mars Reconnaissance Orbiter (Shirley et al., 2020) and from the surface by the Curiosity rover (Guzewich et al., 2019) (and initially also by the Opportunity rover). As a direct consequence of the dust storm, the Opportunity rover mission came to an end (Voosen, 2019), after having been active since 2004.
Mars General Circulation Models (MGCMs) provide valuable insight into the Martian climate by filling the spatial and temporal gaps in observational coverage and providing predictions for atmospheric fields that have not been observed directly (e.g., winds). MGCMs may be constrained by observational data to various degrees: through data assimilation techniques where the most guidance is provided to the MGCMs (e.g., Greybush et al., 2012;Guzewich et al., 2013), with simulations where only the total column amount of dust is constrained (prescribed dust) or in simulations where the model freely lifts and advects the dust in a self-consistent manner (interactive dust). As summarized in Montabone et al. (2015), analytical expressions for the horizontal and vertical dust distribution in MGCMs have been increasingly replaced with observation-based data from the Thermal Emission Spectrometer (TES; Christensen et al., 2001) and other remote sensing instruments. Either concurrent climatological data or prescribed dust are commonly employed. The Mars dust climatology is described comprehensively in Montabone et al. (2015). Prescribed dust includes scenarios for specific past MYs, scenarios representative of background dust (i.e., the hypothetical atmospheric dust content of a year with no major dust storm occurring) throughout the entire year, and dust storm year scenarios (e.g., Guzewich et al., 2013;Haberle et al., 2017;Natarajan et al., 2015). Such an approach, although computationally cheap, has a major drawback. It neither offers insight into the processes governing surface dust gain/loss nor allows for a long-term prediction of the dust cycle (Newman & Richardson, 2015).
Alternatives to the integration of the observed dust optical depth are MGCM simulations with interactive dust. In this set up, the surface dust lifting and deposition are self-consistently generated by the model. Newman and Richardson (2015) ran the Martian implementation of the Planet Weather Research and Forecasting (WRF; Skamarock et al., 2008) model, MarsWRF (Richardson et al., 2007), with interactive dust, considering both infinite and finite surface dust availability. A first set of simulations, run with infinite surface dust, resulted in shortcomings in the simulated onset regions, times, sizes, and durations of GDEs. In particular, the dust cycle was not realistic (e.g., just over a third of the lifted dust came from three source grid points, and the interannual variability was not totally in line with observations), and there was a lack of replenishment of surface dust for the primary source region. For example, in some dust source regions, there is a net loss of dust at a high rate of~1 cm every 30 MYs. Although surface dust in these simulations is assumed to be inexhaustible, this high rate of dust depletion is unlikely to be sustainable on real Mars. Subsequently, simulations with finite dust were proven to significantly improve the dust cycle characteristics, in particular the onset times of GDEs (early for L s ∼ 210-220°; late for L s ∼ 252°). More recent studies suggest that orbit-spin coupling, a theory linking Solar System dynamics to interannual variability in wind patterns on Mars, may be a potential mechanism for correctly predicting the occurrence of dust storms in some MYs and their absence in others (e.g., Newman et al., 2019;Shirley et al., 2019). Both of the latter studies are based on infinite surface dust. In Newman et al. (2019), the use of finite surface dust reservoirs is pointed out to have the potential to improve the onset times of the dust storms. For all three studies mentioned above, a model grid of 5.625°× 5°in longitude/latitude, with 20 vertical levels was used.
Specific cases of dust storms have been studied using MarsWRF at a similar resolution of 5°× 5°with 52 vertical layers and a model top at 80 km. Chow et al. (2018) ran the model with interactive dust, assuming an infinite dust source, to investigate the dynamics behind the southern hemisphere spring equinox dust storms that develop in the Hellas basin (Hellas storms). The authors concluded that those storms are triggered by strong nighttime downslope winds in the Hellas basin, arising from the large contrast in surface temperature between the southern edge of the basin and the CO 2 ice-covered south polar cap. Using the same configuration, Xiao et al. (2019) investigated the dynamical processes of dust lifting at the northern midlatitude region of Mars during the dust storm season (L s 180-360°). They reported three different modes of dust lifting activities: seasonal mode (time scale of 50-200 Sols), driven by the annual path of the Sun, synoptic mode (time scale of 3-7 Sols), associated with midlatitude eddies, and diurnal mode (time scale of 0.5-1 Sol), linked with the thermal tides and their modulation by the topography. Finally and in semiidealized experiments where the zonal asymmetry of the terrain poleward of 15°N was removed, it was concluded that the topography of the northern midlatitude region (40-50°N) was important for dust lifting activities. Both of these studies acknowledged that the coarse spatial resolution (5°) of the model grid was a limitation as the topography and, for Chow et al. (2018), the boundary of the ice cap edge were not sufficiently represented in the model. Fonseca et al. (2018) identified the use of a prescribed dust scenario as a possible cause for discrepancy between the diurnal cycle of meteorological variables observed by Curiosity's Rover Environmental Monitoring Station (Gomez-Elvira et al., 2012) and predicted in nested mesoscale simulations. Similarly, Fonseca et al. (2019) showed that the choice of prescribed dust scenarios is critical to reproduce the entry, descent, and landing (EDL) vertical profiles measured by past missions to the planet. In fact, when forced with dust scenarios specific to the particular MY of each EDL event, Natarajan et al. (2015) concluded that MarsWRF gives a more accurate simulation of the observed temperature and density profiles. These results highlight the need to properly represent the observed dust distribution and its spatial and temporal variability, something that is attempted by setting up the model with interactive dust lifting.
In this study, we present MarsWRF runs with interactive dust lifting and an infinite surface dust reservoir at a refined model grid of 2°× 2°to allow for a better representation of surface topography and other surface properties. This, in turn, results in changes in larger mesoscale wind features and associated surface dust lifting. The output of this simulation is compared with another simulation at 5°× 5°, using an otherwise identical experimental setup. This is motivated by the fact that 5°× 5°is a similar spatial resolution to that employed in previous studies (e.g., Chow et al., 2018;Newman & Richardson, 2015;Xiao et al., 2019). Toigo et al. (2012) refined the model configuration from a spatial resolution of 5°× 5°to 2°× 2°and down to 0.5°× 0.5°. Comparing the total amount of dust lifted from the surface during one MY of simulation, they found that the dust lifting in the 2°× 2°grid exceeded that of the 5°× 5°grid for threshold wind stresses greater than 0.014 Pa and by a factor of~2.3 for a threshold wind stress of 0.049 Pa. This is explained by the finer grid better representing the locations with high wind stress governed by topographic slopes. Our comparison between 2°× 2°and 5°× 5°thus offers an opportunity for evaluating the potential of increased horizontal model resolution for improving the representation of the dust cycle/GDEs. It is important to note that Toigo et al. (2012) did not model the dust cycle interactively as done in this article. Instead, they employed a prescribed dust scenario and calculated the global amount of dust lifted for different values of wind stress lifting threshold. (This was done without considering radiative dynamic feedback mechanisms, accordingly.) This paper is structured as follows. In section 2, the model configuration is provided and the sensitivity experiments for the tuning of the dust lifting parameters are described. In section 3, the final selection of the model's tunable dust parameters and the resulting dust storms are presented. In section 4, the long-term characteristics of the Martian dust cycle for the optimal model configuration are analyzed and discussed. The main findings are summarized in section 5.

Model Configuration and Sensitivity Experiments
This study is based on using the MarsWRF model version 3.3.1, set up in a 2°× 2°and 5°× 5°configuration without a water cycle and with the model topography given in Figure 1. As can be seen and, in particular, at 2°× 2°, important surface features, such as Hellas basin, Valles Marineris, Olympus Mons, the Tharsis volcanoes, Elysium Mons, and their associated topographic slopes, are reasonably captured. This relatively high resolution (and the performing of multiple model runs, as described in the following) is computationally affordable by assuming a quasi-infinite surface dust reservoir or, briefly, "infinite dust." This means that any change of surface dust is tracked, but the initial surface dust cover is set to be so large (1,000 kg m −2 ) that no region will ever be exhausted. The tracking of dust implies that the total mass of dust on the surface and in the atmosphere is conserved. The main advantage of infinite dust is the fact that the surface dust distribution does not require any spin-up. Note that the model spin-up may take up to hundreds of MYs in the case of a finite surface dust cover (Newman & Richardson, 2015). In the vertical, 45 levels are used, given in Table 1, with a finer resolution in the PBL (planetary boundary layer) region and a model top at 120 km.
The model setup is largely similar to the global grid of Fonseca et al. (2018Fonseca et al. ( , 2019. The topography is interpolated from the 1/64°Mars Orbiter Laser Altimeter (MOLA) data set (Smith et al., 2001), and MY26 albedo and TES apparent nightside thermal inertia maps are used (Putzig & Mellon, 2007). The simple CO 2 microphysics scheme employed in MarsWRF assumes that if any layer in the vertical column becomes supersaturated, CO 2 directly deposits onto the surface, while if the near-surface atmosphere over CO 2 ice becomes subsaturated, CO 2 sublimes back to the atmosphere. The total CO 2 inventory and polar cap properties have been adjusted so that the model generates a surface pressure cycle that closely matches the Viking Landers 1 and 2 observations (Guo et al., 2009). The PBL scheme used is the Medium Range Forecast (MRF; Hong &

Journal of Geophysical Research: Planets
Pan, 1996) PBL scheme. The radiation scheme employed is a kdistribution radiative transfer model (Mischna et al., 2012). At latitudes poleward of 70°, a zonal filter damps high-frequency waves, in order to mitigate the need for smaller time steps, given the decreasing zonal grid spacing with latitude (Richardson et al., 2007). In the top 4 km, Rayleigh damping is applied to the u, v, and w wind components and potential temperature on a time-scale of 5 s (Skamarock et al., 2008).
The model time step is 1 min, and the output is stored every sol at 0 UTC. (For producing Figures 6 and 7, some additional few-sol model runs were performed with saving the model output at every full hour.) At the beginning of each model run, a spin-up takes place. The model spin-up ensures that transient motions are damped and a stable circulation is initiated in terms of the representation of the CO 2 cycle and dynamics of the CO 2dominated Martian atmosphere. For reasons of model spin-up, it is common practice to discard the first 1-2 years of any model run from scientific analysis. In line with that, we discarded the first year of our 5°× 5°m odel run and denoted subsequent years as Year 1, Year 2, etc.
In the case of our 2°× 2°run, we decided to focus on the 10-year period between the 6th year, L s 180°, and the 16th year, L s 180°. This period is particularly interesting given that it includes practically everything from nonmajor dust storm years to GDEs, but no runaway dust storms. It is given by mid-Years 5 to 15 using the same denotation of years as introduced before. Reminiscent of this 10-year period, the model run had also started promisingly. Year 1 was a nonmajor dust storm year, Year 3 was at the boundary between a nonmajor dust storm and potential GDE year, and Year 2 had a GDE, albeit a relatively weak one, according to the dust storm classification scheme of Newman et al. (2019). In Year 4, however, this was followed by an almost runaway dust storm that extended beyond L s 0°and persisted until Year 5 L s~9 0°. The model eventually produced a runaway dust storm and crashed in late Year 15. Confining our study to mid-Years 5 to 15 excludes any (almost) runaway dust storm and satisfies the below explained minimum requirement of a 10-year model run. Of course, the fact that the same 2°× 2°model run produced classical GDEs decaying toward the end of the dust storm season and, likewise, (almost) runaway dust storms deserves attention. It is another manifestation of the well-known too slow decay of dust storms (e.g., Newman & Richardson, 2015).
A substantial argument for excluding (almost) runaway dust storms that occurred in the 2°× 2°model run is that they are clearly unrealistic compared against any observation-based knowledge. This is thus a reasonable data quality criterion for presenting our 2°× 2°model run in the following. Comparisons between the 2°× 2°and 5°× 5°model run are based on 10 years of data each. As described before, these are given by the mid-Years 5 to 15 of the 2°× 2°run and Years 1 to 10 of the 5°× 5°run. Note that this is practically equivalent to a more symmetric comparison (i.e., using mid-Years 5 to 15 of both model runs instead of Years 1 to 10 of the 5°× 5°run), given that there are no year-to-year variations in climate model forcing such as solar variability and orbit-spin coupling for the model setup used here. The intrinsic year-to-year variability of both model runs will be significantly different anyway with using different dust lifting parameters as explained below.
In addition to an infinite surface dust supply, the model is run with interactive dust. The latter is based on the combined use of one scheme for dust gravitational sedimentation (Forget et al., 1999) and two dust lifting schemes as described in Newman and Richardson (2015): a convective or "dust devil" scheme (Newman et al., 2002a;Rennó et al., 1998) and a near-surface wind stress lifting scheme (Basu et al., 2004;Kawamura, 1951;Newman et al., 2002a;White, 1979). The former accounts for the dust lifted by warm-core vortices that develop during the day as a result of the strong heating of the surface by the Sun and provides the annual background variability of the dust as explained below. The latter represents dust lifting by saltation and is responsible for the local/regional-to-global dust storms (Basu et al., 2004(Basu et al., , 2006 Kahre et al., 2006;Newman et al., 2002aNewman et al., , 2002b. Following Newman and Richardson (2015) and precursor studies, we assume a dust particle diameter of 4 μm and parametrize the background cycle of dust solely by the dust devil scheme. The two dust lifting mechanisms are effectively controlled by three tunable parameters, one for the dust devil scheme and two for the wind stress scheme. These are (1) the dust devil lifting rate constant α D (kg J −1 ), which is a constant of proportionality between the amount of dust lifted by dust devils per area and time, the surface sensible heat flux, and the PBL thickness; (2) the wind stress lifting threshold τ (Pa), which is the critical value of wind stress which must be exceeded for saltating particles to start ejecting surface dust; and (3) the wind stress lifting rate constant α N , which is the constant of proportionality between the vertical flux of dust and the horizontal saltation flux of sand (which depends on air density, drag velocity, and threshold wind stress).
The above interactive dust lifting parameters sensitively depend on the model resolution and have to be determined individually for any model configuration by a trial-and-error approach. In order to optimize these tunable parameters, it is common practice to match the model to observation-based atmospheric properties, such as the global T15 temperature curves (e.g., Newman & Richardson, 2015;Newman et al., 2019). The T15 temperature refers to temperatures measured in the 15 μm CO 2 absorption band, representing roughly the altitude range from 10 to 40 km, with the maximum measurement sensitivity at ca. 25 km. The weighting function used to convert temperature at different pressure levels to T15 is given in Figure 2b. The MarsWRF T15 estimates are calculated as the weighted average of the model temperatures over all except the topmost vertical model level. The latter is omitted because of possible artifacts resulting from applying Rayleigh damping at the model top, as explained before. The average weight of each model level is given by w(p)/p, with w(p) being the vertical weighting function linearly interpolated from Wilson and Richardson (2000;their Fig. 1, bottom panel, solid curve for an emission angle of 0°), given in Figure 2b, and p the vertical pressure coordinate. Average weighting by both w(p) and 1/p ensures an equal representation of all vertical layers, as our model grid has equidistant vertical pressure levels, which result in the vertical resolution increasing from the model top to bottom in terms of altitude. In order to produce global T15 temperatures, model data on T15 are averaged over 40°S to 40°N (with area weighting). As stressed in Newman and Richardson (2015), it makes more sense to match the model-predicted temperatures to those observed than to match retrievals of the dust optical depths, as a correct simulation of the optical depth does not guarantee a good model performance given that the dust radiative properties in the model may be different from those in the real atmosphere on Mars. Following Newman and Richardson (2015), the tuning of our α D , α N , and τ parameters is accomplished by comparisons between T15 temperature curves simulated by MarsWRF and those observed, with the latter shown in Figure 2a.
α D is the first interactive dust parameter to undergo tuning. A reasonable agreement between the observation-based background global T15 temperature and that simulated by MarsWRF is obtained for a dust devil lifting rate of 1.0 × 10 −9 kg J −1 , as shown in Figure 2c. Given this, α D ¼ 1.0 × 10 −9 kg J −1 is used in all subsequent 2°× 2°model runs. In the model and in addition to α D , the dust lifted by this parameterization scheme is controlled by the surface sensible heat flux (F s ; heat input to the base of the vortex) and by the parameter η, which gives the fraction of the heat input turned into work and which is a function of the PBL depth (Newman & Richardson, 2015). The amount of dust lifted by the dust devil scheme peaks around the southern spring/summer (Figure 2d) when the PBL depths and sensible heat fluxes are generally large in the southern hemisphere where solar insolation peaks.
For a 5.625°× 5°grid, Newman and Richardson (2015) obtained the most realistic infinite dust source simulation of GDEs in terms of interannual variability, onset times, and locations for dust lifting parameters τ ¼ 0.047 Pa and α N ¼ 7.5 × 10 −5 . However, for the 2°× 2°model grid considered here and consistent with Toigo et al. (2012), the wind stress lifting parameters τ and α N have to be set to different values. This meant that we had to explore a clearly different range of possible parameter values.
The tuning of τ and α N required a considerable number of sensitivity experiments. All in all, a list of τ values were tested for the 2°× 2°model grid (including 0.060, 0.065, 0.0675, and 0.070 Pa; considerable interannual dust storm variability was obtained for τ larger than 0.060 Pa) and, for each, different values of α N were considered, in the overall range 1 × 10 −5 to 2 × 10 −4 . It is worthwhile noting that when τ is increased beyond a certain value that is never exceeded at any of the model grid points, no dust storm activity can be generated by the model. When reaching a critical value of τ, there will be excess at a sufficient number of model grid points for generating a GDE in a few, but not all, years of the simulation. While a simulation time of 5 MYs (first year discarded for reasons of model spin-up) is considered in the sensitivity experiments, an in-depth assessment of the dust storm interannual variability requires running the model for at least 10 MYs, in order to consolidate statistics on GDE occurrence. The GDE classification scheme of Newman et al. (2019) is used in this work. According to this scheme, a GDE year is defined as having a T15 that exceeds the background values by more than 16 K, a possible GDE year if background values are exceeded by 10-16 K; otherwise, the MY is deemed as a nonmajor dust storm year. The respective background values are given by the dust-devil-only T15 curves in Figure 2c.

Final Selection of Parameters and Resulting Dust Storms
For the fine tuning of model parameters at 2°× 2°resolution, we did intensive testing of τ values between 0.065 and 0.066 Pa. In this range, the model has very large sensitivity to only slight parameter variations. As we found, whether the model produces a GDE in a certain year may depend on changing τ by just approximately 0.02%. This is a clear indication of having reached a critical set of dust parameters for this particular model and modeling setup. Ultimately, the final values of τ and α N were narrowed down to τ ¼ 0.0657625 Pa and α N ¼ 8.5 × 10 −5 . Note, however, that this choice of parameters has no real significance in terms of better understanding Mars, and the "best fit" parameters would certainly change if this version of MarsWRF were run with slightly different conditions (e.g., roughness map, vertical mixing scheme, and vertical grid) or if a different model were utilized altogether. On real Mars, τ and α N values are likely to vary considerably from place to place, depending on local soil cohesion, particle size distribution (of both dust and sand), etc. And, even if Mars surface dust might be taken as infinitely available, there are likely also multiple feedbacks and processes that would prevent critical behavior like the (almost) runaway dust storms, described before, from being obvious (Pankine & Ingersoll, 2002. Out of all the simulations performed, this model run best matches the occurrence of GDEs in some MYs and their absence in others. An analysis of Figure 3a reveals that there are pronounced dust-storm-related peaks consistent with observations, with a maximum temperature between 205 and 210 K in Years 9 and 14. By contrast, the maximum values of T15 are close-to-background in Years 5, 8, and 12. It is, thus, interpreted that GDEs occurred in Years 9 and 14 and were absent in Years 5, 8, and 12. Years 7 and 10 are possible GDE years, while Years 6, 11, and 13 are at the boundary between possible GDE years and GDE years. (For Years 6, 11, and 13, the maximum temperature is close to but still not clearly exceeding, 200 K.) Overall, 2-3 GDEs occurred during the 10 MY simulation presented in Figure 3a. This interannual variability is consistent with GDEs that occur, on average, once every 3-4 MYs . Hence, the model, in its present configuration, has generated a realistic rate of GDE occurrence. However, the onset times (and locations) of GDEs are far more uniform in the model than in reality.
An interesting feature is that there are no two annual T15 temperature curves in Figure 3a exactly the same. The curves have departures of several degrees K from one another beginning early in the dust storm season (e.g., the spike of the orange curve (Year 14) around L s 190-195°). It is worthwhile noting that sensitivity experiments showed that such departures also occur early in the dust storm season but tend to be enhanced, for larger values of τ (e.g., for τ ¼ 0.070 Pa). Starting at L s 250-260°, all curves in Figure 3a depart from the dust-devil-only curves given in Figure 2c by at least a few K. While some have rather gentle increases, others undergo steep increases of several K at, or shortly after, L s 250-260°. Some of these temperature increases cease at, or below, ≃200 K, which produces the possible GDE/nonmajor dust storm years according to the Newman et al. (2019) classification scheme described above, while others continue after a short pause to maximum values close to, or somewhat below, 210 K (i.e., GDEs). The obtained (possible) GDEs are thus rather uniform in time and temperature until reaching ca. 200 K but may diverge considerably after that point.

10.1029/2019JE006253
Journal of Geophysical Research: Planets Newman and Richardson (2015) can be explained by the slightly different experimental setup used here, in particular, the higher vertical resolution (45 instead of 20 levels) and a grid spacing of 5°instead of 5.625°in longitude. Still, our Figure 3b is reminiscent of Fig. 3a of Newman and Richardson (2015), their most realistic infinite surface dust simulation.
Our 5°× 5°simulation in Figure 3b resembles a binary system, that is, it produces either GDEs or is at, or below, the boundary between nonmajor dust storm and possible GDE years, which is likely because of the lower value of τ. Positive feedback effects made T15 temperature rise quickly in some MYs to maximum values of 210 K or more, while remaining clearly below 200 K in other years. Intermediate dust storm activity at the boundary between a possible GDE and GDE year was therefore only generated in the 2°× 2°model run. Although such behavior is also seen in the finite dust simulations of Newman and Richardson (2015), it may not be consistent with reality. At least in the era of continuous global spacecraft observations since the late 1990s, there has never been any doubt whether a certain event qualifies as GDE or not (Forget & Montabone, 2017, their Fig. 3).

Analysis of the Characteristics of the Martian Dust Cycle
It is important to note that any model output on surface dust lifting is the sum of all dust lifted per area during one model output time step, which, in this case, is 1 Sol (0 to 0 UTC) as described in section 2 and not instantaneous values. The long-term average of the surface dust lifting by the wind stress and dust devil parameterization schemes is shown in Figure 4. As expected, the dust lifting by the wind stress scheme for 2°× 2°( Figure 4d is as Figure 4a but for the 5°× 5°simulation. As follows, the magnitude and spatial extent of dust lifting have significant departures between the 2°× 2°and 5°× 5°model run. In the latter run, dust lifting tends to occur over a wider fraction of the planet (e.g., dust is lifted in a nearly circumglobal banded region, which extends from Tharsis to Syrtis Major and then from Isidis Planitia to the northwest of Tharsis). This is not surprising, as the wind stress lifting threshold, τ, is lower in the coarse resolution simulation. In particular, in regions where the topography is already well captured at 5°× 5°resolution, this implies additional lifted dust. Indeed, in the 5°× 5°run, N-HB and N-OM appear to be stronger dust lifting regions than in the 2°× 2°run. Regions that are activated in the 5°× 5°but not in the 2°× 2°run include (i) Terra Cimmeria and south of Elysium Planitia, (ii) the southeastern and southern parts of the Tharsis plateau, and (iii) Isidis Planitia.
S-CP and S-HB are essential dust lifting regions in the 2°× 2°model run, as discussed below, but only minor dust lifting regions in the 5°× 5°run. Such contrast is explained by differences in the modeled atmospheric circulation and topography. Figure 5 shows the near-surface wind vectors and the column-integrated dust optical depth at 0 UTC for the 2°× 2°and 5°× 5°runs around the solstices for a typical non-GDE year. As follows from comparing Figures 5b and 5d, at L s 270°, the western boundary current to the east of the Tharsis Plateau is stronger in the 2°× 2°run, where the large topographic slopes are also better represented. This is consistent with S-CP being one of the most essential dust lifting regions in this experiment, in contrast to the 5°× 5°run. As follows from comparing Figures 5a and 5c, at L s 90°, the more intense near-surface winds around S-HB in the 2°× 2°simulation are in line with the dust lifting in this region being significantly higher than in the 5°× 5°run (as shown before in Figures 4a and 4d and 4b and 4e). In interpreting Figure 5, it should be noted that the local time is different at different longitudes, which will be reflected in the prevailing local atmospheric circulations.
In order to highlight the diurnal variability of the surface wind field and how it differs between the two simulations, Figure 6 shows the near-surface wind vectors and column-integrated dust optical depth at 0, 8, and 16 UTC for Hellas Planitia, where the topographic slopes are better represented in the higher-resolution model run, and around L s 90°. Figure 7 shows the same for Isidis Planitia, where the representation of the topography is roughly the same in the 2°× 2°and 5°× 5°model runs, and around L s 270°. As expected, the steeper slopes at Hellas Planitia in the 2°× 2°simulation give rise to stronger winds than for 5°× 5°, whereas the wind strength at different times of the day at Isidis Planitia is about the same for both model runs. At Hellas Planitia and taking a reference longitude of 60°, the local time is UTC + 4 and the times of 0, 8, and 16 UTC correspond with 4, 12, and 20 local time. At the southern parts of the basin, downslope flows prevail at day and night, likely driven by the temperature gradient between the CO 2 ice-covered terrain to the south and the ice-free regions inside Hellas (Chow et al., 2018). On the northern side, on the other hand, during daytime, upslope flows are predicted by the model, with speeds larger than 10 ms −1 and much weaker downslope winds (<2.5 ms −1 ) at night. At Isidis Planitia and taking a reference longitude of 90°, the local time is UTC + 6 and 0, 8, and 16 UTC correspond with 6, 14, and 22 local time. As at northern Hellas Planitia, during daytime, upslope flows with speeds larger than 10 ms −1 prevail, possibly leading to a slight increase in the dust optical depth over the high terrain, whereas at night, clearly weaker (<5 ms −1 ) downslope winds are predicted.
In addition to the tuning of the wind stress lifting, dust devils are another critical component when setting up a model with interactive dust, as described in section 2. In Figure 4c, the spatial pattern of surface dust lifting by dust devils is provided for the 2°× 2°model run. There is a belt of dust lifting regions along the latitude circle at ≃30°S and dust lifting areas at and near prominent orography such as Olympus Mons, the Tharsis Montes, Alba Patera, and Elysium Mons. A comparison of the dust devil lifting of the high and coarse resolution simulations (Figures 4c and 4f) shows an overall similar spatial pattern but with reduced magnitudes in the 2°× 2°run, typically by 10-30% but up to 60% in some places, in particular near Olympus Mons. This is mostly because of a higher value of α D in the coarser resolution run (1.5 × 10 −9 kg J −1  Figures 8a and 8b show the change in surface dust from the first to last sol in units of kg m −2 (assuming infinite dust as described in section 2, the initial surface dust cover was set to 1,000 kg m −2 ) of the 2°× 2°and 5°× 5°simulation, respectively. As can be seen in Figure 8a, the dust lifting region S-CP is accompanied by a pronounced dust deposition region to the south, E-EM to the east and south, and N-HB and S-HB inside the Hellas basin. For Figure 8b and 5°× 5°, there is relatively strong dust gain at values larger than 0.04 kg m −2 over large parts of the planet. This is in line with the pattern of dust lifting regions being spatially more widespread and partly stronger for 5°× 5°than 2°× 2°, as described before. Although there are more regions that experience a net deposition than a decrease in the amount of surface dust, the total mass of dust at the surface and in the atmosphere is constant throughout the simulation, as explained in section 2. This is

Journal of Geophysical Research: Planets
At least in the Northern Hemisphere, the patterns of dust gain and loss of both 2°× 2°in Figure 8a and 5°× 5°in Figure 8b have a certain degree of similarity with the dust cover according to the Dust Cover Index (DCI) map of Ruff and Christensen (2002;their Fig. 14). The DCI map clearly indicates dust cover over large parts of the Northern Hemisphere from the tropics to the midlatitudes with only few exceptions such as much of Syrtis Major and Acidalia Planitia. The pronounced dust gain regions in the northern low to middle latitudes, given by red to purple color in our Figures 8a and 8b, resemble, to a certain extent, the dust covered regions in the Northern Hemisphere, given by orange to red color in Figure 14 of Ruff and Christensen (2002). A close look reveals that the before-mentioned exceptions of much of Syrtis Major and Acidalia Planitia are better matched by the dust loss regions of the 2°× 2°model run. By contrast, the 5°× 5°run has a strong and spatially widespread dust loss region to the north and east of Alba Patera. The latter is much reduced in spatial extent and partly turns into dust gain for 2°× 2°. This is consistent with the DCI map not showing clear evidence for a dust-free region to the north and east of Alba Patera. In the Southern Hemisphere, however, the agreement between the model and DCI map is more ambiguous: The model has a particularly strong dust gain region inside the Hellas Basin. That is somehow reminiscent of the run. Note that the contour intervals are not evenly spaced.

Journal of Geophysical Research: Planets
DCI map having intermediate values in and around the Hellas Basin, which are likely to indicate a partially dust covered surface (Ruff & Christensen, 2002). Likewise, the DCI map shows dust-free surfaces over large parts of the Southern Hemisphere that extend clearly further to the south than the already mentioned belt of dust lifting and the related dust loss region of the model along~30°S. This is potentially improved by the 2°× 2°model run having less dust gain at the southern midlatitudes than the 5°× 5°run. Of course, the obtained pattern of surface dust gain and loss of our model simulations is governed by the assumption that surface dust is infinitely available in source regions. If the surface dust in primary source regions was to be completely exhausted, as in the finite dust simulations of Newman and Richardson (2015), the patterns of surface dust lifting and deposition would change significantly. In fact, Newman and Richardson (2015) obtained a reasonable agreement between the modeled surface dust cover and the surface albedo map in the Southern Hemisphere in their finite dust simulation.
The patterns of dust lifting and deposition in Figure 8a are consistent with the low-level atmospheric general circulation on the planet: As evidenced by Mitchell et al. (2015, their Fig 2d), the Mars northern polar jet stream is tilted away from the pole with decreasing altitude and comes close to the surface at latitudes of around 50°N. For similar latitudes, Joshi et al. (1997) found a low-level jet that is strongest around N-OM, N-AT, and E-EM during the boreal winter/austral summer seasons (their Fig. 1, lower panel). This low-level jet, which can also be seen in Figure 5b, has a general direction from west to east and follows the wave number two topography, which is consistent with the fact that N-OM has a dust deposition region to its east to southeast and E-EM to the east and south. Joshi et al. (1997) further reported a westerly low-level jet at latitudes of ≃30°S in the same season that shows a certain degree of zonal asymmetry, and is also seen in Figure 5b. This matches the latitude of the dust lifting region N-HB and the aforementioned belt of dust lifting regions distributed along this latitude circle (given by green color in Figure 8a).
In the 2°× 2°model run, surface dust lifting by wind stress occurred at 708 out of the 16,200 points of the model grid, which is less than 5%. (This figure rises to ≃20% in the 5°× 5°run, with lifting at 535 out of the 2,592 points.) For effectively comparing different points, all of their wind-stress-lifted dust over the 10-year simulation was area weighted and the total global sum over all points was calculated. As shown in Figure 9a, the red point, alone, makes up ≃16% of the global sum. There are six yellow points that, together with the red point, account for almost 50% of the global sum of surface dust lifted by wind stress. This is followed by 189 green points that are, together with the red and yellow points, responsible for ≃99% of the global sum. By implication, the blue points are responsible for not more than ≃1% of the global sum. The red point is part of S-CP and the yellow points of S-CP, N-HB, S-HB, and E-EM. All regions from S-CP to W-AB include a number of green points. Table 2, middle column, provides the percentage contribution to the global sum of dust lifted by wind stress for each of the seven regions, as also highlighted in Figure 9a. S-CP, N-HB, S-HB, N-OM, and E-EM are the major dust lifting regions, each contributing more than 10% to the global sum. At the other end are N-AT and W-AB that do not exceed 0.3%. In the 5°× 5°s imulation, more than 75% of the global sum of dust lifted by wind stress are provided by N-HB and N-OM, which, thus, dominate other source regions (S-CP, E-EM, S-HB, N-AT, and W-AB) (see Table 2, right column). The top 9 (234) dust contributing grid points account for almost 50% (≃99%) of the global sum. The top 9 points are located on the northern slopes of the Hellas basin and north of Alba Patera (see red and yellow points in Figure 9c).
It is worthwhile noting that the dust lifting characteristics of the 5°× 5°model run are reminiscent of Newman and Richardson (2015). In their most realistic infinite dust simulation, their 6 (100) top ranked grid points contributed over 50% (99.6%) of the entire surface dust lifted by wind stress. The top 6 points are part of their most pronounced dust lifting regions, which are north and east of Alba Patera and on the northern slopes of the Hellas Basin (their Fig. 5).
The fact that there are several major dust lifting regions of similar strength at low latitudes is a key factor governing (possible) GDEs of variable strength in the 2°× 2°model run. We found that in the case of (possible) GDEs, regional dust storm activity is typically triggered both around S-CP and in the Hellas basin (N-HB). This can be seen in Figure 10, which shows a time series of all dust lifted by the wind stress scheme at S-CP, N-HB, and S-HB, which are major dust lifting regions according to Figure 9a and  6, 7, 9-11, 13, and 14), both areas exhibit significant dust lifting, with that at N-HB at times exceeding the one at S-CP. If these two regional activities connect, a GDE occurs, as highlighted in Figures 11 and 12. If not, the T15 temperature reaches values of ≃200 K but does not increase beyond that. This finding is consistent with "Hellas storms" playing a role in the development of GDEs (e.g., Chow et al., 2018;Wang & Richardson, 2015). Figures 11 and 12 show maps of the dust optical depth at six sols for a GDE year (Year 14) and a non-GDE year (Year 10), respectively. In Figure 11, the dust lifting at N-HB leads to optical depth values in excess of 1 over a large portion of the southern hemisphere by L s 281°, with the dusty conditions having extensions up to roughly 40°N. For this year, there is a clear link between the dust storm activities at N-HB and S-CP, in particular for L s 286-288°. Dust lifting at N-OM is also evident in the optical depth plots for L s 281-283°. Overall, the dust optical depth maps for Figure 9. (a) Distribution of all grid points of the 2°× 2°model run with surface dust lifting by wind stress. Overall, these are 708 out of the 16,200 grid points. The dust contribution of each such point was area weighted and the total global surface dust lifted was calculated. The red point alone accounts for ≃16% of the total global surface dust lifted. Together with the yellow (both yellow and green) points, almost 50% (≃99%) of the total global dust lifted is originated. By implication, the blue points are responsible for no more than ≃1% of the total global surface dust lifted. (b) Is the same as (a) but for the aphelion season only (i.e., Years 6 to 15, each Sols 1 to 373, corresponding to L s of 0°to 179.758°). The yellow (yellow and green) points make up less than 50% (≃99%) of the total global dust lifted. (c and d) As (a) and (b) but for the 5°× 5°simulation. Here dust lifting by the wind stress scheme occurs at 535 out of 2,592 grid points.

Journal of Geophysical Research: Planets
Year 14 resemble those in Fig. 14 of Newman and Richardson (2015). The latter is based on their simulations with finite surface dust and has similarity with the observed MY28 global storm, apart from the early onset.
In other words and despite not capturing all aspects of the real Mars dust cycle, the dust storms generated by our 2°× 2°MarsWRF configuration exhibit a certain degree of realism, in particular, regarding the associated global dust loading and interannual variability. For non-GDE years and as seen in Figure 12, the model simulated optical depth is quite different: values clearly greater than 1 occur over S-CP, N-HB, and E-EM (and nearby dust deposition regions known from Figure 8a) for short periods of time but do not obviously connect.
As shown in Figure 10, dust lifting at S-HB occurs in the aphelion season and also early in the perihelion season, before the onset of dust lifting at S-CP and N-HB. Regardless of that, the dust lifting at S-HB can reach large magnitudes, at times exceeding that at S-CP and N-HB later in the year. Moreover, the dust lifting at S-HB is generally confined to single-to-few-sol events, as opposed to that at S-CP and N-HB where it occurs in a rather continuous manner. An analysis of the model-predicted near-surface winds and dust optical depths indicates that the time difference of subsequent dust lifting events at S-HB is mostly in the rangẽ 2-7 Sols. This is the dominant time scale of midlatitude high-frequency baroclinic eddies (e.g., Haberle et al., 2018), with some additional similarity with low-frequency eddies (i.e., time difference larger than 10 Sols). In other words, dust lifting at S-HB appears to occur when the local downslope winds, visible in Figures 5a and 6, are reinforced by the winds from midlatitude synoptic-scale transients (with the latter having their origin in the southern hemisphere).
For revealing the contrasts between the perihelion and aphelion season, Figures 4a and 9a were split up into separate figures for the perihelion and aphelion season only. This was done by running the used dust lifting Figure 10. Time series of the total daily surface dust lifting by the wind stress scheme (in tons) for the 2°× 2°model run and the dust source regions of S-CP (red), N-HB (green), and S-HB (blue), respectively. Each time series was calculated as the sum over all grid points of the respective dust source region on a sol-by-sol basis. The covered time span is mid-Years 5 to 15. On the bottom horizontal axis, the time is provided as a decimal number (e.g., 5.5 is Year 6 and L s 180°). In addition, each year is classified according to the dust storm classification scheme of Newman et al. (2019): "GDE" is a "GDE year"; "P+" is a year at the boundary between a "GDE year," and "potential GDE year"; "P" is a "potential GDE year"; "NON" is a "nonmajor dust storm year". Vertical solid lines denote the start and end of a year, whereas vertical dashed lines mark L s 270°of Years 9 and 14, which are the approximate GDE onset times in these years. The solar longitude is added on the top horizontal axis. The vertical axis has a logarithmic scale (zero values are missing in the plot; consecutive nonzero values are connected by a line). statistics over the same MYs, but each only L s 0-180°(aphelion season) and L s 180-360°(perihelion season), respectively. Figures 4a and 9a and the perihelion-season-only plots are very similar; thus, the latter are not shown. This is not surprising given that the perihelion season is the time of the year when most of the dust storm activity occurs. The aphelion-season-only plots are shown in Figures 4b and 9b. They exhibit major dust lifting regions on the southern slopes of the Hellas basin (S-HB) in similarity with Figures 4a and 9a. This is because dust lifting at S-HB is connected with dust storm activity early in the MY around L s 180°( compare Figure 10, blue curve). This contributes both to the aphelion-season-only and full-year plots. In principle, related regional dust storms are similar to the "pre-spring Hellas" regional storms of Newman and Richardson (2015; see their Table 3 and corresponding text). No such regional dust storm reached the size of a possible GDE. The most promising candidate is a temperature rise of more than 5 K in Year 14, around L s 190-195°(see spike in orange curve in Figure 3a). Still, it is not unreasonable to consider such Journal of Geophysical Research: Planets events as precursors of (possible) GDEs. In line with that, the 2°× 2°model runs of this study are a step forward in the ability of reproducing the interannual variability in the onset times of observed GDEs.

Summary
In this paper, we presented a MarsWRF simulation at a horizontal resolution of 2°× 2°with 45 vertical levels using fully interactive dust lifting schemes and compared it with a 5°× 5°simulation using the same experimental setup, which is similar to the spatial resolution considered in previous studies. As explained in section 2, the high computational cost of performing this study could be stemmed by assuming an infinite surface dust reservoir. In the considered "interactive dust" configuration, there are three tunable parameters: the dust devil lifting rate constant α D (kg J −1 ), the wind stress lifting threshold τ (Pa), and the wind stress lifting rate constant α N . Out of the multiple sensitivity experiments performed, the optimal results for the 2°× 2°model run were obtained with τ ¼ 0.0657625 Pa, α N ¼ 8.5 × 10 −5 , and α D ¼ 1.0 × 10 −9 kg J −1 . In the 5°× 5°simulation, these parameters were set to 0.046 Pa, 1.5 × 10 −5 , and 1.5 × 10 −9 kg J −1 , respectively. The higher value of τ and the smaller value of α D in the 2°× 2°simulation are consistent with a better representation of the topographic slopes and hence of the associated atmospheric circulation (Toigo et al., 2012).
In the 5°× 5°simulation, surface dust lifting by wind stress takes place at roughly 20% of the model grid points. As expected, the 2°× 2°model resolution can overcome some, but not all, of the shortcomings of an infinite dust simulation at a coarser spatial resolution. In particular, dust lifting by wind stress occurred only at roughly 5% of the model grid points. In addition, the onset times of (possible) GDEs are quite uniform, at about L s 260°. Nonetheless, the 2°× 2°model simulation exhibits a number of interesting and partly innovative features. Some of the most relevant findings are highlighted below: • Major wind stress lifting regions are better defined in the 2°× 2°resolution model run. They are located mostly along the meridional corridors Acidalia-Chryse, Utopia-Isidis, and Arcadia-West of Tharsis. Dust lifting by the dust devil scheme occurs mostly along the latitude bands of ≃20-40°N/S, as well as around volcanoes such as Tharsis Montes, Olympus Mons, Alba Patera, and Elysium Mons. In the 5°× 5°simulation, wind stress lifting occurs across much broader areas, covering a wider fraction of the planet, with a nearly circumglobal strip at roughly 50°N. The contrast in the wind stress lifting regions between the two simulations is explained by the differences in the wind stress lifting model parameters, atmospheric circulation, and model topography. In the 5°× 5°resolution model run, the spatial pattern of dust-devil lifting has a high degree of similarity with 2°× 2°, but the amplitudes are increased, mostly because of the higher value of α D . • In the 5°× 5°simulation, the model produces either GDEs or nonmajor dust storm years. By contrast, in the 2°× 2°resolution model run, there is a different variability in terms of the T15 magnitude. In addition to GDEs and nonmajor dust storm years, intermediate years with a peak T15 temperature of ≃200 K are observed. The different interannual variability in the 2°× 2°model run is governed by the interaction of major dust lifting regions at low latitudes. GDEs are likely to develop if regional dust storm activity south of Chryse Planitia and northeast of Valles Marineris (S-CP) connects with that on the northern slopes of the Hellas basin (N-HB). Otherwise, the peak T15 never clearly exceeds 200 K. • In the aphelion season, L s 0-180°, there are major dust lifting regions in the 2°× 2°resolution model run, such as on the southern slopes of the Hellas Basin (S-HB). The dust lifting at S-HB can reach large magnitudes, at times exceeding that at S-CP and N-HB in the perihelion season. Interestingly, the dust lifting at S-HB shows a predominant 2-7 Sol variability, which is consistent with the time-scale of midlatitude baroclinic eddies (Haberle et al., 2018). This suggests that the interaction between the high-frequency transients and the local downslope flow drives the dust lifting in the region. MarsWRF produces significant dust storm activity in some MYs near L s 180°at the 2°× 2°resolution. This is further in line with winds in the vicinity of the southern polar cap edge and dust lifting on the southern slopes of the Hellas basin (Chow et al., 2018). By contrast, S-HB is not a major dust lifting region in the 5°× 5°resolution simulation.
It is important to note that the finite surface dust 5°× 5.625°simulation of Newman and Richardson (2015) was able to reproduce far more widespread dust lifting and distribution of primary dust source regions than the infinite surface dust 2°× 2°model run presented here, as well as far more variability in dust storm onset time and location. In addition, the 5°× 5.625°simulations were able to generate very realistic Chryse "flushing" dust storms, Noachis storms, and also southwest Hellas storms around L s~1 80°. It is interesting that several of the improvements associated with finite dust are similar to those obtained here with increased spatial resolution. This suggests that further realism, in particular with respect to the global dust lifting and interannual variability of the dust cycle, would be obtained if repeating the 2°× 2°model run with finite dust. Enhancing the horizontal model resolution beyond 2°× 2°thus appears an even more interesting outlook for the future.