Some features of effective radius and variance of dust particles in numerical simulations of the dust climate on Mars

Airborne dust is an important constituent in the Martian atmosphere because of its radiative interaction with the atmospheric circulation. Dust size is one crucial factor in determining this effect. In reality dust sizes are varied; however, in numerical modeling of dust processes, dust size has usually been described by choice of a particular size distribution function, or by use of fixed values of effective radius (ER) and effective variance (EV). In this work, we present analytical expressions that have been derived to specify ER and EV for N‐bin dust schemes, based on a model‐calculated dust mixing ratio. Numerical simulations based on this approach thus would consider the effects of variable ER on the atmospheric radiation and their interaction. Results have revealed some interesting features of the dust distribution parameters, such as seasonal and spatial variation of ER and EV, which are generally consistent with some previous observational and modeling studies. Compared with the usual approach of using a fixed ER, simulation results from the present approach suggest that the variability of ER can have significant effects on the simulated thermal field of the Martian atmosphere.


Introduction
Dust particles suspending in the atmosphere of Mars have an important effect on the Martian climate due to their radiative heating and cooling effects in the atmosphere Goody, 1968, 1972). The radiative properties of dust particles are determined by their microscopic properties, such as composition (i.e. complex index of reflection), shape, and size. While both composition and shape are usually assumed to be homogenous or uniform, the size distribution can vary both spatially and temporally, so do the radiative effects. The effects of dust on radiation are basically dependent on the quantities Effective Radius (ER) and Effective Variance (EV). These two quantities generally can represent the essential effects of a atmospheric distribution of dust particles of different sizes, and thus can be derived from the underlying size distribution if it is known (Hansen and Travis, 1974). The importance of both the values and the variations of dust ER and EV to our understanding of the general circulation of Mars has attracted the attention of some previous studies (e.g. Murphy et al., 1993;Kahre et al., 2008). ER and EV have been measured by various remote-sensing observations since Mariner 9, Viking, and several subsequent missions (see reviews by Pollack et al. 1995;Tomasko et al., 1999;Dlugach et al., 2003;Smith, 2008). These retrievals of ER have required the assumption of particular types of dust size distribution, such as the Gamma distribution (e.g. Lemmon et al., 2004;Wolff et al., 2006Wolff et al., , 2009, the modified Gamma distribution (e.g. Clancy et al., 2003;Wolff and Clancy, 2003), and the Lognormal distribution (e.g. Fedorova et al., 2009Fedorova et al., , 2014. In these retrievals, in order to determine ER, EV has been fixed and its value prescribed a priori. This approach has led to consistent values for ER of approximately 1.5 μm (a canonical value), when EV is chosen to lie between 0.2 and 0.5. Moreover, it is not surprising that seasonal or regional variation of ER has also been observed, first in Clancy et al. (2003) and Wolff and Clancy (2003), and recently in Smith et al. (2016) and Vicente-Retortillo et al. (2017). ER has been reported to be correlated to atmospheric dust loading (see also Chen-Chen et al., 2019), with larger values of ER occuring mostly in northern fall and winter. Vertical variations of ER have also been observed Fedorova et al., 2009;Clancy et al., 2010;Määttänen et al., 2013;Guzewich et al., 2014), indicating a strong gravitational segregation effect in most of the year except for events of intense dust lifting, during which large particles may reach higher altitudes and the effective dust size can become more uniform (see also Montmessin et al., 2017 andKahre et al., 2017).
In many general circulation models (GCMs) of Mars, the dust cycle is usually simulated interactively: dust lifting, transportation, and sedimentation processes are parameterized by some model-resolved parameters (such as temperature, wind stress, and surface wind stress) that may be affected by the radiative effects of these simulated dust parameters.
In some numerical models, the dust cycle is simulated by a twomoment scheme (e.g. Schulz et al., 1998;Morrison and Gettelman, 2008). By assuming a constant value of EV and a prescribed dust size distribution, an ER can be obtained from the model-calculated mass mixing ratio. The most common functions that have been used to specify the dust size distribution include the Gamma distribution (e.g. Lee et al., 2018) and the lognormal distribution (e.g. Madeleine et al., 2011;Wang C et al., 2018).
In this work, we investigate the values and variations of ER and EV by using the Mars GCM MarsWRF with a two-particle scheme in which dust population is simulated by using two size bins. This kind of dust scheme is generally called an N-bin scheme, in which dust particles are represented by a finite number of tracers of different particle size. Although reasonable simulations of the global dust cycle and its thermal effects have been achieved by use of even a small number of bins (e.g., Basu et al., 2004Basu et al., , 2006 for a twobin case; Kahre et al., 2005and Neary and Daerden, 2018 for three-bin cases), so far there has been no explicit formulation for the calculation of ER and EV in this kind of scheme; fixed values for these important quantities have had to be specified a priori. But when a fixed value of ER is used in radiation calculations, the interactive feature between the radiation and the dust mixing ratio cannot be well represented. One objective of this study is to provide an analytical formulation for these two quantities that can be used in the usual N-bin schemes. Compared with the two-moment scheme, the present approach has some distinctive features. First, a priori assumptions about the dust size distribution, such as a functional form and an assumed value of EV, are not needed. Second, seasonal and spatial variations of EV can be obtained since EV is not assumed to be fixed. Third, the approach presented here can be readily extended to cases of arbitrary N.
The following study is focused on the two-particle scheme since it is the simplest case of a varying ER; and the analytical expressions derived for this simplest case can readily be extended to the Nparticle case. The two-particle scheme may indeed provide more insight into dust loading and its relationship with dust particle sizes than a more complex case. It is also worth mentioning that recent observations suggest the prevalence of a bimodal size distribution of dust (Montmessin et al., 2002Määttänen et al., 2013;Fedorova et al., 2014), challenging the traditional assumption of a continuous monomodal distribution form. These considerations les us to focus this study on the two-particle scheme, as a suitable approach for a preliminary study of this topic.
In Section 2, we derive the analytical expressions of ER and EV for the N-particle scheme, with particular attention to the case of N = 2 (two-particle scheme). Section 3 describes numerical simulations performed by the general GCM MarsWRF. The results of these simulations are discussed in Section 4. Finally, in Section 5 we summarize the study's main results and discusses their implications.

General Case (N-Particle)
In a general N-particle scheme, dust is simulated by N particle types corresponding to radius with . The corresponding mass mixing ratios of these types are traced individually and there are no interactions between different types of dust particles. At any grid point, the number density of particle type is related to the mixing ratio , where and are densities of the dust particle and the environmental CO 2 gas, respectively. All particles are assumed to be spherical in shape. Recall that for a continuous size distribution with particle radius , the kth-moment is given by the integral . In the N-particle scheme, takes discrete values and so the kth-moment can be given by a summation: where the last equality is justified by the fact that . As we can see, the kth-moment at any grid point is determined by the total number density and the tracers . It is thus straightforward to express ER and EV for a N-particle scheme, as both are given by the moments of the size distribution (Hansen and Travis, 1974): r eff n tot where and v eff are ER and EV, respectively. Notice that the factor is always cancelled out in the above expressions.

Two-Particle Case
In the two-particle scheme, using Equation (2a), Equation (2b), and taking we have In these expressions both and v eff can be fully determined by and . Unlike the approach in which a continuous size distribution is implemented in a two-moment scheme (v eff has to be fixed for solving from tracers), both and v eff in the above twoparticle scheme can be solved independently and hence may provide a way to evaluate the spatial and temporal variation of v eff .
The functional properties of and v eff can easily be seen by noti- cing that both Equation (3a) and Equation (3b) can be converted to single variable functions. Let us define a dimensionless variable that can be interpreted as the relative abundance among particle types at a particular grid point. For a fixed total mass mixing ratio of dust , dust particles of the finer mode ( ) may be more responsible for a large value of R while dust particles of the coarser mode ( ) will be relatively unimportant. In the extreme cases, when and , may approach the limit of +1 and -1 respectively. In fact, and can also be converted to single variable functions by using the variable instead of . However, the range of this variable will then become 0 to . Now, in terms of , we have We can see that is a monotonically decreasing function of , while v eff behaves as a parabola centered at . The minimum and maximum values of are respectively and , corresponding to and respectively. This implies that is close to the radius of the dominant mode at the grid point when either or is negligible there. Also, when dust particles are generally larger, will likely take a larger value. In this case, the minimum and maximum values of v eff are respectively 0 and , corresponding to and 0 respectively.

Numerical Model and Simulations
To illustrate the quality of our approach, we utilize a robust numerical model MarsWRF to generate global circulation of the Martian atmosphere, especially its dust cycle. MarsWRF is the Marsdedicated version of PlanetWRF (Richardson et al., 2007;Guo X et al., 2009;Toigo et al., 2012), and was developed on the basis of the Weather Research and Forecasting (WRF) model for terrestrial weather and climate studies. MarsWRF is a grid point model utilizing Arakawa C-grid in horizontal directions, with 36 and 72 grids in latitudinal and longitudinal directions respectively, corresponding to the resolution of 5° squared or about 300 km squared in the equatorial region. The model has 52 terrain-following hydrostatic pressure layers, defined by surface pressure and a fixed model top pressure (taken as 0.0057 Pa). The basic model configuration and the physics schemes used are effectively the same as those used in Chow et al. (2018) and Xiao J et al. (2019).

ΔT
As mentioned in the previous section, dust is simulated numerically in this work by two size bins with particle radii of 1 and 3 μm. Processes of dust lifting, transportation, and sedimentation for each bin are controlled by model-resolved conditions. Dust lifting is parameterized by two schemes used in Newman and Richardson (2015) with some slight modifications (see below), and dust is assumed to be available everywhere and at all times over the whole planet surface except those surfaces with ice cover. In the first scheme the lifting of dust is a function of surface wind stress. Dust lifting occurs over the surface when the local near-surface stress exceeds a particular threshold value (a constant value 0.043 N•m −2 ). The second scheme parameterizes dust lifting due to thermal convection similar to dust devils. Dust lifting occurs when the temperature difference between surface and the ΔT/T surf surface air exceeds a certain threshold value (27 K in this case). The amount of dust lifting is determined by thermodynamic efficiency , as well as by the sensible heat flux. Dust lifted to the atmosphere is transported by the model-resolved wind, and then settles down under gravity according to the size-dependent Stokes-Cunningham relation.
The short-wave and long-wave radiations are evaluated by the Wide Band Model (WBM) as described in Richardson et al. (2007), which considers the radiative interaction of dust with the CO 2 atmosphere; dust in the atmosphere may change the atmospheric radiation and thus the circulation. The contribution of dust to long-wave and short-wave radiation is evaluated following Haberle et al. (1982) and Briegleb (1992), respectively. The change in the shortwave extinction opacity due to the suspended dust is given by (Madeleine et al., 2011) where is the gravitational constant, is the density of dust particles, and is the effective radius given by Equation (3a). The extinction coefficient is considered to be sizedependent. Its values are evaluated by the Python package miepython, which solves the Mie scattering theory numerically when a complex index of reflection is given. The values are then smoothed by a very narrow Gamma distribution of (Hansen and Travis, 1974, Fig. 8). In this work, the observational data from Wolff and Clancy (2003) have led us to adopt and for the particle sizes of 1 μm and 3 μm respectively. For simplicity, other radiative parameters (e.g. single scattering albedo and asymmetric factor) are assumed to be size-independent. This setting is reasonable as WBM considers only averaged effects over a wide range of wavelengths, and so is basically insensitive to differences in dust size. With a more sophisticated radiation scheme (e.g. correlated-k scheme), a full set of size-dependent parameters will be needed.
In this study, our primary goal is to explore the values and variations of ER and EV with Equations (3a) and (3b). To evaluate the implementation of this approach in a Mars GCM, we perform a MarsWRF simulation (referred as SimMain) with the aforementioned setting for 16 Martian Years (MYs); our main discussions will focus on results from this simulation. Moreover, we investigate whether the variability of ER is important in a Mars GCM simulation. For this purpose, we perform another 16-MY simulation (referred to as SimRef) as a reference, which has the same configuration as SimMain except for its use of the fixed values of 1.5 μm and for evaluating the dust opacity (Equation (5)). A comparison of SimMain and SimRef should reveal more information about the effect of including variable ER in a Mars GCM simulation.

Results of Simulations
Ref results are discussed and compared with those of SimMain.

Dust Climate
In the 16 years of the SimMain simulation, the dust climates in 13 years show a regular pattern (Figure 1a) in which the simulated zonal-mean and column-integrated dust extinction optical depth (CDOD) are consistent with some derived observations (e.g., Madeleine et al., 2011;Montabone et al., 2015). For example, the dichotomy of the dust loading seasons can be capturedlow-dust-loading (LDL) season in the northern spring and summer and high-dust-loading (HDL) season in the northern fall and winter. The simulation also captures other prominent features of the dust climate including the two-episode feature in the equatorial region and the northern mid-latitude region (Xiao J et al., 2019), and the corresponding solsticial pause around solar longitude L s = 270° (Lee et al., 2018). A regular year climate has been obtained by averaging the results of these 13 years.

Dust Size
In Regular years, values and variations of zonal-mean and columnmean (defined by for variable , where is surface pressure and the integration is from the surface to the top of the atmosphere) ER and EV, can be seen in Figure 1b and 1c respectively. In the beginning of the year, ER can be as little as 1.1 μm (near polar regions), increasing gradually to 1.3 μm around L s = 60° (Figure 1b). Its value rises to 1.5 μm around L s = 100°, and reaches its annual peak of 1.75 μm around L s = 240° in the equatorial region. This range of ER is consistent with the observed annual average of 1.5 μm (e.g., Pollack et al., 1995;Tomasko et al., 1999;Clancy et al., 2003;Wolff and Clancy, 2003;Lemmon et al., 2004;Wolff et al., 2006Wolff et al., , 2009. Moreover, spatial and temporal variations of ER are at levels consistent with those reported by, for example, Wolff and Clancy (2003), Vicente-Retortillo et al. (2017), and Chen-Chen et al. (2019), who report that ER could be as small as 0.6 μm and as large as 2 μm. In the beginning of the year, ER is about 1.1 μm and it reaches to about 1.75 μm in the dusty season (northern fall and winter), indicating that the contribution of large particles is mainly from a wind stress lifting process. Since smaller particles can suspend in air for a longer duration, one could expect that even smaller values of ER could be reached if more bins corresponding to smaller particles are introduced. Also, we can see an approximate correlation between dust size and dust loading. There is another kind of dichotomy found in the seasonal variation of ER, in which large ER usually happens in HDL seasons. In contrast to the CDOD, ER presents only one peak at the equatorial region, while a two-peak pattern can still be found in the northern mid-latitudes.   weaker. For most of the year, EV is found to be greater than 0.2 (except in regions near the poles) and reaches a still greater value (about 0.29) from L s = 90° to 360°. It reaches its peak, around 0.32, from L s = 260° to 310°. In general, EV is greater in this period. Unlike ER, for EV there is neither an apparent dichotomous pattern nor a good correlation to the CDOD episodes. This finding is indeed consistent with the assumption of some observations (e.g., Clancy et al., 2003;Wolff and Clancy, 2003), in which dust size is retrieved by prescribing a single-valued EV for the underlying size distribution. Also, a relatively uniform value of EV is qualitatively consistent with the previous model study (Kahre et al., 2008) in which the EV seems to be uniform. From Figure 1c, a very small EV (as small as 0.1) can be found in the polar regions in the northern spring. This could be, again, due to the limited number of bins presented in our model. When dust loading is dominated by just one bin at a model grid point, variance of particle size is in principle close to zero. This situation happens easily in northern spring since wind stress lifting is weak in this period (Kahre et al., 2006). With more bins filled with smaller particle sizes, EV could increase and become smoother since the smaller particles have significant suspension ability.
Besides column-averaged quantities, vertical profiles may tell more about the spatial distribution of ER and EV. We consider two periods of L s = 90° and L s = 255°, which respectively represent the clearest and most dusty time in the Regular year climate. The zonal-mean vertical profiles of ER and EV at these times are shown in Figure 2. Again, we observe that ER is generally smaller in the LDL season ( Figure 2a) and greater in the HDL season (Figure 2b). In both periods, ER is larger near the ground and decreases as altitude increases. Large particle can extend to a greater height in low latitude regions compared with that in high latitude regions. These patterns show a strong gravitational segregation, and are qualitatively consistent with the findings of Madeleine et al. (2011) in which dust size is simulated by a two-moment scheme.
For EV, the general patterns in the vertical profiles are similar to those for ER. However, although values of EV in the LDL season ( Figure 2c) are generally smaller than those in the HDL season (Figure 2d), the difference is not as apparent as that of ER. This result is indeed consistent with that in Figure 1c. In both periods, greater values of EV can extend to a greater height in the low latit-0.32 ude region. At L s = 90°, the EV maximum occurs near the ground around the latitudes of 10°N to 50°N. However, at L s = 255°, the EV maximum ( ) occurs at the height of 25-30 km at the equatorial region. The altitudes of maximum EV decrease northward and southward from the equatorial region to form a bell shape (see Figure 2d).

Effect on Radiation Process
It is interesting to ask what effect the variability of ER will bring about in a Mars GCM. The effect can be inspected by considering differences between the simulated temperature fields of SimMain and SimRef. The present discussion is based on results from the Regular year climate in SimMain and in SimRef.
From the vertical profiles of the zonal-mean dust mixing ratio, averaged over the latitudes from 40°S to 40°N (Figure 3a), we can see that the results at L s = 90° are almost identical for SimMain and SimRef. On the other hand, for the corresponding values of ER ( Figure 3b) and EV (Figure 3c), SimRef values are slightly larger than SimMain values. At L s = 255°, the situation is different. ER and EV are almost identical between the two simulations while the dust mixing ratio in SimMain (Figure 3a) is apparently larger than that in SimRef.
Recall that dust opacity is proportional to the factor , where is given by Equation (3a) in SimMain, and is taken as 1.5 μm in SimRef. As height increases, the decreasing ER in SimMain pushes this factor up, thus generating the difference between SimMain and SimRef. Moreover, the difference in dust mixing ratios between the two simulations would also contribute to . Figure 3d shows the vertical profile of at the two chosen time periods, calculated directly with the zonal-and latitudinal-mean values shown in Figures 3a-3c. From Figure 3d we can see that reaches its maximum at the heights of about 20 and 30 km at L s = 90° and 255° respectively. Moreover, the maximal values at L s = 255° are larger than at L s = 90°, mainly due to the difference in dust mixing ratio of the two simulations at that time.
q tot /r eff A difference in implies a difference in opacity, thus in radiative heating rate, and thus in the temperature field. The vertical profiles of the zonal-mean temperature differences between Sim- Similarly, at L s = 255°, even more apparent positive temperature anomalies occurs at heights from about 20 km to at least 60 km, with the maximum at between 35 and 40 km for the same latitudinal region. These are warming effects generated by differences in indicated in Figure 3d. The heights of maximums of are consistent with the heights of maximums of temperature anomalies at both times. This pattern demonstrates that including dust size variation can be expected to affect the calculated thermal structure, revealing that the use of fixed or variable ER could significantly alter the modeled thermal structure.

Discussion and Conclusions
Effective radius ER and effective variance EV are two important quantities concerning the distribution of dust particle sizes, and have been extensively studied as reported in many previous studies. In most Mars GCMs, these quantities are commonly evaluated either by a two-moment scheme in which ER is related to the particle concentration through a prescribed distribution function, or by an N-bin dust scheme in which dust particles of N sizes are assumed in modeling dust processes. The value of ER is generally chosen to be fixed in calculating the radiation associated with dust in N-bin schemes, although ER is known to be variable in reality. In this study, an approach to evaluating ER and EV has been proposed when an N-bin dust scheme is used in a Mars GCM. In this approach, ER and EV are given by analytical functions of the model-resolved dust mass mixing ratios. The functions are derived from first principles and details of the underlying size distribution and parameters are not required. In this approach, the effect of ER variation is taken into account during the calculation of radiation.
To illustrate the application of these formulations, we have performed numerical simulations with a GCM that considers dust particles of two sizes (a two-particle scheme) and the dust particles are interactive with the radiation. Based on the present formulation, we have obtained values and variations of ER that are consistent with previous observations and modeling studies.
Based on results of the present approach's simulations, the variation of EV with different dust assumptions has been evaluated. The values of EV are usually assumed to be constant in many   GCMs. We find that the seasonal variation of EV exhibits a pattern generally similar to that of ER, but is relatively less apparent. In most of the year, EV is not significantly changed at some observed locations on Mars. This uniformity of EV is usually assumed in most of the literature, and this property has been further supported by the present study.
To investigate the significance of variable ER in numerical simulations, we have compared two simulations. One adopts a variable ER (from the present formulation) for the opacity calculation while another adopts a fixed valued of ER. The results suggest that, due to the vertical variation of ER (and increase of dust loading), an additional heating is present at heights between 10 and 30 km in the LDL season of L s = 90° (heights between 20 and 50 km in the HDL season at L s = 255°) in the case with a variable ER. This effect is similar to that reported in the literature.
The present study is a preliminary investigation of the problem of dust size distribution and its interactive effects on the dust radiative process. In this preliminary work, we have treated as size-independent some parameters in the radiative process (such as the scattering albedo, asymmetric factor, and extinction coefficient) that could be dust-size dependent. A more sophisticated investigation in the future should consider the size-dependent nature of these parameters.
Finally, we have used only a two-particle scheme to illustrate the formulations derived in this study, and have not extended the scheme to more dust bins. Nevertheless, we believe that the simple two-particle scheme is illustrative and should not significantly degrade the main results of the present study. The two-bin model is the simplest scheme possible to implement the present formation, and the simplest scheme to reproduce the observed Martian dust climate. We do not intend to suggest that such a simple scheme is a realistic model. We hope that future studies will apply the present formulations to a GCM with more dust bins to obtain more realistic results.