Modelling coarse and giant desert dust particles

. Dust particles larger than 20 µm in diameter (0.2 μm < D < 100 µm) have been regularly observed to remain airborne during long-range transport. In this work we extend the parameterization of mineral dust cycle in the GOCART-AFWA dust 15 scheme of WRFV4.2.1, to include also such coarse and giant particles. The initial particle size distribution in our parameterization is based on observations over desert dust sources and the Stokes’ drag coefficient has also been updated to account for dust particles of all sizes (Re < 10 5 ). The new code is applied to simulate dust transport over Cape Verde during the August 2015 AER -D campaign. Model results are evaluated using both airborne dust measurements and the CALIPSO-LIVAS pure dust product. The results show that the modeled lifetimes of the coarser particles are shorter than those observed. 20 Various processes are proposed to explain such inaccuracies, such as the electric field inside dust plumes and non-spherical aerodynamics. Additional sensitivity runs are performed by artificially reducing the settling velocities of the particles to compensate for such underrepresented processes in the model. Our simulations show that particles with diameters of 5-17 μm and 40-100 μm are better represented assuming 80% reduction in settling velocity (UR80) while particles at the range 17-40 μm are better represented in the UR60 scenario. The overall statistical analysis shows that the UR80 experiment presents the 25 closest agreement with the airborne in situ measurements both in Cape Verde and over the sources. The experiment improves also the vertical distribution of dust in the model, as compared to the CALIPSO-LIVAS pure dust product. Further research is requested in order to understand the physical processes behind the reduction of settling velocity.

Where is the particle diameter in (μm) and is the particle number size distribution in number of particles per cm -3 . The parameters of each size bin are shown in Table 3. Henceforward, references about the size of the particle correspond to particle volume equivalent diameter, unless mentioned otherwise. 130 In the default GOCART-AFWA of WRF, the total emitted vertical dust flux is estimated at each grid point prone to dust emission, when favorable conditions are met. The dust flux is then distributed over five transport size bins, based on the fragmentation theory of Kok, (2011), although limited to diameters up to 20 μm. Since our goal is to include larger dust particles than those commonly used in current atmospheric dust models, we redefine the five transport model bins to include particles with diameters up to 100 μm (Table 1). We use a prescribed PSD for emitted dust particles at the source based on in 135 situ measurements from the FENNEC campaign (Ryder et al., 2013a). Ryder et al., (2013a) made airborne in situ measurements of dust PSDs at various altitudes near dust sources in the Sahara Desert. The emitted dust PSD used in our work is derived from measurements of fresh upwelling cases at the lowest available altitudes from aircraft profiles representative of 1 km and is hereafter referred to as the "observed FENNEC-PSD". The observed FENNEC-PSD is shown in Fig. 2(a) with red squares, and the shaded areas show the size range of the individual bins. In Sect. 2.2.1 more information are provided about the 140 FENNEC campaign and the instruments used for the measurements.
The distribution of emitted mass over the redefined size range is obtained by calculating the mass fraction resulting from the weighting factors ( ) for each transport bin, as shown in Eq. 2. The weighting factors are also shown in Fig.2 145 Where is the particle diameter, dlnD is the volume size distribution in μm 3 cm -3 , , and , are the margins of each size bin in μm.

Updated gravitational scheme
In WRF-GOCART-AFWA, the forces acting on a dust particle moving along the vertical direction, are the gravitational force 150 and the aerodynamic drag force , which are mathematically expressed in Eq.3 and Eq.4, respectively. 155 The constant velocity that a particle builds up, as it falls vertically in the Earth's atmosphere, is defined as the terminal settling velocity , and it can be estimated by solving the 1-D equation of motion in the steady state limit, where ΣF is assumed to be equal to zero: https://doi.org/10.5194/acp-2022-94 Preprint. Discussion started: 7 April 2022 c Author(s) 2022. CC BY 4.0 License. 160 Where is the particle density in 3 , is the gravitational acceleration in 2 , is the particle volume in 3 and is the particle projected area normal to the flow in 2 , is the atmospheric air density in 3 and is the aerodynamic drag coefficient (unit less). For each size bin it is assumed that the particles are spherical with diameter in (as defined in Sect. 2.1.1), thus their volume and projected area are defined by the following equations for spheres: 165 The drag coefficient is that of Stokes' Law and is defined as: 170 Where is the Reynold's number (unit less) given by the following equation: Where is the air dynamic viscosity in • defined as a function of air temperature in o by the following equation: Equation 7 has been derived with the simplification of no-slip boundary conditions, thus a correction factor , proposed by Davies, C. N. (1945), is applied to the Stokes' relationship to account for velocity slip at the particle's surface. The corrected drag coefficient become is: 190 Where = 1.1•10 −3 •√ is the air mean free path in and is the air dynamic viscosity in • , as defined by Eq.10. is the air temperature in o and the air pressure in ℎ .
Substituting Eq. 6-9 in Eq. 4 we end up with the relationship for the terminal velocity of the dust particles, as shown in Eq. 12.
Using Eq.5, 6, 7, 9 and 14 we calculate the terminal velocity for each model size bin. Since the resulting equation is not linearly 205 dependent by we apply the bisection method to solve the equation.
In the default code the slip correction is applied unconditionally, as mentioned above, for all the values of . However, slip correction is defined in Stokes' regime (Mallios et al., 2020). Thus, in the updated drag coefficient, only when < 0.1 (Stokes' regime), we recalculated the settling velocity using the corrected drag coefficient . = ′ / ′ , where ′ = ( ′ ) with ′ the mean free path adopted by (Jennings, 1988): 210 here μ is air dynamic viscosity in • , as defined by Eq.10, and the atmospheric pressure is in Pa.

Model experiments 215
Using the WRF-L code, we run a simulation that serves as a CONTROL experiment. Our simulation period coincides with the AER-D experimental campaign and covers the days from July 29, 2015 to August 25, 2015 for a region extending in latitude  (Fig. 3). The simulation area is located over the major Saharan sources and also includes the downwind areas in the eastern sector of the tropical Atlantic. We use an equal-distance grid with a spatial grid spacing of 15 km x 15 km that includes 550 × 300 points and 70 vertical sigma pressure levels up to 50 hPa. For 220 each run, 84-hour forecast cycles are performed and reinitialized every 3 days using the 6-hour Global Forecast System Final Analysis (GFS -FNL) reanalysis product, available at a 0.25 o x0.25 o model grid, to initialize the model and set boundary conditions. The first week of the simulation served as a spin-up run for the accumulation of the background dust loading and is excluded from the analysis. The simulation run in a dust-only mode, without including radiative feedback from aerosols, to avoid (in this first case) a more complicated analysis that would include the radiative effect on dust transport. The scaling of 225 the dust source strength is chosen to best match the modeled DOD with the AERONET measurements (RMSE=0.34, bias=-0.07) from the desert stations: Banizoumbou, Dakar, El_Farafra, Medenine-IRA, Oujda , Tizi_Ouzou, Tunis_Carthage ,Ben_Salem). We only use the measurements where DOD is higher than 0.75 and the Angstrom exponent is lower than 0.2 to ensure that contamination by aerosols other than dust is negligible. The complete configuration options for the run are listed in Table 2. 230 In addition, we investigate the implications of a possible mechanism to counteract gravitational settling in order to reduce the differences between the CONTROL run calculations and the in situ observations (shown in Sect. 3.4). To this end, we perform additional sensitivity tests by reducing the settling velocity by 20 to 80%, with a step size of 20%. The experiments are referred to as the "URx experiment", using the percentage (x%) by which the settling velocity is reduced. With this artificial tuning, we aim to reproduce the net force acting on dust particles falling into the atmosphere and overcome the current shortcomings 235 of the model (i.e., the absence of all real forces that determine the lifetime of dust particles in nature).
It should be noted that several studies have pointed out the importance of fine-resolution dust simulations ( Solomos et al., 2012;Basart et al., 2016;Gu et al., 2021;), which, among other things, help the model resolve small-scale dynamics and account for possible interactions between different scales. Given the complicated meteorological conditions during the study period (i.e., August 2015), the fine resolution increases the accuracy of the dust simulations and provides a good estimate of 240 the magnitude of the missing mechanism. The reduced deposition of particles can be attributed to either an updraft counteracting gravity or a reduction in particle settling velocity, both of which slow dust deposition rates. In the first case, this can be attributed to either as yet unresolved meteorological conditions (e.g., small-scale haboobs, dunes) or atmospheric feedbacks due to dust-radiation interactions (i.e., atmospheric heating due to absorption of solar radiation by mineral particles).
where , , , and , are the dust mass concentration in g/m 3 , the particle density in g/cm 3 , and the effective diameter in μm of size bin . 550, is the extinction efficiency at 550 nm, calculated using the Mie scattering code (Mie, 1908), 260 considering a spherical shape for the dust particles, and a refractive index of 1.55 + i0.005, which is representative of dust (e.g. Dubovik et al., 2002). For simplification of the computations, we assume that the particles in each size bin have the same size  300μm, while during the AER-D campaign, the PSD was provided for particles with diameters up to 200μm. Full details of instrumental measurements and processing are given by Ryder et al., (2013b) and Ryder et al., (2018), for FENNEC and AER-D, respectively. In Sect. 2.1.1 we describe the way that FENNEC campaign measurement used in this study. The LIVAS pure-dust product has been used to a variety of dust-oriented studies including the investigation of the dust sources and the seasonal transition of the dust transport pathways Proestakis et al., 2018); the evaluation of the performance of atmospheric and dust transport models (e.g. Tsikerdekis et al., 2017;Solomos et al. 2017;Georgoulias et al., 2018;Konsta et al., 2018), the evaluation of new satellite-based products (e.g. Georgoulias et al., 2016;Chimot at al. 2017;300 Georgoulias et al., 2020;Gkikas et al., 2021), and on dust assimilation experiments (Escribano et al., 2021). Herein, the LIVAS pure-dust extinction product is used for the assessment of the simulated dust vertical patterns. In the geographical region of our study, the uncertainty of the product is estimated to be less than 20% in altitudes up to 6km ). Sahara and eastern Tropical Atlantic Ocean), MIDAS DODs covariate (R ~0.90) very well with AERONET-derived DODs, although they are slightly overestimated by <0.04 (see Fig. 4 in Gkikas et al., (2021)). Moreover, the intercomparison of MIDAS, LIVAS and MERRA-2 DODs show a remarkable consistency in reproducing the seasonal cycle of dust loads over the W. Sahara and the eastern segment of the Tropical Atlantic Ocean. Overall, the MIDAS dataset is quite useful for the 315 current study, due to the high reliability of the derived DOD product and the product availability at fine spatial resolution, on a daily basis.

MSG-SEVIRI-DUST RGB product
We use the Meteosat Second Generation -Spinning Enhanced Infrared and Visible Imager (MSG-SEVIRI) DUST RGB product, which is produced by the RGB colors (Red-Green-Blue), corresponding to the three infrared channels of the MSG-320 SEVIRI instrument. The functionality of the geostationary SEVIRI sensor in the infrared area of the electromagnetic spectrum, and the combination of the different sensitivities of the three channels, enables both daytime and nighttime continuous observations, along with the discrimination between land, clouds and aerosols, making the Dust RGB product very useful for monitoring intense dust and volcanic ash plumes. Dust particles are depicted on images as bright magenta (during day) or purple color (during night) over land, and as a magenta color over the sea. 325 Figure 4 shows the altitude profiles of the settling velocities for each size bin from the CONTROL run, averaged over the simulation domain, and the simulation period. As the size gets bigger the settling velocity is increased. Terminal velocities of particles of bin 5 are two orders of magnitude greater than the particles of bin 2 and bin 3, and one order of magnitude greater 330 than the particles of bin 4. The terminal velocities increase with height following the temperature lapse rate, and are sensitive to the thermodynamic condition of the atmospheric air, increasing as temperature or air density drops, based on Eq.10, 13 and the relationships of air viscosity. The average settling velocities for the CONTROL run near the surface differ from those at 6 km height, by approximately 10 %, which is a significant reduction, especially for coarser and giant particles where velocities are greater. 335 A direct comparison between the modelled and the observed PSD for the dust concentration above the sources, cannot be conducted for the AER-D campaign, since the measurements were only performed over the ocean. Figure 5 shows a more 345 qualitative comparison, using the observed FENNEC-PSD at 1km (red squares). The modelled and observed PSD differ. The modelled volume concentrations have larger values for bins 1-3, and lower for bins 4 and 5. The maximum concentration of the modelled PSD is at bin 3, whereas for the observed PSD is at bin 4, suggesting that the model underestimates the concentrations at bins 4 and 5, already from the initial transport stage, near the dust sources. Those differences can be attributed either to an underestimation of the contribution of the coarser particles on the emission, to an overestimation of their loss 350 during transport from the surface to 1 km, or to higher updrafts that remain unresolved in the simulation of this study.

Mean dust load and spatiotemporal distribution of dust
In order to further demonstrate the distribution of the total dust mass at the different sizes, Fig. 6 shows the simulated fields of the total columnar dust load, along with the corresponding concentrations at each size bin. The simulations in Fig. 6 are performed using the parameters of the CONTROL experiment, and the calculated concentrations are averaged over the period 355 of interest (5-25/8/15). For the first three bins, the spatial patterns of dust load are very similar, showing the dust sources in the Western Sahara, and the advection of the particles towards the Atlantic Ocean. The mass increases from bin1 to bin 3 (5.5 -17 μm), which has maximum values for the whole size range. Dust particles with diameters between 17 μm and 40 μm (bin 4) are found mainly over land, and are subjected to short-range transport westwards (i.e., off the Moroccan coast). Giant particles (bin 5) are found at very low concentrations (< 0.5 gr m -2 ), at isolated areas over/near dust sources, probably due to 360 their quick gravitational settling.
The comparison of the model simulations with satellite retrievals shows that, in general, there is a good agreement on the spatiotemporal distribution of dust during the days and times of the AER-D flights. Deviations between the simulations and the observations are found for flight b920, due to a shift of the center of the simulated dust mass towards the south (Fig. 7(a)).

Moreover, the observations show that the dust plume traveled towards Morocco and Canary
Islands, whereas the model shows 365 that it traveled mainly towards Cape Verde (see Dust RGB image of MSG-SEVIRI, during the time-of-flight b920 in Fig.7(b), and the MODIS DOD and corrected reflectance in Fig. 7(c) and 7(d), respectively). This difference results in an overestimation of the simulated dust mass in the area of South West Africa and West Mauritania, affecting the transport towards the area of flight b920. The main cause for this discrepancy is the difference in the intensity by which the various dust sources in Northern Algeria were activated during the previous days. As it is depicted in Fig. 7(e, f) there are sources in the model that have been 370 strongly activated in circles A and B, although in RGB-Dust MSG-SEVIRI images they are depicted with much less intensity (fewer pink colors). That deficiency of the model could be attributed to various reasons, such as underrepresentation of the meteorological conditions mostly in cases of haboobs (Ryder, 2018a(Ryder, , 2021 Fig. 5(a). Considering the different experiments for reduced settling velocities (Table 3), we see that the 385 reduction mainly affects the simulations for the coarser particles (bins 3, 4 and 5), with the effect increasing with the size of the particles. The simulated concentrations of giant particles (bin 5) become significant when the reduction in settling velocity is greater than 60%. The comparison of the observed and modelled mean average PSDs in Fig. 9(a) shows that UR60 and UR80 experiments are closer to the observations, with UR80 to better reproduce the observed values of bins 3 and 5, whereas UR60 better reproduces the values of bin 4. This reduction results in settling velocities of ~0.066 m/s for bin 3 (D=5.5-17 μm), 390 ~0.32 m/s for bin 4 (D=17-40 μm) and ~1,88 m/s for bin5 (40-80μm). In general, UR80 simulations of the mean PSD provide the best agreement with the observations. In terms of total volume, the UR80 simulations have the smallest relative difference with the observations for most flights, providing a ~50% improvement in relative difference, as it is depicted in Fig. 9(b).
UR80 also provides better agreement with the observed FENNEC-PSD above the dust sources, by shifting the maximum of the PSD to bin 4 (Fig. 5). 395

Dust vertical distribution
Figure 10(a) shows the profile of the mean extinction coefficient at 532 nm, provided by the LIVAS pure-dust product (black line), and the profile of the mean extinction coefficients at 550 nm, provided by the CONTROL, UR20, UR40, UR60, and UR80 experiments. This comparison is an initial demonstration of the good initialization and performance of the different model experiments, with respect to capturing the vertical distribution of dust. The intercompared profiles are in good 400 agreement, with the simulations falling well-within the variability of the dust observations, although discrepancies are also present, especially close to the dust sources, in the nighttime boundary layer (Fig.10(b) region I), and within the upper free Troposphere (Fig. 10(b) region III). The discrepancies close to the dust sources are attributed to the complex topography, in terms of geographical characteristics (Proestakis et al., 2018), weighting effects, surface returns, and representativeness issues related to the aggregation of CALIPSO L2 profiles to LIVAS 1 o x1 o grid resolution (Amiridis et al, 2013. The discrepancies in the upper free Troposphere (above 6 km) are attributed to the presence of tenuous aerosol layers which fall below the CALIOP layer detection threshold. Thus, the assessment of the different model experiments with the LIVAS pure-dust product, is performed in the region between 1.5 km and 6.4 km a.m.s.l. (Fig. 10region II).
According to the comparison of observations and simulations of the mean extinction coefficient (Fig. 10(a)), the statistical overall analysis reveals that the UR40 experiment demonstrates a better performance compared to LIVAS, reducing the mean 410 bias close to zero. However, the UR80 experiment provides a more constant (positive) bias with height, which suggests a better distribution of the dust mass in the vertical.

Discussion and Conclusions
The frequent presence of large desert dust particles (D>20 μm) far from their sources, is well established by numerous observational studies over the last decade (van der Does et al., 2018;Liu et al., 2018;Ryder et al., 2013bRyder et al., , 2018aRyder et al., , 2019a415 Weinzierl et al., 2009415 Weinzierl et al., , 2011415 Weinzierl et al., , 2017b. However, the processes that result in the particle retainment in the atmosphere, and subsequently their travel at greater distances than predicted, remains unrevealed. In this study we extend the particle size range acknowledged in WRF-GOCART-AFWA transport code, to include particles with diameters up to 100 μm. The evaluation against airborne in situ observations of the size distribution shows that larger particles, are underestimated, both above their sources and far from them. This suggests that there are atmospheric processes that are not taken into account in the model 420 simulations. We investigate the effect of reducing the settling velocity of the dust particles due to these unknown processes, and we see that for a reduction of 60% and (especially for) 80%, the simulations of the PSD in Cape Verde are improved with respect to the observations. The reduction of 80% results in settling velocity of 0.066 m/s for particles with D<25 μm, which is double than the value reported by Maring et al. (2003) for similar sizes. It should be noted though that Maring et al. (2003) derived this settling velocity using observations that were taken with a five-year difference. Ginoux (2003), has also reported 425 an improvement in model simulations for a reduction in settling velocity of approximately 45% and 60%, for particles with diameters 10 to 30 μm. Though, the differences in the model resolution, the dust scheme and the drag coefficient in Ginoux (2003) compare to this study, could cause the different values of the required corrections in the settling velocities. The difference with the values suggested herein, can mainly be attributed to the different drag coefficient used in Ginoux (2003), which results in lower settling velocities for the spherical particles. 430 One of the processes proposed in the literature to explain the longer atmospheric lifetimes of large mineral dust particles is the particle asphericity. Huang et al. (2020) used globally averaged shape distributions of particle aspect ratio and height to width ratio and provided a correction to the spherical particle settling velocity, which is valid for ellipsoidal particles. According to their empirical expression, there is a 20% reduction of particle settling velocity in the case of ellipsoidal particles compared with their spherical counterparts of the same volume. Among the limitations of their methodology (see Huang et al. 2020), is 435 that it is valid only in the Stokes' regime (Re<<1), which limits the applicability of the study for particles with sizes less than 20 μm, and that the ellipsoidal particles are randomly oriented. in the Earth's atmosphere in the presence of electrical and gravitational forces. They found that the random orientation assumption is, in principle, valid only for particles with size (two times the particle major semi-axis) less than 2 μm, reducing 440 even more the applicability of the methodology by Huang et al. (2020). As the size increases, the gravity or the electrical force tend to create sufficient torque to rotate the particle horizontally or vertically with respect to the ground, respectively (depending on the particle net electrical charge and the large-scale atmospheric electric field).
Moreover, Mallios et al. (2020) provided new expressions for the drag coefficient of prolate spheroids that are valid beyond the Stokes' regime (specifically for Re≤100) and that take into account the orientation and the aspect ratio of the particle. They 445 showed that in the case the aspect ratio ranges between 1.4 and 2.4, prolate spheroids fall faster than their spherical counterparts of the same volume. This is attributed to the projected area of the prolate spheroids, which depends strongly on the particle orientation (although on average it is larger for ellipsoids than spheres (Vickers, 1996), the projected area of ellipsoids becomes smaller than the projected area of spheres of the same volume as the particle becomes vertically oriented), and the aerodynamic properties due to the impact of the prolate spheroid shape factors on their drag coefficients. They also showed that when 450 comparing prolate spheroids with spherical particles of the same maximum dimension the conclusions are different. In the case of particles with aspect ratio equal to 1.4, the settling velocity of prolate spheroids is on average 6% (in the case of horizontal orientation) up to 23% (in the case of vertical orientation) less than their spherical counterparts (of the same maximum dimension). As the aspect ratio increases to 2.4, the difference becomes 20% (for horizontal orientation) and 52% (for vertical orientation). 455 Another process that can influence mineral dust settling has to do with the electrical properties of dust particles. The dust particles are charged in the atmosphere either due to the attachment of atmospheric ions on them (Mallios et al. 2021b) or/and due to collisions, a process known as triboelectric effect (Ette, 1971, Eden and Vonnegut, 1973, Mills, 1977, Jayaratne, 1991. Moreover, there is a large-scale atmospheric electric field, due to the potential difference between the lower part of the Ionosphere and the Earth's surface (Rycroft et al., 2008). The electric field is modified by ion attachment process (Mallios et 460 al. 2021b) or by the charge separation caused by updrafts (Krauss et al., 2003). Therefore, electrical forces are generated that might influence the particle settling process by balancing the gravity or changing the particle orientation. The quantification of the particles electrical properties is still an open question.
Triboelectric effect is able to modify the particles saltation process at the emission sources right above the ground due to the generation of very high values of electric charge caused by the large collision frequency which is a consequence of the wind 465 and the large particle number density (Kok andRenno, 2006, 2008). As the particles are aloft, the collision frequency decreases (Rahman et al., 2008) and the ion attachment process can modify the acquired particle charge, because the electric field of the charged particles tend to attract ions of opposite polarities (Mallios et al., 2021b). It is still unknown if the acquired charge of the particles remains or is neutralized. Toth III et al., (2020) estimated that if the particle net charge persists, then the ambient electric field is sufficient to generate electrical forces that can keep particles suspended at higher elevations and enrich the 470 concentration of larger particles at higher elevations. Mallios et al. (2021)