Characteristics of Dust Devils in Two Pre-Selected Landing Regions of the Tianwen-1 Mission—Comparing Observations and Predictions Using Numerical Model

: The spatial and temporal distribution of dust devils (DDs) in the two pre-selected landing regions (ZA and ZB) of the Tianwen-1 mission in southern Utopia Planitia have been investigated by using images from the Context Camera (CTX) of the Mars Reconnaissance Orbiter (MRO). From the images of the regions in 8 Martian years, no DD was found in ZA, while 77 DDs were found in ZB. The observed DDs are mainly distributed in the northeastern part of ZB. The temporal variation in the observed DDs shows a prominent two-peak pattern in their local early spring and late summer. The size and height of the observed DDs have also been evaluated from the images, and they show a similar temporal variation as the occurrence. To investigate the possible conditions pertinent to these observed patterns of DD distribution, some analysis based on the thermodynamic theory of heat engines was performed using the output of the Mars climate model, MarsWRF. The spatial and temporal distribution of the simulated DDs are generally consistent with the observation, with signiﬁcantly more DDs in ZB. Analysis of the model results suggests that the spatial distributions of the predicted DDs are mainly related to the distribution of sensible heat ﬂux, which, in turn, is mainly determined by the surface-to-air temperature difference. The difference in DDs between ZA and ZB (more DDs in ZB) is dominated by the difference in sensible heat ﬂux, which, in turn, is mainly related to the spatial variation of surface albedo. shows little effect on the spatial distribution of surface temperature. The near-surface air to southeast, be caused thermal ad- vection. results show distribution of surface-to-air temperature the solar radiation and thermal advection, in the distribution of sensible heat flux and so the DLIF. In ZB the simulation results show that sensible heat flux in the western part is mostly contributed by the surface-to-air temperature difference, while the horizontal variation in near-surface air density seems to be In the southwest- part, the sensible heat flux is mainly contributed by the relatively large drag velocity there, which to surface roughness and that, the and sensible heat flux are relatively large in the southeastern part, the there is nearly zero. in our model, the occurrence of dust lifting by DDs requires three conditions as the surface sensible heat flux the southeast of ZB the for the occurrence of the surface-to-air temperature difference does the surface-to- air temperature is the most significant in determining the distribution of sensible heat flux and, the DLIF. spatial distribution of sur- face temperature, near-surface air temperature, net short-wave flux, and surface proper- distribution of surface temperature controlled by the short-wave flux the


Introduction
Dust devils (DDs) are thermally driven columnar vortices and the smallest-scale dustlifting events on Mars. Surface dust can be lifted into the atmosphere due to the strong surface wind and the suction effect of the low-pressure vortex core of a DD [1]. The dustlifting process of DDs has been parameterized in many Mars climate models. Model results and observations suggest that DDs may make a significant contribution to the dust loading in the Martian atmosphere and so the observed annual dust cycle [2][3][4][5][6], which, in turn, will affect the thermal structure and atmospheric circulation of Mars.
DDs were first identified on Mars by the Viking Orbiter [7] and were first measured by the Viking Lander [8], and subsequently observed by other missions, such as Mars Global Surveyor (MGS)/Mars Orbiter Camera (MOC) [9,10], Mars Express (MEX)/High-Resolution Stereo Camera (HRSC) [2], Mars Reconnaissance Orbiter (MRO)/Context Camera (CTX) and High-Resolution Imaging Science Experiment (HiRISE) [11][12][13]. In general, the appearance of DDs can be different in different types of images. From a rover (e.g., Opportunity, Curiosity), DDs always show a narrow base and a broader top. From orbiters, a typical active dust devil (ADD) is basically composed of a small bright cloud and a long The main job of the CTX camera is to help provide the wider context for the highresolution observations of HiRISE and CRISM (Compact Reconnaissance Imaging Spectrometer for Mars) at the MRO. Images from CTX are used in this study because of the relatively large spatial coverage (a typical image can be as wide as 30 km and as long as 160 km or more), long orbital time (from 2006 to now), and a pretty good resolution (~6 m/pixel). However, due to the limitation of the orbit, all images were taken at a local time (LT) around 15

Dimension of Dust Devils
In addition to the occurrence of DDs, we also surveyed their sizes. Observed ADDs by orbiters are usually tens of meters to hundreds of meters in diameter, and hundreds of meters to thousands of meters in height, with a narrow base and broader top. A typical ADD is usually composed of a small bright cloud with long tapered shadows [7]. Although some small terrains may also show this feature, it can be ruled out by the images in the same location at a different time. Topography (shading) of regions ZA and ZB and the spatial distribution of ADDs (white dots) from CTX observations. Data of topography obtained from https://astrogeolo$-$gy.usgs.gov/search/map/Mars/Topography/HRSC_MOLA_Blend/ Mars_HRSC_MOLA_BlendDE$-$M_Global_200mp_v2 (accessed on 20 July 2021).
From all available CTX images, a total of 258 images in ZA and 266 images in ZB can be obtained, which were then used to identify DDs from September 2006 to February 2021 (spanning~8 MYs). As shown in Table 1, no ADD was found in the 258 images of ZA. On the other hand, 77 ADDs were found in ZB, and are mainly distributed in the northeast part of ZB (white dots in Figure 1). This result is interesting because, by intuition, ZA and ZB should have a similar number of DDs, since they are at similar latitudes and so receive a similar amount of downward solar radiation. In addition, the altitudes of ZA and ZB are similar (~4.28 km). In the light of some previous literature, it is speculated that different atmospheric conditions, topography, or surface properties of these two areas may explain the different distribution of ADDs. The terrain of ZA is relatively flat, and there are no large terrain undulations around it. In contrast, ZB is located at the northwestern foot of the Elysium Mons, and the terrain is more complex, and the surface roughness and slope in the southeast of ZB are higher. It is interesting to note that the geographical situation of ZB is similar to that of the Amazonis Planitia. Both regions have a mountain in their east. As mentioned before, Amazonis Planitia is the region particularly active in DD occurrence.

Dimension of Dust Devils
In addition to the occurrence of DDs, we also surveyed their sizes. Observed ADDs by orbiters are usually tens of meters to hundreds of meters in diameter, and hundreds of meters to thousands of meters in height, with a narrow base and broader top. A typical ADD is usually composed of a small bright cloud with long tapered shadows [7]. Although some small terrains may also show this feature, it can be ruled out by the images in the same location at a different time.
The CTX images are stretched by the ESRI ArcGIS software in ArcMap to increase the contrast, so that the ADDs can be easier to identify manually. To measure the diameters and shadow lengths of the ADDs, a simple cylindrical coordinate system of Mars is applied to all images imported into the ArcMap. Their measurement is not straightforward because of their variations in morphology during movement. At the base of the ADD, there is always a "dust skirt" that is usually larger than the vortex itself and so may overestimate the diameter of an ADD [21,22]. In addition, sometimes it is difficult to determine the edge of the shadow because of the variations in solar incidence angle. Besides, the DDs tend to lean in the direction of the ambient wind, which will affect the accurate measurement of diameter and shadow length.
Because of these uncertainties, a uniform standard is set for the measurements to minimize the subsequent analysis error (averaged error of~20% in diameter, and~11% in the height of ADDs). That is, the standard should be to draw a circle on the brightest part of the dust cloud and try to include the area where the brightest pixels are densely gathered. The diameter of the circle is then considered to be the diameter of the ADD. The distance from the center of the circle to the edge of the shadow line is then considered as the length of the shadow (see the example in Figure 2a,b). If the dust cloud is in a cone shape (probably due to the inclination of ADD by ambient wind [23]), the intersection of the shadow and the bright dust cloud is considered as the starting point for shadow measurement (see the example in Figure 2c,d). The height of an ADD is calculated by dividing the shadow length by the tangent of the local solar incidence angle.
the contrast, so that the ADDs can be easier to identify manually. To measure the diameters and shadow lengths of the ADDs, a simple cylindrical coordinate system of Mars is applied to all images imported into the ArcMap. Their measurement is not straightforward because of their variations in morphology during movement. At the base of the ADD, there is always a "dust skirt" that is usually larger than the vortex itself and so may overestimate the diameter of an ADD [21,22]. In addition, sometimes it is difficult to determine the edge of the shadow because of the variations in solar incidence angle. Besides, the DDs tend to lean in the direction of the ambient wind, which will affect the accurate measurement of diameter and shadow length.
Because of these uncertainties, a uniform standard is set for the measurements to minimize the subsequent analysis error (averaged error of ~20% in diameter, and ~11% in the height of ADDs). That is, the standard should be to draw a circle on the brightest part of the dust cloud and try to include the area where the brightest pixels are densely gathered. The diameter of the circle is then considered to be the diameter of the ADD. The distance from the center of the circle to the edge of the shadow line is then considered as the length of the shadow (see the example in Figure 2a,b). If the dust cloud is in a cone shape (probably due to the inclination of ADD by ambient wind [23]), the intersection of the shadow and the bright dust cloud is considered as the starting point for shadow measurement (see the example in Figure 2c,d). The height of an ADD is calculated by dividing the shadow length by the tangent of the local solar incidence angle.  (Table 1). (c,d) are from CTX image p19_008503_2 058_xn_25n227w.

Seasonal Variations in the Observed DDs
As can be seen from Figure 3 (left), the number of ADDs shows a prominent twopeak pattern in the temporal variation, with peaks in the local early spring (Ls: 0-30°) and late summer (Ls: 150-180°). This two-peak pattern of ADDs occurrence has also been found in other regions in some previous studies, such as the Argyre Planitia [17], and the  (Table 1). (c,d) are from CTX image p19_008503_2058_xn_25n227w.

Seasonal Variations in the Observed DDs
As can be seen from Figure 3 (left), the number of ADDs shows a prominent two-peak pattern in the temporal variation, with peaks in the local early spring (L s : 0-30 • ) and late summer (L s : 150-180 • ). This two-peak pattern of ADDs occurrence has also been found in other regions in some previous studies, such as the Argyre Planitia [17], and the Amazonis Planitia [16]. Combining with the results from the previous research of [20], based on dust devil tracks (DDTs) in the same study area (Figure 3, right), it can be seen that the sum of ADDs and DDTs still show a two-peak pattern, as before. However, the first peak came a little bit late.
Amazonis Planitia [16]. Combining with the results from the previous research of [20], based on dust devil tracks (DDTs) in the same study area (Figure 3, right), it can be seen that the sum of ADDs and DDTs still show a two-peak pattern, as before. However, the first peak came a little bit late.

Numerical Simulations
From the observational results discussed above, some questions are worth exploring. For instance, why are DDs greater in ZB compared with ZA? Why are there two peaks for the occurrence of DDs in the temporal distribution? What are the possible factors that determine the spatial distribution of DDs in ZA and ZB (e.g., ADDs are mainly in the northeastern part of ZB)? To explore these questions, we conducted some analysis in ZA and ZB based on numerical simulations with the Mars climate model MarsWRF. The approach is basically similar to that in Newman et al. [19] but the analyzed area is much larger compared with the Gale crater in Newman et al. [19]. based on dust devil tracks (DDTs) in the same study area (Figure 3, right), it can be seen that the sum of ADDs and DDTs still show a two-peak pattern, as before. However, the first peak came a little bit late. The median and mean values of the measured dimension of ADDs ( Figure 4) show a similar pattern of seasonal variation as the ADDs number ( Figure 3). During the two-peak periods, the sizes of ADDs reach their maximum when the number of occurrences reaches their maximum. In the second half of the year, it is clear that the size of ADDs is suppressed.

Numerical Simulations
From the observational results discussed above, some questions are worth exploring. For instance, why are DDs greater in ZB compared with ZA? Why are there two peaks for the occurrence of DDs in the temporal distribution? What are the possible factors that determine the spatial distribution of DDs in ZA and ZB (e.g., ADDs are mainly in the northeastern part of ZB)? To explore these questions, we conducted some analysis in ZA and ZB based on numerical simulations with the Mars climate model MarsWRF. The approach is basically similar to that in Newman et al. [19] but the analyzed area is much larger compared with the Gale crater in Newman et al. [19].

Numerical Simulations
From the observational results discussed above, some questions are worth exploring. For instance, why are DDs greater in ZB compared with ZA? Why are there two peaks for the occurrence of DDs in the temporal distribution? What are the possible factors that determine the spatial distribution of DDs in ZA and ZB (e.g., ADDs are mainly in the northeastern part of ZB)? To explore these questions, we conducted some analysis in ZA and ZB based on numerical simulations with the Mars climate model MarsWRF. The approach is basically similar to that in Newman et al. [19] but the analyzed area is much larger compared with the Gale crater in Newman et al. [19].

Numerical Model MarsWRF
MarsWRF is a Martian atmospheric model, which can perform multi-scale climate simulations of Mars [24,25]. MarsWRF includes most of the physical processes for simulating the Martian climate, such as the Carbon Dioxide cycle [26]. Short-wave and long-wave radiation is evaluated by the Wide Band Model (WBM), as described in Richardson et al. [24] and Toigo et al. [25]. For dust radiative forcing, the so-called interactive approach is implemented in the simulations. In this approach, dust lifting is dependent on the atmospheric conditions (e.g., surface wind stress), while the suspending dust particles may change the atmospheric radiation and so the atmospheric conditions. The configuration in the MarsWRF is basically similar to that used in our previous studies (e.g., [27,28]), and it has Remote Sens. 2022, 14, 2117 7 of 17 been shown in these studies that the interactive dust approach could reasonably reproduce the observed dust climate on Mars.
Two simulations with different horizontal resolutions were performed in this study. The relatively low-resolution simulation (global domain, resolution: 3 • ) was used to study the seasonal variations in a Martian year. Another simulation has a higher resolution, which is used to study the spatial variations in the predicted DDs. In this simulation, a nested domain is used. The inner domain (Domain 2) has a resolution of about 1.67 • (about 100 km in the equatorial region), which covers the two study regions, ZA and ZB, while the parent global domain has a resolution of about 5 • (about 300 km). Because of the expensive computing cost, nested simulations were not performed for a whole Martian year but for a particular period (high season of DDs).
There are two parameterization schemes of dust-lifting processes in MarsWRF, one for DDs and the other for wind stress lifting [29]. Wind stress lifting occurs in cases when stress is greater than the threshold (0.043). See Richardson et al. [24] and Xiao et al. [28] for more details.
For DDs, most previous studies (e.g., [3,19]) interpreted them as convective vortices and the convective heat engine theory [30], which is usually used to describe or parameterize the DD process. In this theory, the strength of DDs (and so the associated dust-lifting flux) is proportional to the sensible heat flux and the thermal efficiency, as follows: where the strength Λ is the dust-lifting flux by dust devils (DLIF), η is thermal efficiency, given as Equation (2), and F s is sensible heat flux, given as Equation (4).
In the above, T top is the temperature at the top of the planetary boundary layer (PBL), T sur f is the surface temperature. The PBL top is defined by the PBL depth. There are certain conditions for DDs to occur in the model (F s > 0, T surf − T air > 27 K, T surf > 270 K), generally referred to as the condition of atmospheric instability near the surface.
The sensible heat flux is calculated by the bulk formula: where C d is the exchange coefficient, depending on the stability of the near-surface atmosphere, ρ is the near-surface air density, u * is drag velocity, T sur f is surface temperature, and T air is the near-surface air temperature (above the surface about 1.5 m).

Seasonal Variations in Predicted DDs
The predicted DDs from MarsWRF are reflected by DLIF. The time series of DLIF in ZB ( Figure 5) shows that the temporal pattern of predicted DDs is generally consistent with observations. Obviously, the DLIF in ZB is higher than that in ZA in the first half of the year. In addition, in both ZA and ZB, there are two peaks of DLIF during their local early spring and late summer, which are in good agreement with the observations. After entering the dust season (the second half of the year), the generation of dust devils is very rare, which is also consistent with our observations and similar to those in previous studies [31]. To further understand why there are more DLIFs in ZB and the two-peak structure in temporal distribution, we analyzed the corresponding temporal variation in sensible heat flux and thermal efficiency, based on Equation (1). heat flux, and thermal efficiency in Equation (1)). After entering the dust season in the second half of the year, the thermal efficiency and sensible heat flux drop sharply and so does the DLIF. Apparently, their net short-wave flux at surface and surface temperature ( Figure 8) both reach their maximum at the two peak times (early spring and late summer). The results indicate that the two-peak structure in temporal distribution is related to the seasonal variation of solar radiation, and the pause between the two peaks corresponds to the period with relatively weaker solar radiation.  Overall, in the first half of the year, the sensible heat flux and thermal efficiency in ZB ( Figure 6) are always higher than those in ZA, and the difference in sensible heat flux between ZA and ZB (~5 W m −2 ) is generally larger than that of the thermal efficiency. Based on the DD theory of Rennó et al. [30], it can be inferred that the difference in DLIF between ZA and ZB is mainly determined by sensible heat flux, which, in turn, is mainly related to the surface-to-air temperature difference (~2 K, Figure 7). The drag velocity plays a minor role in the difference (~0.2 m s −1 ). Besides, the mean albedo in ZA is about 0.229, while the mean albedo in ZB is around 0.205. This difference in albedo is also reflected in the net short-wave flux at the surface (short-wave flux at surface minus the reflected short-wave flux), which in ZB, is always stronger than ZA in their local spring and summer ( Figure 8). By intuition, ZA and ZB are at similar latitudes and so should have similar solar radiation at the surface. The above results suggest that more DDs in ZB are mainly due to the relatively low albedo there.
The two-peak pattern can generally be seen in DLIF and its contributors (i.e., sensible heat flux, and thermal efficiency in Equation (1)). After entering the dust season in the second half of the year, the thermal efficiency and sensible heat flux drop sharply and so does the DLIF. Apparently, their net short-wave flux at surface and surface temperature ( Figure 8) both reach their maximum at the two peak times (early spring and late summer). The results indicate that the two-peak structure in temporal distribution is related to the seasonal variation of solar radiation, and the pause between the two peaks corresponds to the period with relatively weaker solar radiation.

Spatial Distribution of Predicted DDs
This section focuses on the spatial distribution of DDs during their local summer (L s = 150-180 • ), which is the time when solar forcing and the number of occurrences of DDs both reach their maximum, as discussed above. In the analysis of spatial distribution, we considered the day-time mean of the model results over the time period betweeñ 11:00 and~15:00 local time in the corresponding regions of ZA and ZB. This is the period when the dust lifting by dust devils occurs in the model simulations. Figure 9 (first column) clearly shows the completely different spatial distributions of predicted DDs in ZA and ZB. Predicted DDs are relatively active in ZB and are basically concentrated in the northwest part of ZB. There is a slight difference between the predicted and observed spatial distribution of DDs. This difference may be caused by the insufficient sample size or the inaccurate prediction by the model, which will be further discussed. Remote Sens. 2022, 14, x 9 of 18

Spatial Distribution of Predicted DDs
This section focuses on the spatial distribution of DDs during their local summer (Ls = 150-180°) , which is the time when solar forcing and the number of occurrences of DDs both reach their maximum, as discussed above. In the analysis of spatial distribution, we considered the day-time mean of the model results over the time period between ~11:00 and ~15:00 local time in the corresponding regions of ZA and ZB. This is the period when the dust lifting by dust devils occurs in the model simulations. Figure 9 (first column) clearly shows the completely different spatial distributions of predicted DDs in ZA and ZB. Predicted DDs are relatively active in ZB and are basically concentrated in the northwest part of ZB. There is a slight difference between the predicted and observed spatial distribution of DDs. This difference may be caused by the insufficient sample size or the inaccurate prediction by the model, which will be further discussed.
According to Equation (1), we plotted the spatial distribution of sensible heat flux ( Figure 9, second column) and thermal efficiency (Figure 9, third column) in their local summer (Ls: 150-180°). In ZA, the predicted DDs are mainly distributed in the north and this distribution is mainly determined by the sensible heat flux. The thermal efficiency in ZA has an opposite trend in the south-north direction. As mentioned above, the predicted DDs in ZB are mainly distributed in the northwest. In addition, the relatively large DLIF in the northwest of ZB is dominated by the sensible heat flux. The thermal efficiency may also play a positive role in the greater DLIF in the northwest of ZB, although the effect is much lower than that of the sensible heat flux. Notice that the DLIF in the southeast of ZB is nearly zero, while the sensible heat flux and thermal efficiency in the southeast are relatively large. This could be caused by the conditions for the occurrence of DD in the DD parameterization (discussed in Section 3.1), which will be further explained later. Overall, from the model results, we can see that the spatial distribution of DLIF both in ZA and ZB is dominated by sensible heat flux. To further understand the possible environmental conditions that contribute to the spatial distribution of sensible heat flux, we analyzed their contributors in ZA ( Figure 10) and ZB (Figure 11), according to Equation (4).

Net shortwave flux at surface
Surface Temperature (Tair)   According to Equation (1), we plotted the spatial distribution of sensible heat flux ( Figure 9, second column) and thermal efficiency (Figure 9, third column) in their local summer (L s : 150-180 • ). In ZA, the predicted DDs are mainly distributed in the north and this distribution is mainly determined by the sensible heat flux. The thermal efficiency in ZA has an opposite trend in the south-north direction. As mentioned above, the predicted DDs in ZB are mainly distributed in the northwest. In addition, the relatively large DLIF in the northwest of ZB is dominated by the sensible heat flux. The thermal efficiency may also play a positive role in the greater DLIF in the northwest of ZB, although the effect is much lower than that of the sensible heat flux. Notice that the DLIF in the southeast of ZB is nearly zero, while the sensible heat flux and thermal efficiency in the southeast are relatively large. This could be caused by the conditions for the occurrence of DD in the DD parameterization (discussed in Section 3.1), which will be further explained later. Overall, from the model results, we can see that the spatial distribution of DLIF both in ZA and ZB is dominated by sensible heat flux. To further understand the possible environmental conditions that contribute to the spatial distribution of sensible heat flux, we analyzed their contributors in ZA ( Figure 10) and ZB (Figure 11), according to Equation (4).

Sensible Heat Flux and Its Contributors
The predicted sensible heat flux in ZA ( Figure 10) decreases from northwest to southeast, and so does the surface-to-air temperature difference (Tsurf − Tair). The variation in near-surface air density (~0.001 kg m −3 ) and drag velocity (~0.2 m s −1 ) shows a similar trend but the variation is very small. In addition, the wind direction basically corresponds to the Hadley circulation and thermal tide during their local summer. This indicates that the spatial distribution of sensible heat flux is mainly determined by the surface-to-air temperature difference. To better understand the surface-to-air temperature difference, the surface temperature, near-surface air temperature, net short-wave flux, and surface properties are further analyzed (Figure 12). The surface temperature decreases with increasing latitudes, and so does the net short-wave flux at the surface. On the other hand, the variation in the spatial distribution of albedo and thermal inertia is very small and, thus,

Sensible Heat Flux and Its Contributors
The predicted sensible heat flux in ZA ( Figure 10) decreases from northwest to southeast, and so does the surface-to-air temperature difference (T surf − T air ). The variation in near-surface air density (~0.001 kg m −3 ) and drag velocity (~0.2 m s −1 ) shows a similar trend but the variation is very small. In addition, the wind direction basically corresponds to the Hadley circulation and thermal tide during their local summer. This indicates that the spatial distribution of sensible heat flux is mainly determined by the surface-to-air temperature difference. To better understand the surface-to-air temperature difference, the surface temperature, near-surface air temperature, net short-wave flux, and surface properties are further analyzed (Figure 12). The surface temperature decreases with in-creasing latitudes, and so does the net short-wave flux at the surface. On the other hand, the variation in the spatial distribution of albedo and thermal inertia is very small and, thus, shows little effect on the spatial distribution of surface temperature. The near-surface air temperature increases from northwest to southeast, which may be caused by thermal advection. These results show that the distribution of surface-to-air temperature difference is mainly determined by the solar radiation and thermal advection, which, in turn, determines the distribution of sensible heat flux and so the DLIF. although surface temperature strongly controls the distribution of near-surface air temperature via radiative heating and sensible heat flux, the near-surface air could be mixed with the local wind (e.g., wind associated with local topography). Overall, our results suggest that the spatial distribution of sensible heat flux in ZB is mainly determined by the distribution of albedo, and the local topography may have a secondary effect.

Thermal Efficiency and Its Contributors
Although DLIF is mainly determined by the sensible heat flux in both ZA and ZB, it is still worth analyzing the thermal efficiency. According to Equations (2) and (3), thermal efficiency is proportional to the temperature difference between the surface and the top of the convective vortex. Considering the spatial distribution of thermal efficiency ( Figure  14), we can see that its pattern is dominated by the surface temperature in both ZA and ZB, and the distribution of surface temperature has been discussed before. Therefore, the distribution of thermal efficiency in ZA is mainly related to solar radiation, while it is mainly related to the distribution of albedo in ZB. In ZB (Figure 11), the simulation results show that sensible heat flux in the northwestern part is mostly contributed by the surface-to-air temperature difference, while the horizontal variation in near-surface air density seems to be secondary. In the southwestern part, the sensible heat flux is mainly contributed by the relatively large drag velocity there, which is proportional to surface roughness and altitude. Notice that, although the drag velocity and sensible heat flux are relatively large in the southeastern part, the DLIF there is nearly zero. This is because, in our model, the occurrence of dust lifting by DDs requires three conditions (Section 3.1), as mentioned above. Both the surface temperature and sensible heat flux in the southeast of ZB meet the conditions for the occurrence of DDs, while the surface-to-air temperature difference does not. Therefore, the surface-to-air temperature difference is the most significant factor in determining the distribution of sensible heat flux and, hence, the DLIF. We also analyzed the spatial distribution of surface temperature, near-surface air temperature, net short-wave flux, and surface properties in ZB ( Figure 13). The distribution of surface temperature is mostly controlled by the net short-wave flux at the surface, which is inversely proportional to the albedo. The surface temperature in the southeastern part of ZB is mostly related to the distribution of thermal inertia. Lower thermal inertia implies a faster heating rate at the surface. Besides, although surface temperature strongly controls the distribution of near-surface air temperature via radiative heating and sensible heat flux, the near-surface air could be mixed with the local wind (e.g., wind associated with local topography). Overall, our results suggest that the spatial distribution of sensible heat flux in ZB is mainly determined by the distribution of albedo, and the local topography may have a secondary effect.

Thermal Efficiency and Its Contributors
Although DLIF is mainly determined by the sensible heat flux in both ZA and ZB, it is still worth analyzing the thermal efficiency. According to Equations (2) and (3), thermal efficiency is proportional to the temperature difference between the surface and the top of the convective vortex. Considering the spatial distribution of thermal efficiency ( Figure  14), we can see that its pattern is dominated by the surface temperature in both ZA and ZB, and the distribution of surface temperature has been discussed before. Therefore, the distribution of thermal efficiency in ZA is mainly related to solar radiation, while it is mainly related to the distribution of albedo in ZB.

Tsurf − Tair Tsurf Tair
Net short wave flux at surface Albedo Thermal inertia Figure 13. Similar to Figure 12, but for ZB.

Thermal Efficiency and Its Contributors
Although DLIF is mainly determined by the sensible heat flux in both ZA and ZB, it is still worth analyzing the thermal efficiency. According to Equations (2) and (3), thermal efficiency is proportional to the temperature difference between the surface and the top of the convective vortex. Considering the spatial distribution of thermal efficiency (Figure 14), we can see that its pattern is dominated by the surface temperature in both ZA and ZB, and the distribution of surface temperature has been discussed before. Therefore, the distribution of thermal efficiency in ZA is mainly related to solar radiation, while it is mainly related to the distribution of albedo in ZB.

Summary
This study focuses on the distribution and characteristics of DDs in the two pre-selected landing sites (ZA and ZB) of the Tianwen-1 mission (ZA was finally chosen for landing). A significant difference in the number of DDs occurring in these two regions was found in observations, which was also supported by the results of the numerical simulations with the MarsWRF model. Some questions have been put forward based on the

Summary
This study focuses on the distribution and characteristics of DDs in the two pre-selected landing sites (ZA and ZB) of the Tianwen-1 mission (ZA was finally chosen for landing). A significant difference in the number of DDs occurring in these two regions was found in observations, which was also supported by the results of the numerical simulations with the MarsWRF model. Some questions have been put forward based on the observational results. In order to answer these questions, some analysis based on numerical simulations was performed from the perspective of atmospheric conditions and surface properties.
DDs are mainly concentrated in the first half of the year and there are more DDs in ZB than that in ZA. Based on the DD theory of Rennó et al. [30], our analysis suggests that the difference in dust lifting by DDs (DLIF) between ZA and ZB is mainly dominated by the difference in sensible heat flux, which is mainly caused by the surface-to-air temperature difference. Thermal efficiency plays a relatively minor role. Furthermore, ZB has a stronger net short-wave radiative flux at the surface due to the lower albedo (~0.205) at the surface than ZA (~0.229). Overall, the results of the analysis suggest that more DDs in ZB are mainly caused by the lower albedo in ZB. However, it is worth mentioning that some of the differences could also be due to the difference in surface dust availability between the two regions, but there is a lack of observation and this is beyond the scope of the present study.
Both the observed and predicted DDs have a two-peak feature in their local early spring and late summer, and are generally coincident with the time of maximum solar radiation in their local time. This suggests that the temporal variation in solar radiation due to the revolution of Mars is primarily responsible for the two-peak pattern of temporal distribution. After entering the dust season in the second half of the year, the thermal efficiency and sensible heat flux drop sharply and so do the DD activities.
The spatial distributions of DDs, both in ZA and ZB, are mainly related to the distribution of sensible heat flux and surface-to-air temperature difference. In ZA, solar radiation at the surface dominates the spatial distribution of sensible heat flux, thus, causing the particular spatial distribution of DLIF. In ZB, our results suggest that the spatial distribution of surface-to-air temperature difference and sensible heat flux are mainly determined by the albedo, and so the resulting DLIF is concentrated in the northwest. The surrounding topography has little effect on the distribution of the predicted DDs. In the southeast corner of ZB, the surface-to-air temperature difference is relatively small, which is not favorable to the occurrence of DDs.

Discussion
As mentioned in Section 3.3, the predicted DDs are mainly concentrated in the northwest of ZB, while the observed DDs are mainly distributed in the northeast of ZB. The model's prediction ability may be responsible for this discrepancy. First, the predicted DDs in the southeast are almost 0, but there are quite a few observed DDs here. According to the analysis of the sensible heat flux and its contributors, the threshold of surface-to-air temperature difference for DD occurrence is responsible for the no DLIF in the southeast. In fact, the distribution of sensible heat flux is more similar to that of the observed DDs. Besides, the distribution of surface temperature seems to be more similar to that of the sensible heat flux and so the predicted DDs, compared with the distribution of the surfaceto-air temperature difference. Therefore, the role played by surface temperature should be more emphasized in the parameterization of DDs in the model.
From the simulation results, it is worth noting that there are more predicted DDs in the south of the Elysium Mons compared with that in ZA and ZB ( Figure 15). However, from 518 CTX images (~8 MY), there are only 12 ADDs found in the south of Elysium Mons. The list of CTXs with ADDs in the south of Elysium Mons is shown in Table 2. In a similar geographical situation, there are only two DDs found in the south of the Olympus Mons, and a large number of DDs (9827) were found in the Amazonis Planitia from the global survey of DDs [10]. However, in our numerical simulations, DLIF in the south of Olympus is also greater than that in the Amazonis Planitia ( Figure 16). Nevertheless, our simulations did very well in the time series of DDs in the Amazonis Planitia ( Figure 17). The above discrepancy in DDs suggests that the present configuration of MarsWRF is likely not good enough to predict DDs in mountainous areas. There are two possible reasons for the discrepancy. One is the use of the unlimited dust source in our model, while in reality, there may be little dust on the hillside (most of which is blown away by the wind). The other one is that the parameterization of dust-lifting processes is not good enough for mountainous areas. Improving the simulation of DDs in mountainous areas is one important task of our future work. similar geographical situation, there are only two DDs found in the south of the Olympus Mons, and a large number of DDs (9827) were found in the Amazonis Planitia from the global survey of DDs [10]. However, in our numerical simulations, DLIF in the south of Olympus is also greater than that in the Amazonis Planitia ( Figure 16). Nevertheless, our simulations did very well in the time series of DDs in the Amazonis Planitia ( Figure 17). The above discrepancy in DDs suggests that the present configuration of MarsWRF is likely not good enough to predict DDs in mountainous areas. There are two possible reasons for the discrepancy. One is the use of the unlimited dust source in our model, while in reality, there may be little dust on the hillside (most of which is blown away by the wind). The other one is that the parameterization of dust-lifting processes is not good enough for mountainous areas. Improving the simulation of DDs in mountainous areas is one important task of our future work.