Characterizing Dust‐Radiation Feedback and Refining the Horizontal Resolution of the MarsWRF Model Down to 0.5 Degree

In this study, three simulations by the Mars Weather Research and Forecasting Model are compared: two 10 Martian year (MY) 2° × 2° simulations with (i) fully radiatively active dust and (ii) a prescribed dust scenario, and a (iii) 1 MY 0.5° × 0.5° simulation with prescribed dust as in (ii). From comparing (i) and (ii), we found that the impact of dust‐radiation feedback is individually different for any region. The most striking evidence are major dust lifting activities to the south of Chryse Planitia (S‐CP) seen in (i) but not in (ii). By contrast, dust lifting and deposition on the southern slopes and inside the Hellas Basin are similar in both simulations. The latter, in turn, points toward a similar near‐surface atmospheric circulation. In (iii), the total global amount of wind stress lifted dust is by a factor of ∼8 higher than in (ii), with S‐CP being a major lifting region as in (i). Nonetheless, the surface dust lifting by wind stress in (iii) may be also reduced regionally, as seen at the peak of Elysium Mons because of its unique topography. The zonal mean circulation in (i) is generally of a comparable strength to that in (ii), with exceptions in global dust storm years, when it is clearly stronger in (i), in line with a dustier atmosphere. The differences in the zonal mean circulation between (ii) and (iii) are mostly at lower altitudes and may arise because of differences in the representation of the topography.

develop as a result of the steep surface-air temperature gradient in response to the strong heating of the surface by the Sun (Newman & Richardson, 2015). Parameterization schemes are an attempt to represent the subgrid-scale processes that cannot be resolved explicitly by the model and may be thought of as an approach to estimate the cell-averaged behavior. Furthermore, it is common practice to use globally uniform values for certain parameters, whereas their equivalents on real Mars likely exhibit pronounced spatial variability. One day, it may even be possible to quantify relevant surface properties and processes comprehensively through networks of orbital and/or in situ observations. This would be invaluable given that small changes in the surface parameters that describe the dust lifting process, such as the wind stress lifting threshold and efficiency, and particle size distribution, etc., may result in drastic changes in the resulting model simulations.
To date, the values of model dust parameters are estimated to the best of our knowledge empirically and/or following practical considerations from the modeling side. It is possible to calibrate or "tune" such parameters by comparing model simulations with orbital observations of certain measurable quantities, such as satellite data of the atmospheric temperature and others. This approach is followed by a number of studies such as Basu et al. (2004), Newman and Richardson (2015), Newman et al. (2019), and Gebhardt et al. (2020a), just to name some of them.
Since there is pronounced feedback between atmospheric dust and radiation, it is known that vertical mixing and the spatial and temporal patterns of surface dust distribution depend strongly on whether a model simulation includes such feedback. In addition, previous studies demonstrated that surface dust lifting may vary strongly with the spatial resolution of a model simulation (Toigo et al., 2012). To illustrate and quantify the impact on the dust cycle arising from dust-radiation feedback and varying spatial resolution, we perform here a comparison study of three model simulations: two with 2° × 2° and one with 0.5° × 0.5° horizontal resolution.

The Dust Cycle and Dust-Radiation Feedback
The Mars atmosphere and climate is subject to a pronounced annually repeating seasonal variability. The latter manifests itself as the CO 2 , water, and dust cycles (Zurek et al., 1992). The importance of mineral dust for the Mars atmosphere and climate is explained as follows. Dust heats the atmosphere by absorbing incoming solar radiation and outgoing thermal infrared radiation (e.g., Madeleine et al., 2011;Newman & Richardson, 2015, and references therein) and therefore has essential effects on the atmospheric temperature and dynamics.
The time around perihelion (solar longitude, L S , from 180° to 360°), which coincides with Mars southern hemisphere spring/summer, is the Mars dust storm season. It is characterized by a relatively warm and dusty atmosphere from the southern hemisphere to the northern midlatitudes. The main driver is the Martian hemispheric dichotomy in topography (M. Richardson & Wilson, 2002). The higher terrain of the southern hemisphere enhances the global overturning circulation, which is accompanied by stronger surface dust lifting and transport. While local/regional dust storms may occur at any time of the year, global dust storm events (GDEs) have only been observed during the dust storm season (Shirley, 2015;Wang & Richardson, 2015).
In general, Martian dust is lifted from the surface into the atmosphere, transported by the atmospheric circulation, and deposited back onto the surface. While the latter is controlled by gravitational sedimentation (Forget et al., 1999), major dust lifting mechanisms are given by dust devils (M. Balme & Greeley, 2006;Neakrase et al., 2016) and near-surface winds (Kawamura, 1951;White, 1979). As described below, the latter is closely linked to the occurrence of dust storms from local to global size. Modeling and observation-based work point toward dust devils and local/regional dust storms (the latter are initiated by surface wind stress) lifting a similar amount of surface dust throughout the Martian year (MY) Kahre et al., 2006;Whelley & Greeley, 2008).
Dust devils are warm-core vortices that develop during daytime as a result of the solar heating of the surface (Kahanpää et al., 2016). On Mars, they are known to have a maximum size of several hundred meters in diameter and few kilometers in height (M. Balme & Greeley, 2006;Cantor et al., 2006, and references therein). From Earth, dust devils are known to have durations of up to several minutes (M. Balme & Greeley, 2006).
The fact that dust devils do not become larger in area and/or duration is explained by a negative feedback mechanism (Newman & Richardson, 2015;Newman et al., 2002a). The absorption of radiation by the lifted dust results in a warming of the atmosphere and subsequent cooling of the surface. As a result, when the dust loading increases, the vertical surface-atmosphere temperature gradient and convection of warm air are thus reduced, which inhibits further dust devil growth. It was found that dust devil activity can drop to near zero under very high atmospheric dust loads observed during dust storms (Cantor et al., 2006;Greeley et al., 2010;Guzewich et al., 2019;Lemmon et al., 2015). Likewise, dust devils are an important contributor to the annually repeating background cycle of Mars dust, particularly outside the dust storm season (e.g., Basu et al., 2004). Dust devils are found not to be correlated or even anti-correlated with the initiation and growth of dust storms (M. R. Balme et al., 2003;Cantor et al., 2006;Greeley et al., 2006Greeley et al., , 2010Lemmon et al., 2015;Newman et al., 2002b;Whelley & Greeley, 2008). Both observations and modeling corroborate that the main and direct driver of dust storms is dust lifting by near-surface wind stress, and accompanying positive feedback mechanisms (Newman & Richardson, 2015;Read & Lewis, 2004, their Section 7.5.3.3).
It is important to note that our following model runs at a horizontal grid of 2° × 2°, which is synoptic resolution, are not one-to-one comparable to the idealized model experiments of Rafkin (2009), with their mesoscale resolution. The latter, nonetheless, offers a compact overview of possible mechanisms of radiative-dynamic feedback at the local/regional scale and related dust storm growth. Interestingly, the introductory part of Rafkin (2009) outlines a possible chain of processes, into which such local feedback may be broken down: (i) uplift of surface dust into the atmosphere by the wind; (ii) the lifted dust leads to enhanced radiative heating (less radiative cooling) during the day (night); this results in a relatively warm region in which the surface pressure is decreased hydrostatically; (iii) an increased near-surface pressure gradient force and resulting stronger winds that give rise to enhanced dust lifting. Also, the summary part of Rafkin (2009) illustrates the fact that pronounced topographic gradients are connected with enhanced positive local feedback: the surface temperature is mainly governed by thermal inertia, albedo, and solar insolation, whereas the atmospheric temperature is rather a function of altitude. As a consequence, geolocations that are elevated with respect to their surroundings thus have a tendency toward higher lapse rates and stronger updrafts/ adiabatic cooling of air/conversion of thermal into mechanical energy. The dust's radiative properties may also allow for its own vertical transport, as the heating which results from the absorption of solar radiation increases its buoyancy. Such "rocket dust storms," as denoted by Spiga et al. (2013), can explain the existence of detached dust layers up in the Martian troposphere.
In addition to local scales, the wind stress lifting has a strong positive feedback response on a global scale, as the overturning circulation intensifies for higher dust loadings (Newman & Richardson, 2015) via direct heating and often also amplification of tides (Wilson, 1997), which results in generally stronger winds over a very large region. Among others, this is based on Newman et al. (2002a) who found that the modeling of the atmospheric dust transport depends sensitively on the vertical mixing of dust in the planetary boundary layer (PBL). Positive feedback mechanisms enhance the local vertical velocity and hence are crucial to dust entering the planetary-scale circulation. As a consequence, surface dust may be lifted and transported from one hemisphere to the other by the cross-equatorial dust storms within days (Haberle et al., 2017, their Section 7.3.3;Newman et al., 2002b). Toigo et al. (2018) examined the thermal and dynamical response of the global Martian atmosphere to local and regional dust storm radiative heating. Their approach was based on idealized experiments, which approximated dust storms by increasing the atmospheric dust loading for certain regions and duration to emulate the effect of different local and regional storms, in a 2° × 2° simulation with prescribed dust loading. This included perturbations of local/regional size, 1.5 and 10.5 sols duration, and dust optical depth values between 1 and 250, at different locations. These experiments were designed to identify the conditions required to trigger self-induced dynamical growth of a storm or precondition the atmosphere for a follow-on dust storm. The authors concluded that local-scale storms, which cover an area of the order of 10 5 km 2 , and regardless of 1 or 10 sols, do not change the thermal or dynamical state of the atmosphere beyond the forcing region by amounts larger than the background environmental variance. The atmospheric state is only significantly perturbed for regional storm length scales, which cover an area of the order of 10 6 km 2 , and that are long lasting (10 sols). If the heating is placed in a region that comprises the upwelling branch of the Hadley Cell (e.g., Hellas Planitia in the austral summer), a stronger overturning circulation and subsequent increase in optical depth is predicted to occur. However, when placed in an area that includes the descending branch of the Hadley Cell, it will act to weaken it. In other words, the feedback of a regional dust storm on the overturning circulation will depend on its location, in addition to the spatial and temporal scales.
We use in the following an approach that is complementary to the artificial construction of dust storms by idealized model experiments. In brief, we configure the Mars Weather Research and Forecasting (Mars-WRF) Model to self-consistently produce dust storms, which originate from different major source regions. This allows us to capture synoptic-to-planetary-scale feedback in its entirety. As an integral part, this includes the impact of feedback on the initiation of new dust storms. Moreover, we encompass complex interactions between different dust storms and the large-scale atmospheric circulation as well as the self-induced growth of dust storms. At least when it comes to GDEs, several active dust lifting centers are involved. All this is interesting because if more dust storms are initiated at a certain source region with radiatively active dust, a higher number of them may, in turn, pass through the cascade described in Toigo et al. (2018).
By implication, our approach means that feedback on local/regional scale is convolved with that on global scale. Such mixed feedback is the basis of this article. It is analyzed quantitatively by comparisons of surface dust loss between model runs with and without fully radiatively active dust. Differences between the former and latter are an indicator for changes regarding dust storm initiation/sustainment/amplification, taken as a whole. It is part of the nature of this approach that local/regional-and global-scale feedback as well as their impact on dust storm initiation, sustainment, and amplification typically cannot be totally disentangled from one another. As an integral part, this includes the dynamic response to a dustier atmosphere in general (or, at least, in some regions), associated with regional to global-scale dust storms. The global atmospheric circulation is enhanced in the presence of a dustier environment, which is strong positive feedback (Leovy et al., 1973;Newman et al., 2002a). Also, there may be negative feedback leading to weaker local flow regimes. This may arise from a reduced local thermal contrast and is a possible mechanism for the self-shut down of dust storms as noted in Newman et al. (2002a).

Concept and Structure of This Article
This article provides quantitative study on dust-radiation feedback mechanisms based on the simulation of the Martian dust cycle by the MarsWRF model (M. I. Richardson et al., 2007;Toigo et al., 2012). This is accomplished by repeating a model run with fully interactive dust (i.e., fully radiatively active dust) and a horizontal resolution of 2° × 2° presented in Gebhardt et al. (2020a), but with a prescribed dust scenario. Moreover, we rerun this model configuration with prescribed dust at the even higher horizontal resolution of 0.5° × 0.5°. The model starting from the scratch requires two full MYs of simulation (the first MY is for model spin-up, and the second MY for scientific analysis). This is followed by comparing patterns of surface dust lifting by wind stress between all the referred model simulations. Surface dust deposition is further compared between the 2° × 2° simulations with interactive and prescribed dust. The implications for the dynamics of dust storms are analyzed and discussed.
This article is organized as follows. Section 2 is an overview of relevant methods of dust cycle simulation, such as "prescribed dust" and "interactive dust," and the effect of horizontal model resolution. In Section 3, the model runs with interactive and prescribed dust are compared, both having a horizontal resolution of 2° × 2°. This is followed by presenting and discussing the model run with 0.5° × 0.5° resolution in Section 4. The findings of this study are summarized in Section 5. A particular focus is the implications for the dynamics of dust storms.

The Dust Cycle in Mars General Circulation Models
Mars General Circulation Models (MGCMs) offer valuable insight into the Martian climate by filling gaps in time and space of observational data sets and characterizing atmospheric variables that are challenging to quantify on an observational basis (e.g., winds; Newman et al., 2017). There are various techniques of constraining MGCMs by observational data, such as data assimilation methods (e.g., Greybush et al., 2012;C. Lee et al., 2011;Lewis & Barker, 2005) and prescribing certain data fields, such as the three-dimensional dust distribution (Guzewich et al., 2013) or total column amount of dust. An overview of the significance of the Mars dust climatology for constraining MGCMs is given in Montabone et al. (2015).
Alternatively, MGCMs may be run in interactive dust mode. In this unconstrained method, the model freely lifts, transports, and deposits surface dust in a self-consistent manner. The advantages with respect to prescribed dust scenarios are given by (i) the explicit simulation of surface dust lifting/deposition, (ii) the ability to predict the actual dust amount at a specific location and time of the year on the planet, and (iii) capturing dust-radiation feedback. Detailed information on running the MarsWRF MGCM with interactive dust is given in Newman and Richardson (2015), C. Lee et al. (2018), and Gebhardt et al. (2020a). In a nutshell, there are two dust lifting schemes. The first is the dust devil lifting, which accounts for the dust lifted by daytime convective vortices that develop in response to the strong heating of the surface by the Sun. The amount of dust lifted by this scheme is proportional to the surface sensible heat flux (i.e., the heat input to the base of the vortex) and the thermodynamic efficiency (i.e., the fraction of the heat input that is converted into work driving dust devils). The constant of proportionality is α D . The second scheme is the wind stress lifting, which has an assumed surface wind stress threshold, τ. If the latter is exceeded, sand particles begin saltating on the surface. The horizontal sand flux is a function of the surface wind stress, surface wind stress threshold, and the atmospheric density and proportional to the vertical dust flux. The constant of proportionality, α N , is together with τ and α D a tunable parameter. Interactive dust is therefore controlled by three user-specified dust lifting parameters.
The basis of Gebhardt et al. (2020a) were MarsWRF simulations with fully interactive dust at a horizontal resolution of 2° × 2° and 5° × 5° over a period of 10 MYs each. These used a vertical grid with 45 levels, having the greatest vertical resolution in the PBL, and a model top at 120 km (in units of midlayer eta coordinates, the uppermost model layer is ca. 0.01; compare with Table 1 of Gebhardt et al. [2020a]). The thermal inertia, albedo, and topography map of the model is based on observational data from the Mars Global Surveyor (MGS), such as the Mars Orbiter Laser Altimeter 1/64° topography (Smith et al., 2001) and Thermal Emission Spectrometer (TES) albedo and thermal inertia maps (Putzig & Mellon, 2007). A simple CO 2 microphysics scheme, the Medium Range Forecast model as PBL scheme (Hong & Pan, 1996), and a k-distribution radiative transfer model (Mischna et al., 2012) were employed. A zonal filter at latitudes poleward of 70° (M. I. Richardson et al., 2007) and Rayleigh damping at the topmost vertical model level (Skamarock et al., 2008) were applied to ensure the stability of the model run. An infinite reservoir of surface dust was considered in those runs. After intensive testing, for the 2° × 2° simulations, the following combination of parameters was eventually selected for in-depth scientific analysis: α D = 1.0 × 10 −9 kg J −1 , τ = 0.0657625 Pa, and α N = 8.5 × 10 −5 . With this set of values, MarsWRF was found in Gebhardt et al. (2020a) to predict a midlevel atmospheric temperature (denoted as T15 temperature) that best matches observation-based knowledge.
Building-upon this work, this study will address the Martian dust cycle at 2° × 2° and at an even higher spatial resolution. First, the aforementioned 2° × 2° model run with interactive dust of Gebhardt et al. (2020a) is repeated with prescribed dust and otherwise identical model configuration. This means that surface dust lifting by dust devils and wind stress still takes place but is ignored by the model's radiative transfer scheme (i.e., the associated dust-radiation feedback is not accounted for). Instead, a prescribed dust scenario is inserted into the radiative transfer scheme, namely the widely used Mars Climate Database (MCD) MGS dust scenario (Montmessin et al., 2004;Toigo et al., 2012). This dust scenario was originally developed to fit thermal profiles observed via MGS Radio Occultation measurements during MGS years without a major dust storm. As a test, we calculated the atmospheric T15 temperature as in Gebhardt et al. (2020a). The T15 produced by the MCD MGS dust scenario is consistent with the observed background atmospheric temperature. This prescribed dust simulation is then repeated at a horizontal resolution of 0.5° × 0.5° with an otherwise identical experimental setup, except for using a model integration time step of 30 s, which is half that used in the 2° × 2° simulation.
Taken together, three model runs are the basis of this study: (1) A simulation with fully interactive dust at a horizontal resolution of 2° × 2°, set up in the same way as in Gebhardt et al. (2020a), and referred to in the following as simulation IC (interactive coarse resolution), (2) a 2° × 2° model run with prescribed dust with all other model parameters identical to those in (1) and referred to in the following as simulation PC (prescribed coarse resolution), and (3) a 0.5° × 0.5° model run referred to in the following as simulation PH (prescribed high resolution). All model parameters other than the horizontal resolution and model integration time step are identical to simulation PC.
Differences between simulations IC and PC are not only an indicator of local/regional dust-radiation feedback mechanisms. As an integral part, they also include the effects of a dustier atmosphere in general (i.e., in a global sense), as noted at the end of Section 1.2. On the other hand, the contrasts between simulations PC and PH give clues regarding the impact of the model resolution. Patterns of surface dust lifting and deposition are inferred over a period of 10 MYs for simulations IC and PC and 1 MY for PH. For simulation IC, this period is identical to that selected for in-depth scientific analysis in Gebhardt et al. (2020a), where also more information on the symmetry of comparing model runs can be found. For simulations PC and PH, the first MY of the model run is discarded for reasons of model spin-up and the following 10 MYs and 1 MY, respectively, is analyzed in this study. Using just 1 MY for simulation PH is unavoidable by the factor four increase in horizontal resolution, and the associated far higher computational cost. Just to give readers an idea, 1 MY comprises 669 (sols) × 24 (h per sol) × 60 (min) × 60 (s)/30 (model integration time step used in the PH run, in s) = 1,926,720 integration steps. In addition, the 0.5° × 0.5° model grid has 720 × 360 = 259,200 horizontal grid points and 45 vertical levels, which gives a total number of 720 × 360 × 45 = 11,664,000 model grid boxes. For stemming this computational load, we employed four computing nodes at 36 processors each (i.e., 144 processors in total). Note that the use of even more nodes/ processors does not necessarily increase the speed of the model run (Laux et al., 2013). This resulted in an elapsed computation time of about 38 days per MY! The effect of model resolution on the surface wind stress and related dust lifting was systematically investigated in Toigo et al. (2012). Using the MCD MGS prescribed dust scenario, also considered here, maps of the surface wind stress were compared between resolutions of 5° × 5°, 2° × 2°, and 0.5° × 0.5°. It was found that the surface wind stress is governed by a linear superposition of the large-scale circulation and topography-driven small-scale circulations. The authors concluded that surface wind stress maps for coarse horizontal resolution are practically spatial averages of those at higher resolution. When increasing the horizontal resolution, regions of peak surface wind stress that are highly correlated with topographical gradients (e.g., canyons, basins, volcanoes, and craters) come to the fore. Compliant with that, the total global mass of wind stress lifted dust does not change significantly between different resolutions, if a close-to-zero wind stress lifting threshold is considered. For larger thresholds, however, surface dust lifting is clearly more effective for finer than for coarser resolutions. For example, for the largest threshold considered in Toigo et al. (2012) (i.e., 0.049 N m −2 ), the total mass of lifted surface dust in the 2° × 2° simulation exceeds that of the 5° × 5° simulation by a factor of more than 2 (or 100%). This is explained as follows. In many places, a large surface wind stress threshold is hardly exceeded in the coarse resolution run, while in the high-resolution simulation it is, with the dust lifting proportional to a function that includes the surface drag velocity cubed. When considering fully radiatively active dust, the above contrasts are expected to be enhanced further by positive feedback mechanisms. The simulations performed for this study and described in the following will allow for a quantitative analysis of this effect.

Interactive Versus Prescribed Dust at 2° × 2° Horizontal Resolution
The model topography of simulation IC and PC is shown in Figure 1a, which will be revisited in Section 4. In infinite dust simulations, surface dust lifting by the wind stress is typically confined to a number of spatially isolated but pronounced dust lifting regions. As explained in Newman and Richardson (2015), this is because the surface dust reservoir is never exhausted, and hence lifting will occur continuously in areas where the wind stress typically exceeds the threshold for dust lifting. Here, these regions are assigned acronyms based on their geographic position, following Gebhardt et al. (2020a): S-CP is south of Chryse Planitia; N-HB/S-HB are on the northern/southern slopes of the Hellas Basin, respectively; E-EM is to the east of the peak region of Elysium Mons; N-OM extends from north of Olympus Mons to Acidalia Planitia; N-AT is north of Arabia Terra from the Acidalia to Utopia Planitiae; W-AB is on the western slopes of the Argyre Basin. Figure 2 shows long-term averages of the daily sum of the surface dust lifting by wind stress and dust devils for simulations IC and PC, which are overall similar to one another. For simulation PC, the wind stress lifting is much weaker at S-CP and similar but somewhat weaker than simulation IC at the other dust lifting regions. On the contrary, the dust devil lifting has the overall tendency to be similar to somewhat stronger for simulation PC than for simulation IC. The latter is expected, as simulation IC produces regional to global dust storms as described in Gebhardt et al. (2020a). As generally known, dust storms warm the atmosphere and cool the surface. The resulting reduced surface-air temperature gradient acts as a negative feedback on GEBHARDT ET AL.
10.1029/2020JE006672 7 of 23 the dust devil growth (Kahanpää et al., 2016, and references therein). In that respect, simulation IC differs from PC, as the latter uses dust forcing representative of a year without any major dust storms. To summarize, both the smaller wind stress lifting and the similar or larger dust devil lifting are consistent with the absence of respectively positive and negative dust-radiation feedback mechanisms in simulation PC. This is a sound basis for comparing the wind stress lifting quantitatively between simulations IC and PC, both globally and for individual dust lifting regions.
GEBHARDT ET AL.
10.1029/2020JE006672 8 of 23  The red point alone makes up less than 25% of the total lifted dust. (c) As (a) but for simulation PH. While in (a and b) there is only one red point, in (c) there is a number of red points, which together make up just less than 25% of the total lifted dust. IC, interactive coarse resolution; PC, prescribed coarse resolution; PH, prescribed high resolution.
The percentage contribution to the total global mass of wind stress lifted dust is shown in Figure 3a for simulation IC and in Figure 3b for simulation PC. Again, the dust lifting region S-CP exhibits the most obvious change, from a percentage contribution of 33.2% for simulation IC to close to 0% for simulation PC. At the midlatitudes, the percentage contributions of S-HB and N-OM respectively increase from 20.5% and 12.3% to 48.5% and 28.3%, which for both regions corresponds to more than doubling when prescribed dust is used. By contrast, the percentage contributions of N-HB and E-EM decrease from respectively 21.5% and 11.9% to 17.1% and 5.3%, hence decrease by respectively a factor of 1.3 and 2.2 when prescribed dust is used. All in all, the regions S-HB, N-OM, N-HB, and E-EM show significant dust lifting also in simulation PC, that is, without full dust-radiation feedback. This is evident from the fact that each of them includes a yellow or red point in Figure 3b, and their percentage values are between 5.3% and 48.5%. Compared to other dust lifting regions, N-AT and W-AB are of minor importance for both simulations IC and PC. In Figures 3a-3c, the dust lifting regions defined previously are represented by rectangles. Note that the latter are similar, but not exactly identical, to Figure  9a of Gebhardt et al. (2020a). The difference is that the rectangles used here for representing major dust lifting regions are somewhat wider in order to comprise all relevant model grid points, equally for simulations IC, PC, and PH. Among others, the rectangle for S-CP comprises large parts of the Valles Marineris and that for E-EM also includes Hecates Tholus (the latter is another volcano of the Elysium province; to the northnorth-east of Elysium Mons and of smaller size).
Let p ic (p pc ) denote the percentage dust contribution in Figure 3a and simulation IC (Figure 3b and PC) for a certain region. The model output of dust lifting rates is given in units of kg m −2 and is a cumulative value over the time between subsequent model output time steps, which here is 1 sol. In the following, t ic (t pc ) is the associated total global mass of surface dust lifted by the wind stress from the first-to-last-sol of the model run IC (PC) presented here, in units of kg. In addition, r ic (r pc ) is accordingly the total mass of surface dust lifted by the wind stress at a certain dust lifting region during simulation IC (PC). We found that the global wind stress lifting of simulation IC exceeds that of PC by a factor of 3.2: For a given dust lifting region, the wind stress lifting for simulation IC (r ic ) can be related to that of simulation PC (r pc ) by As follows from inserting the earlier described percentage values into Equation 2, r ic exceeds r pc by a factor of 1.3/1.4/3.9/7.0/>100 for S-HB/N-OM/N-HB/E-EM/S-CP, respectively (see also Table 1). In line with Section 1.2, the increase factors between 1.3 and >100 represent how the dust storm initiation/sustainment/ amplification, as a whole, change from simulation PC to IC.
The obtained increase factors allow for constraining the impact of feedback on the strength, duration, and number of dust storms in certain regions. Both in simulations IC and PC, the dust source regions S-HB, N-OM, and E-EM exhibit apparent similarity and thus merit further investigation. Their dust storm activities consist of frequent dust lifting episodes with intermittent periods of typically a few to several sols, which is overall an annually recurring phenomenon. This is shown in Figure 10 of Gebhardt et al. (2020a) for simulation IC and region S-HB (but not by additional figures here). In general, the drivers for that are regularly repeating strong winds, associated with upslope and downslope flows and/or baroclinic waves, which may reinforce each other.
GEBHARDT ET AL.

S-CP >100
Note. For each, the amount of lifted dust for simulation IC exceeds that of simulation PC. The excess factor f is given in this table.

Table 1 Simulations IC and PC were compared in terms of the Total Mass of Surface Dust Lifted by the Wind Stress, both Globally and for Major Dust Lifting Regions such as S-HB, N-OM, N-HB, E-EM, and S-CP
Interestingly, the region E-EM has just one dust contributing model grid point in simulation PC (see Figure 3b). This increases to three dust contributors in simulation IC, which are the same as in PC and two neighboring grid points (see Figure 3a). There can hence be just one ongoing dust storm at a certain time at E-EM, having a dust lifting center of one to three model grid points in simulation IC, and just one in simulation PC. That is a good basis for a one-to-one comparison of dust storms at E-EM between simulation IC and PC. This and the above factor of 7.0 increase are clear evidence that dust storms at E-EM are stronger, longer, and of larger number in simulation IC than in PC. Moreover, this applies for any year of simulation. As pointed out before, the dust lifting at E-EM is an annually recurring phenomenon. By implication, each of the before, that is, the increase in the average annual dust storm strength, duration, and number cannot be larger than the before factor of 7.0. Of course, this is only a weak upper limit, or so-called constraint. Following this, it is straightforward to assume that dust storms also become longer, stronger, and more frequent at the regions S-HB and N-OM in any year of simulation when full dust-radiation feedback is included. Accordingly, the increase in the average annual dust storm strength, duration, and number from simulation PC to IC (i.e., by adding full dust-radiation feedback) has an upper limit/constraint given by the above factor 1.3 for S-HB and 1.4 for N-OM.
The regions S-CP and N-HB cannot be constrained in the same way as S-HB, N-OM, and E-EM, as they exhibit a nonuniform annual dust storm occurrence. In simulation PC, there is negligible dust storm activity at S-CP, as noted before. In simulation IC, N-HB does not show any significant dust lifting in nonmajor dust storm years (Gebhardt et al., 2020a; their Figure 10). This points toward the dust storm activity at N-HB being governed by sporadically strong winds, which occur in some MYs but not in others. These may, in turn, be triggered by remote dynamical feedback, or so-called teleconnections (Toigo et al., 2018), originating from S-CP. If dust storm activity occurs in simulation IC both at S-CP and N-HB in the same MY, S-CP always precedes N-HB (Gebhardt et al., 2020a; their Figure 10).
It is interesting to note that it is questionable whether simulation PC could self-consistently initiate dust storms at S-CP, reminiscent of the dust clouds imposed to a similar model run in the idealized experiments of Toigo et al. (2018). That is because S-CP is only a dust lifting region of minor importance in PC, as stated before. Also, the fact that dust storm activity showed a response, which is highly variable from region to region, when adding full dust-radiation feedback in simulation IC is in itself worthwhile being reported. This is in contrast to the idealized experiments of Toigo et al. (2018), who did not take into account the impact of feedback on the self-consistent initiation of dust storms by the model, giving similar results for quite different locations on Mars. Because of their degree of similarity, they reported results for only two out of the six Mars locations considered.
The fact that major dust lifting in simulation IC is opposed by close-to-zero dust lifting in simulation PC, makes the region S-CP optimal for studying the nature of feedback processes. A preliminary analysis revealed that the near-surface winds to the immediate south-east of the top dust contributing model grid points of S-CP (given by a red and yellow point in Figure 3a at 7°S and 35°W and 33°W, respectively) are clearly stronger in simulation IC than PC. This can be seen in Figure 4, which shows the near-surface horizontal winds and dust opacity of the lowest model layer in simulation IC at Year 1, sol 517, and at the same for simulation PC without dust storm activity. As follows from Figures 4a and 4b, the near-surface winds to the north of S-CP, in Acidalia and Chryse Planitia, are clearly stronger in simulation PC than IC. As highlighted in a zoom view of S-CP in Figures 4c and 4d, simulations IC and PC have strikingly identical near-surface winds in the northern to north-western corner of S-CP. Figure 4c indicates a dust cloud for IC, which starts with a pronounced maximum centered around about 10°S and 35°W, including the aforementioned top dust contributing grid points of S-CP and thus likely to be an active dust lifting center. From there, it extends to the south-east and has another local maximum at roughly 25°S and 20°W. This dust cloud is governed by enhanced near-surface wind speed in its interior, compared to the nondust storm situation in Figure 4d, and marked wind speed gradients at its edge. This indicates that there is positive local-scale feedback between the dust storm activity and the near-surface wind speed toward the south-east of S-CP. In addition, it demonstrates that the advection of dust is crucial regarding the modeling of dust-radiation feedback.
In case of a GDE year, simulation IC exhibits a stronger global circulation compared to simulation PC. In Figures 5a and 5b Hadley circulation has its ascending branch at around 30°S and descending branch at roughly 30°N, with the expected cross-equatorial return flow close to the surface in between the two. By and large, the Hadley circulation is stronger in simulation IC, as seen in both the horizontal and vertical wind components. The wind speed in the polar vortex increases with altitude. One of the most notable differences between simulations IC and PC is the near-surface zonal jet at about 20°S-30°S (Joshi et al., 1997). It has maximum zonal wind speeds of more than 20 m s −1 in simulation IC, being clearly stronger when compared to PC. This is consistent with the dust source region N-HB, located at roughly the same latitude as the referred jet, undergoing strongest surface dust lifting in GDE years, as opposed to close-to-zero dust lifting in other years (Gebhardt et al., 2020a;their Figure 10).
In Figures 6a and 6b, the global patterns of surface dust gain/loss for simulations IC and PC are provided, respectively. There is a clear correlation between the dust loss for simulations IC and PC, apart from the region S-CP as already noted previously. The hot spots of wind-stress-related dust loss at N-OM, E-EM, N-HB, and S-HB are somewhat weaker in simulation PC, in line with the smaller amounts of dust lifting shown previously in Figure 2. Also, the surface dust gain patterns have notable similarity between simulations IC and PC. A particularly interesting feature is that, both in simulations IC and PC, the Hellas Basin is a strong dust gain region. Regarding this, an interesting case study is provided in Ogohara and Satomura (2008). They simulated the dust transport in the Hellas Basin based on idealized experiments. These were based on dust injection, given by an analytical function, at the south-western rim of the Hellas Basin at L S 180°, which is similar to dust lifting occurring in our dust lifting region S-HB around L S 180° in simulation IC, GEBHARDT ET AL.

10.1029/2020JE006672
12 of 23  as described in Gebhardt et al. (2020a). The resulting dust clouds of Ogohara and Satomura (2008) initially moved northward along the western rim of the Hellas Basin and then eastward. This was followed by a clockwise rotation after reaching the northwest of the Basin. However, the resulting dust cloud was still confined to the Hellas Basin and did not escape it. This illustrates how dust lifting at S-HB may be linked to dust gain that is relatively uniform all over the interior of the Hellas Basin.
In both simulations IC and PC, there is moderate but spatially widespread dust loss at the northern to southern midlatitudes. This includes the latitude bands around 30°N and 30°S, with a certain degree of zonal asymmetry, and large parts of the western boundary current and dust storm propagation corridors GEBHARDT ET AL.

10.1029/2020JE006672
14 of 23 Acidalia/Chryse, Utopia/Isidis, and Arcadia Planitia (given by green color in Figures 6a and 6b). By contrast, this is opposed by pronounced dust gain over Arabia Terra, both for simulations IC and PC. As follows from Newman et al. (2005), regions that are in between western boundary currents are favored dust deposition areas because of their relatively low wind speeds. The fact that other similar regions, such as Elysium Planitia and to the east and west of the Tharsis volcanoes, are stronger dust gain regions for simulation IC (full dust-radiation feedback) than for simulation PC (no full dust-radiation feedback) may be explained as follows. It is known that global dust-radiation feedback mechanisms are of key importance for modeling the vertical and cross-hemispheric transport of dust (Haberle et al., 2017;Newman et al., 2002b), as also described in Section 1.2. A reduction in the vertical transport of dust typically leads to less dust reaching, and subsequently being advected by, the mean meridional circulation. This, in turn, gives weaker correlation between surface dust deposition and topography (Newman et al., 2005).
The dust loss band around 30°S, both for simulations IC and PC, is governed by dust devil lifting, although it encompasses the wind stress lifting region N-HB. Note, however, that N-HB only covers a tiny fraction of the overall spatial extent of this band, everywhere else the dust devil lifting is the driver of dust loss. According to Ogohara and Satomura (2011), two factors are particularly important for the exchange of dust between different model grid cells: shortwave heating and vertical shear of horizontal winds, which give rise to the vertical expansion and tilt of the dust columns, respectively. Here, both may be a main factor. As also pointed out in Gebhardt et al. (2020a), the band of dust loss at about 20°S-30°S follows the low-level westerly near-surface zonal jet described in Joshi et al. (1997) and seen in Figure 4, which is known to be accompanied by strong wind shear (Forget et al., 1999;Mulholland et al., 2016). As follows from Figures 4a  and 4b, the near-surface zonal jet is indeed in overall similarity between IC and PC, with differences being rather of local nature.
The dust loss band around 30°N, both for simulations IC and PC, is zonally asymmetric and confined to within latitudes of about 20°N and 50°N. Near to its northern edge, this band includes the previously identified dust lifting regions N-OM and N-AT and, further to the south, E-EM. And, again here, large parts of the band are governed by dust devil lifting, as it is spatially too widespread to be made up of wind stress lifting only. Indeed, strong vertical wind shear can be expected close to wind-stress-related dust lifting regions such as N-OM, N-AT, and E-EM. Compliant with that, there is a near-surface zonal jet for latitudes of about 50°-60°N in Joshi et al. (1997), their Figure 1. This jet follows a general direction from west to east and has diversions to the south because of the wavenumber two topography. It is strongest around N-OM, N-AT, and E-EM and in the northern hemisphere winter/southern hemisphere summer season. The coincidence of N-OM, N-AT, and E-EM with the near-surface zonal jet was already noted in Gebhardt et al. (2020a). In the latter study, the surface dust gain/loss pattern of simulation IC was also compared against the surface dust cover map of Ruff and Christensen (2002). In addition, it was pointed that the pattern of surface dust lifting and deposition would certainly change significantly, if assuming a finite instead of an infinite surface dust cover.

2° × 2° Versus 0.5° × 0.5° Horizontal Resolution
The model topography of the 2° × 2° and 0.5° × 0.5° simulations is shown in Figures 1a and 1b, respectively. As expected, the cratering of the Mars surface is resolved in great detail for 0.5° × 0.5°. Heavily cratered terrains in the southern part of the dichotomy are opposed by only few and small craters to its north. Also, the steep slopes of the Valles Marineris are remarkably well emphasized at 0.5° × 0.5° resolution. Compared to that, the 2° × 2° resolution gives a relatively coarse representation of the surface cratering and Valles Marineris. Moreover, the topography of volcanoes is resolved in greater detail at the higher resolution. The contrast between the strong slopes at the upper part of the volcanoes and gentler slopes at their lower flanks is more pronounced for 0.5° × 0.5°. For instance, this becomes evident when comparing Arsia Mons, which is the southernmost of the three Tharsis volcanoes, between Figures 1a and 1b.
From the output of the model run with 0.5° × 0.5°, simulation PH, the total global mass of wind stress lifted dust is inferred from the first-to-last-sol of the 1 MY presented here and multiplied by a factor of 10. Thus, a 10-MY-estimate of the global wind stress lifting is obtained for simulation PH, in the following denoted as t ph . This estimate exceeds its counterpart for simulation PC, introduced as t pc in Section 3, by a factor of 8.0: It is worthwhile noting that this is considerably higher than the up to factor 2-3 differences in Toigo et al. (2012), their Table 3 (the latter is based on comparing the dust lifting between their 2° × 2° and 5° × 5° model runs; their 0.5° × 0.5° model run is not part of this table). Despite our above factor 8.0 difference, we compare simulations PH and PC based on their regional percentage contributions to the global wind stress lifting in Figures 3b and 3c. In that sense, it appears that simulation PH has roughly reproduced the dust lifting regions N-OM, N-HB, and S-HB. They have a significant percentage contribution both in Figures 3b  and 3c, not by orders of magnitude different from one another. For simulation PH, the dust lifting region S-CP is much stronger and accounts for about one third of the global wind stress lifting. This is governed by the red and yellow points at the eastern flank of Tharsis and the Valles Marineris outflow channels and slopes in Figure 3c. As follows further from comparing Figures 3b and 3c, the dust lifting regions N-AT and W-AB are of minor importance compared to the other dust lifting regions, both for simulations PC and PH. Interestingly, the percentage contribution of E-EM decays from 5.3% in the low-resolution PC simulation ( Figure 3b) to 0.1% in the high-resolution PH simulation (Figure 3c). Simulation PH spans just 1 MY, but simulations are largely repeatable from year to year when using a prescribed dust scenario. It is thus unrealistic that E-EM may be subject to strong interannual variability with significantly more surface dust lifting in other MYs. Given that, E-EM has transformed from an intermediate dust lifting region for simulation PC to a very minor dust lifting region for simulation PH.
It is important to note that increasing the model resolution may also lead to reduced wind speeds at a given location. While a higher resolution may result in stronger winds by better capturing the topographic slopes, it may also-depending on how these additional slope winds interact with larger-scale flows-result in weaker winds locally. This arises from the fact that destructive interaction between the different wind regimes is also possible in certain cases. For example, the local daytime slope winds at 0.5° × 0.5° may be in the opposite direction to the mean meridional circulation in southern summer, or in the opposite direction to regional slope flows. If so, rather than increasing the wind speed, the addition of the slope flows will be more likely to decrease wind speeds locally. In the case of Elysium Mons, a potential driver for reduced wind stress and dust lifting with increasing resolution are differences in the representation of the volcano topography between simulations PH and PC. As demonstrated below, more of its height and slopes are resolved in simulation PH than in simulation PC. This may have implications for the local topography-induced circulation. The winds at Elysium Mons are known to be due to an interplay between the topography-induced downslope flows and the large-scale atmospheric circulation. The downslope flows dominate at the upper part of Elysium Mons while the background winds have a larger impact at its lower flanks (S. W. Lee et al., 1982).
A detailed comparison of the near-surface wind vectors between simulations PH and PC is presented in Figure 7, for L S ∼ 270° and a subpart of the Elysium region, which includes the peak and north-eastern flanks of Elysium Mons and Hecates Tholus. For both simulations PH and PC, the near-surface wind vector plots are provided for 0, 8, and 16 UTC. Taking as reference a longitude of 150°E, the local time (LT) is UTC + 10, with 0, 8, and 16 UTC corresponding to 10, 18, and 2 LT, respectively. Overall, the near-surface winds are similar between simulations PH and PC in the considered region. This overall consistency is in line with the lower flanks of Elysium Mons, which have rather gentle slopes, being already well captured by simulation PC at 2° × 2° resolution. Nonetheless, there are exceptions from that at (1) the peak region of Elysium Mons and (2) Hecates Tholus, both of which are significantly better resolved at 0.5° × 0.5° resolution. In line with the height contours in Figure 7, the summit of Elysium Mons is at about 8,500 m for simulation PC (i.e., 2° × 2° resolution) and about 12,900 m for simulation PH (i.e., 0.5° × 0.5° resolution), which is about 4,400 m higher. Hecates Tholus is reasonably resolved in simulation PH, but largely uncaptured in simulation PC. At 8 UTC (18 LT), the winds are relatively weak for simulation PC in Figure 7e. For simulation PH and in the right part of Figure 7b, relatively weak local winds are overlaid by the global atmospheric circulation that is known to be from a north to north-easterly direction from wind streak patterns (S. W. Lee et al., 1982). At 0 and 16 UTC (10 and 2 LT) and for simulation PH, Figures 7a and 7c, there is clockwise air flow around Elysium Mons that is known to result from the Coriolis deflection of downslope winds (S. W. Lee et al., 1982). This is accompanied by relatively strong winds (up to 20 m/s) at, north, and south of the summit of Elysium Mons. The region to the immediate east of the peak of Elysium Mons, which makes up the dust lifting region E-EM, exhibits pronounced winds in excess of 10 m/s in simulation PC at 0 and 16 UTC (10 and 2 LT) (Figures 7d and 7f). By contrast, and for simulation PH, there are moderate upslope winds at 0 UTC (10 LT) (Figure 7a), and only moderate downslope winds at 16 UTC (2 LT) ( Figure 7c). As it appears, strong westerly winds are typically seen close to the summit of Elysium Mons for simulation PH, but because of being about 4,400 m higher than in simulation PC, there is a wind shadow to its immediate east. This is in line with the fact that the dust lifting region E-EM includes one of the top dust contributing model grid points (colored yellow in Figure 3b) in simulation PC, but this is not reproduced by simulation PH. In the latter, E-EM is only a minor dust lifting region, as described before. A solar longitude of L S ∼ 270° was selected for this comparison as the dust lifting at E-EM in simulation PC is confined to L S 180°-360°, that is, during the dust storm season, as follows from comparing (and zooming in on) Figures 2d and 2e.
The dust lifting in region E-EM undergoes a transition from making a significant contribution to the total wind stress lifted dust in simulation PC to being of negligible importance in simulation PH. The latter is consistent with the observation that dust storm texture and, relatedly, active surface dust lifting are not favored at Elysium Mons and Olympus Mons, as reported in Guzewich et al. (2015). In the referred paper, the authors present a climatology of textured dust storms, that is, dust storms indicative of ongoing surface dust lifting, based on MGS/MOC imagery. Also, their results are compared against MarsWRF simulations at 2° × 2° horizontal resolution setup with the MCD MGS scenario, following Toigo et al. (2012). This shows that high model wind stress around Elysium Mons and Olympus Mons is opposed by low occurrence of textured dust storms in the region. The authors discuss that this may arise from not favorable local surface conditions and dynamical processes of dust lifting. To repeat, the results of Figure 7 show that the dust lifting region E-EM actually features lower wind speeds in simulation PH compared to PC, due to the shadow effect of a better resolved mountain peak. In this regard, simulation PH outperforms simulation PC in the representation of the dust lifting in this region. Conversely, S-CP makes a negligible contribution to the total wind stress lifted dust in PC but is of key importance in PH, as seen in Figures 3b and 3c. Reminiscent of that, Guzewich et al. (2015) noted higher dust storm frequency in the Chryse Planitia storm track (15°S-30°N, 20°-60°W), which has a clear overlap with S-CP.
The zonal mean circulation and temperature for simulations PC and PH are shown in Figure 8 for L s 45°-135° (sols 94-290) and in Figure 9 for L s 225°-315° (sols 447-589). Around the boreal summer solstice, the winds are weak in the northern hemisphere, while the wintertime polar jet is the dominant feature in southern hemisphere's atmospheric circulation, with zonal wind speeds increasing with altitude. As expected from dynamical considerations, the polar jet coincides with the region of the strongest meridional temperature gradient, with the associated westerly winds extending down to lower levels. Given that both simulations PC and PH are driven by the same prescribed dust scenario, Figures 8a and 8b are rather similar. Still, there are some discernible differences. For example, directly to the north (south) of 60°S, the 5 m s −1 (10 m s −1 ) isoline of the zonal wind speed is further away from the surface in simulation PH when compared to simulation PC. It is estimated that the near-surface zonal wind speed of the southern polar vortex may be lower in simulation PH by up to about 50% compared to simulation PC. This may arise because of differences in the representation of the topography, Figure 1, and associated impact on the low-level atmospheric circulation. In fact, the vertical wind component is of higher strength in simulation PH at lower levels mostly southward of 40°N, where the terrain on the planet is rougher as seen in Figure 1. By and large, the Hadley Cell has comparable strength and structure in the two runs. Figure 9 is as Figure 8 but for the boreal winter solstice. As for the other solstice, the results for simulations PC and PH are broadly similar, with notable discrepancy being at lower levels southward of 30°N, consistent with differences in the representation of the topography.
It is also of interest to compare Figure 9 with the correspondent Figure 2 of Toigo et al. (2012). Due to differences in the vertical model grid, large parts of our model levels are above 1 hPa, and hence the comparison will be restricted to the bottom of the Toigo et al. (2012) results. For the PH simulation, for which Toigo et al. (2012) shows results of the same resolution of 0.5° × 0.5° and around L s 270°, the temperature distribution is basically similar, with the highest values of or in excess of 250 K just above the surface at about 60°S. The return flow of the Hadley circulation also takes place above roughly 4 hPa.

Summary and Conclusion
This study is based on the intercomparison of three MarsWRF model simulations: (1) Interactive dust at 2° × 2° horizontal resolution, simulation IC (interactive coarse resolution); (2) is as (1) but with a prescribed dust scenario, simulation PC (prescribed coarse resolution); and (3) is as (2) but with 0.5° × 0.5° horizontal resolution and a smaller integration time step, simulation PH (prescribed high resolution).
GEBHARDT ET AL.
10.1029/2020JE006672 20 of 23 Gebhardt et al. (2020a) already discussed the ability of simulation IC to produce interannual variability of GDEs. Several major dust lifting regions were identified therein such as south of Chryse Planitia (S-CP), on the northern and southern slopes of the Hellas Basin (N-HB and S-HB), north of Olympus Mons (N-OM), and east of the peak region of Elysium Mons (E-EM). The total mass of the surface dust lifted by the wind stress from the first-to-last-sol is compared between simulations IC and PC, both globally and for major dust lifting regions separately. The surface wind stress and dust lifting in simulation IC are in excess of simulation PC both globally and regionally. This is consistent with the presence/absence of radiative-dynamic feedbacks of dust in simulation IC/PC, both on a local/regional and planetary scale.
The characteristics of region E-EM are clear evidence that there is an increase in dust storm strength, duration, and number from simulation PC to IC. Assuming the same for regions S-HB and N-OM, the excess factors in surface dust lifting of 1.3 and 1.4, respectively, provide an upper estimate/constraint on the increase of dust storm strength, duration, and number from simulation PC to IC, that is, from prescribed to fully radiatively active dust. For region S-CP, we demonstrate that positive local-scale feedback between enhanced dust opacity and wind strength is a key factor. This further stresses that the advection of dust is crucial regarding the modeling of dust-radiation feedback effects. Moreover, in the case of a GDE, there are clear signs of an increased mean zonal and meridional circulation. This is consistent with strong surface dust lifting at region N-HB in some years of the simulation IC.
Regarding the dynamics of Mars dust storms, a certain feature stands out which deserves special attention. The patterns of surface dust change for simulations IC and PC both exhibit strong dust gain inside the Hellas Basin, with a high degree of similarity. In addition, the nearby dust lifting region (S-HB) is of similar strength in both model runs. That is, in turn, reminiscent of the idealized dust lifting experiments by Ogohara and Satomura (2008), and the observed precursory phase of the 2001 GDE (Cantor, 2007;Strausberg et al., 2005). This sheds more light on the fact that simulation IC could not reproduce such major dust storm activity early in the dust storm season. As pointed out in Gebhardt et al. (2020a), the most promising candidate is a midlevel atmospheric temperature (T15) rise of about 5K at L S 190-195 (their Figure 3a, orange curve for Year 14). That is, however, only a precursor of a potential GDE according to the dust storm classification scheme of Newman et al. (2019). In line with the latter, a potential GDE requires a temperature rise of 10 K. The lack of major dust storm activity around L s 180 in simulation IC may arise from the fact that dust clouds inside the Hellas Basin, originating from S-HB, cannot trigger additional dust lifting at N-HB and related expansion of dust clouds, which would be reminiscent of the onset of the 2001 GDE. That is because dust lifting at S-HB and N-HB never coincide in time. In simulation IC, dust lifting at S-HB occurs around L S 180 while that at N-HB starts clearly later in the dust storm season, see Gebhardt et al. (2020a), their Figure 10. This is a valuable insight, as such deficiency might be overcome with future model versions that include additional negative feedback mechanisms such as dust particle scavenging by radiatively active water ice clouds.
In simulation PH, the MarsWRF model was run globally and with prescribed dust at an extremely high horizontal resolution of 0.5° × 0.5° for 2 MYs (the first MY was for model spin-up, the second MY for scientific analysis). We found that the global surface dust lifting by the wind stress for simulation PH exceeded that of simulation PC by factor of about 8.0. Irrespective of that, some major dust lifting regions may also diminish with increasing resolution. That is the case for E-EM, which is a major/intermediate (minor) dust lifting region for simulation IC/PC (PH). This is governed by strong differences in the topography representation of the peak region of Elysium Mons and the associated simulation of local winds between 2° × 2° and 0.5° × 0.5° horizontal model resolution. The more detailed representation of the topography in simulation PH compared to PC may also account for differences in wind vectors at the lowermost vertical model levels. By and large, the zonal mean circulation is of comparable strength and has a similar structure in both experiments at both the boreal summer and winter solstices.

Acknowledgments
Once again, our warmest thanks go to the PlanetWRF development team for providing the MarsWRF model free of charge to us and their proactive attitude in general. We would also like to thank two anonymous reviewers and the Associate Editor Dr Claire Newman for their several detailed and insightful comments and suggestions that helped to significantly improve the quality of the paper. We would like to acknowledge the support of this work by funding from the United Arab Emirates University (UAE University). Also, we are deeply grateful to High-Performance Computing, Division of Information Technology, UAE University, for the valuable access to the computational resources required for this work. We thank IT engineers Asma AlNeyadi, Anil Thomas, and Nithin Damodaran for their professional assistance and support in technical questions. M.-P. Z. has been partially funded by the AEI (MDM-2017-0737, Unity of Excellence "María de Maeztu" -Centro de Astrobiología (CSIC-INTA)) and the Spanish Ministry of Science and Innovation (PID2019-104205GB-C219). Finally, we declare that there are no real or perceived conflicts of interests for any author.