Modeling of Spherical Dust Particle Charging due to Ion Attachment

The attachment of positive and negative ions to settling spherical dust particles is studied. A novel 1D numerical model has been developed to parameterize the charging process in the presence of a large-scale electric field. The model is able to self-consistently calculate the modification of atmospheric ion densities in the presence of the dust particles, and the consequent alteration of the atmospheric electrical conductivity and the large-scale electric field. Moreover, the model estimates the acquired electrical charge on the dust particles and calculates the electrical force that is applied on them. Using observed dust size distributions, we find that the particles can acquire electrical charge in the range of 1–1,000 elementary charges depending on their size and number density. The particles become mainly negatively charged, but under specific conditions giant mode particles (larger than 50 μm radius) can be positive. Moreover, the large-scale electric field can increase up to 20 times as much as the fair weather value. However, our approach shows that the resultant electrical force is not enough to significantly influence their gravitational settling, as the ratio between the electrical force magnitude and the gravity magnitude does not exceed the value of 0.01. This indicates that the process of ion attachment alone is not sufficient to create strong electrical effects for the modification of particle dynamics. Therefore, other processes, such as the triboelectric effect and updrafts, must be included in the model to fully represent the impact of electricity on particle dynamics.


INTRODUCTION
Mineral dust plays an important role in the Earth's atmosphere, and in the Earth System in total. It significantly impacts radiation (Li et al., 2004), as dust particles absorb shortwave and long-wave radiations, potentially causing a net atmospheric warming (Kok et al., 2017). Moreover, mineral dust interacts with liquid or ice clouds, modifying their optical properties and lifetimes (DeMott et al., 2003), and affecting the precipitation processes (Creamean et al., 2013). Finally, the lifetime of the mineral dust particles affects the deposition fluxes over land and ocean (van der Does et al., 2018). Once dust particles are deposited at the surface, they provide micro nutrients to the ocean or to land ecosystems (e.g., Jickells et al., 2005).
All aforementioned processes are strongly affected by the dust Particle Size Distribution (PSD). The dust PSD changes rapidly after emission since the preferential gravitational settling quickly removes the larger particles. However, Earth System Models and transport models fail to properly simulate dust deposition. There have been comparisons of model simulation against measurements, that show a consistent overestimation of large particle removal (e.g., Ginoux et al., 2001;Colarco et al., 2002). On the other hand, dust model inter-comparison by Huneeus et al. (2011) found that, the model estimates of deposition vary over a large range, yet mainly underestimating the observations. Adebiyi and Kok (2020) showed that the mass of coarse dust (Diameter ≥5 μm) in the atmosphere is about four times greater than the simulated in current climate models, resulting in greater total dust mass load. This signifies the importance of proper modeling of the mineral dust deposition.
The conclusion that climate models underestimate the coarse mode dust in the atmosphere, is consistent with observations. Denjean et al. (2016) showed that the effective diameter of the Saharan dust coarse mode does not vary with dust age over the Mediterranean for transport times between 1.5 and 7 days. This conclusion has been drawn by consolidating high-quality aircraft measurements of dust PSDs performed during the last years in the framework of few large-scale experimental campaigns (Ryder et al., 2013a;Ryder et al., 2013b;Ryder et al., 2018;Ryder et al., 2019). Weinzierl et al. (2017) presented results of aircraft and ground-based measurements during the SALTRACE campaign in the tropical Atlantic in 2013/2014. Particles with sizes in the range of 10-30 μm were detected over Barbados, during a Langrangian flight sequence from Cape Verde, far beyond the Stokes gravitational settling expectations. Similar conclusions were also validated by deposition of giant particles (size greater than 75 μm), as depicted by samplers mounted on moored dust-collecting surface buoys, located in the Atlantic Ocean at 2,400 and 3,500 km from the west African coast (van der Does et al., 2018). These collected dust particles were mostly well-rounded quartz particles up to 450 μm in polar diameter, with what appeared to be high aspect ratios. Moreover, Goudie and Middleton (2001) reported that particles larger than 62.5 μm are commonly carried from Sahara to the British Isles. van der Does et al. (2018) discussed four potential different mechanisms that could facilitate long-range transport of large/ giant particles. These mechanisms are: 1) strong winds causing fast horizontal transport, 2) strong turbulence keeping particles in suspension for longer time (Garcia-Carreras et al., 2015), 3) electrical forces that balance the gravitational force (Nicoll et al., 2010;Renard et al., 2018;Toth et al., 2020), and 4) the presence of thunderstorms or tropical cyclones that can create strong uplifts. Additionally, the particle shape and orientation can also influence the dust particle transport time (Huang et al., 2020;Mallios et al., 2020;Mallios et al., 2021b). Specifically, Mallios et al. (2020Mallios et al. ( , 2021b showed that simulating dust particles as prolate spheroids, they can remain much longer in the atmosphere than their spherical counterparts, and that the particle vertical orientation can increase the residual time in the atmosphere significantly, compared to the horizontal one. Vertical orientation becomes possible in the presence of external large scale electric fields of strength in the order of 10 4 V/m (two orders of magnitude greater than the fair weather electric field values) (Mallios et al., 2021a).
The field of Atmospheric Electricity research can provide insights on the electrical processes and their contribution to the transport of dust plumes. The Earth's atmosphere is a conducting medium, attributed to the presence of ions, created by ionization by galactic cosmic rays radiation. Ions attach to dust particles through the processes of ionic diffusion and Coulomb interaction (Gunn, 1954;Klett, 1971) and the polarization due to the presence of an external electric field (Gunn, 1956;Klett, 1971), leading to their subsequent charging (Zhou and Tinsley, 2007;Zhou and Tinsley, 2012). This large scale electric field is created by the potential difference between the lower part of the ionosphere and the Earth's surface (Rycroft et al., 2008), and is modified by the ion density reduction caused by the ion attachment process, which leads to the electrical conductivity reduction and the atmospheric columnar resistance enhancement (Zhou and Tinsley, 2007;Baumgaertner et al., 2014). Moreover, the dust particles can be also charged during collisions, a mechanism known as triboelectric process (Kamra, 1972;Eden and Vonnegut, 1973), that results in large particle being predominantly positively charged and small particles negatively (e.g. Zhao et al., 2002;Lacks et al., 2008;Shinbrot and Herrmann, 2008;Merrison, 2012). This size preferential polarity acquisition in combination with updrafts or shape induced gravitational sorting (Yang et al., 2013;Mallios et al., 2020;Mallios et al., 2021b), can create charge separation that further enhances the electric field magnitude, similarly to processes that take place in thunderstorms (e.g., Chiu, 1978;Helsdon et al., 2002;Williams and Yair, 2006, and references therein).
Motivated by the studies of dust particle electrification and charging, we focus on dust charging due to ion attachment process and we develop an 1D numerical model for the of spherical dust particle settling in the presence of a large scale electric field. The model is able to calculate self consistently the modification of the atmospheric ion densities in the presence of the dust particles, and the consequent alteration of the atmospheric electrical conductivity and the large scale electric field. Moreover, it is able to evaluate the acquired electrical charge on the dust particles, and thus to calculate the electrical force that is applied on them. Finally, the effect of the electrical force on the gravitational settling and the terminal velocity is quantified.
In Section 2, the mathematical formulation of the model is presented and explained. In Section 3, results of the model under different conditions are presented and discussed. Finally, in Section 4 the results of this work are summarized.

System Dynamics
The bottom (z 0 km, Earth's surface) and the top (z 40 km) boundaries are assumed to be perfect conductors of electricity. Although the lower part of the ionosphere is typically at ∼70-80 km (e.g., Rycroft et al., 2008), the perfectly conducting boundary can be placed at 40 km (for simulation purposes), since the electric field at this altitude is mainly vertical, regardless of the electrical activity at altitudes below 12 km (Pasko et al., 1997).
Increasing the upper boundary altitude to 50 km, the relative difference between the results in the area of interest (up to 12 km altitude) is less than 1%. The electrical potential at the bottom boundary is set equal to 0 V, while at the top boundary is set to 250 kV, which can be seen as the average potential along the ionosphere (Rycroft et al., 2000).
In the region under investigation, the ion conductivity is dominant over the electron conductivity, as electron number density exceeds the bipolar ion number densities at altitudes above ∼60 km (Pasko et al., 1997). Therefore, the particle types that are considered in the presented formalism are the atmospheric small ions, and the dust particles in terms of atmospheric aerosol content.
As a first approximation in order to highlight the possible physical mechanisms, we neglect the wind speed, and we consider the atmosphere to be stagnant, which denotes to the particle velocity with respect to the air being the same with the particle velocity with respect to an observer at the ground, ⃗ vp. The continuity equation for positive and negative small ion number densities, n ± , is: where μ ± are the positive and negative small ion electrical mobilities, E ⃗ is the large scale electric field, D ± ion are the positive and negative small ion diffusion coefficients, q is the ion pair production rate, a is the ion-ion recombination rate, and β ± i are the positive and negative ion attachment rates to dust particles with radius r i . The first two terms in brackets are the flux of small ions due to their motion in response to the large scale electric field, and the ionic diffusion. The third term describes the production rate of positive and negative ions, while the last two terms express the ion losses due to recombination, and attachment to dust particles.
The current continuity equation for the charge transport of dust particles with radius r i is: where ρ p,i is the dust particle charge density, ⃗ vp,i is the dust particle velocity, D p,i is the dust particle diffusion coefficient, and q e is the elementary charge. We note that all physical constants appearing in the manuscript, are listed in Table 1. The first two terms in brackets are the charge flux due to the dust particles advection/settling (advection is a valid term for the 3D case where horizontal components of velocity exist, while the settling happens in the 1D case, where the particles fall in the vertical direction), and the diffusion. The third term describes the gain of charge due to the ion attachment.
The continuity equation for the transport of number densities, N p,i , that correspond to dust particles with radius r i , is: where the terms in brackets are the dust particle flux due to the advection/settling and the diffusion. Eqs. 2, 3 have the same particle velocities, and the same particle diffusion coefficients, since they refer to different properties of the same particle. The particle charge density and the particle number density are related through the expression ρ p,i q p,i N p,i , with q p,i being the particle net charge. These three quantities vary over time, but different equations govern their variation. The linkages between their temporal and spatial variations, are the common velocity and diffusion coefficient.
Finally, the large scale electric field, E ⃗ , is calculated from the total charge density, ρ tot , from the Gauss law: Although Eqs. 1-4 are generic, in this work we shall focus on the 1D case, in the vertical direction along the altitude. Applicability of the model can be easily extended to other types of aerosols that can be simulated by spherical shapes and act as good conductors of electricity. Substituting E ⃗ in Eq. 4 with − ⃗ ∇Φ, the Poisson equation is derived, which is solved using a Full Multigrid Algorithm (FMG) (Press et al., 1992, p. 877). The advection terms in Eqs. 1-3 are discretized as an average between the second order Upstream Nonoscillatory (UNO2) advection scheme (Li, 2008) and the third order Quadratic Upstream Interpolation for Convective Kinematics with Estimated Upstream Terms (QUICKEST) method (Leonard, 1979). This approach reduces the artificial numerical diffusion significantly, without using computationally demanding flux limiters or flux corrected transport (FTC) schemes. The diffusion terms in Eqs. 1-3 are discretized as a second order central difference scheme. The simulation domain is divided in 2049 mesh points, and increasing further their number by a factor of two leads to a relative difference less than 5% between the results. The results obtained with the aforementioned approach and the given mesh points, are very close (relative difference less than 2%) to the results obtained with a first order upwind scheme (Press et al., 1992, p. 834) with 8,192 mesh points. This shows that the used scheme is indeed less diffusive than a first order upwind scheme under the same conditions, and also that the averaging of two schemes can be a valid procedure. After all, the principle of a weighted average between a stable but diffusive low order scheme with an unstable but less diffusive higher order scheme is the core of the FTC schemes (Zalesak, 1979). The FMG algorithm ensures fast convergence of the Poisson equation solution, when the number of mesh points becomes very large (as in our case) (Press et al., 1992, p. 871).

Meteorological Conditions
The U.S. Standard Atmosphere 1976 (NOAA/NASA/USAF, 1976) is considered as the standard static atmospheric model, and the vertical distributions of pressure, P, and temperature, T, are constant over time ( Figure 1). The air mass density, 9 air , at a given altitude can be calculated from pressure and temperature, provided that the atmosphere is an ideal gas: where R* is the specific gas constant for dry air.

Ionization Rate
There are three mechanisms that create the atmospheric ions (Tinsley and Zhou, 2006, and references therein): 1) ionization by natural radioactivity originating in the ground, including direct α, β, and γ radiation from the surface layers and dust aerosol, and radiation from radioactive gases (principally 222 Rn but also 220 Rn) and daughter products that are carried up into the troposphere by vertical convection. 2) In the troposphere the dominant ionization mechanism is the flux of galactic cosmic rays (GCR). This does not apply in polar and subpolar regions, where the ionization due to solar energetic particles is important (Mishev, 2013). 3) In the stratosphere there are additional ion pair sources such as the flux of relativistic electrons (of a few MeV energy) precipitated from the radiation belts and peaking at subauroral latitudes, and the solar energetic particle events (mainly protons) in the polar and subpolar regions (Mishev, 2013).
In the current study the ionization rate as formulated by Tinsley and Zhou (2006) is adopted with two simplifications. The first is to neglect the natural radioactivity effect originating in the ground. This leads to a cruder but smoother vertical profile of the ion densities and the electrical conductivity along the altitude. Moreover, such a choice eliminates the variation of the ionization rate between continental and marine areas (as we show in Section 2.10), and thus the only dependence that is left is on the geomagnetic latitude. This leads to an approximation of the ionization rate which is more suitable for highlighting the mechanism of the particle-ions interactions, which is the focus of the work, and interpreting the obtained results.
The second simplification is the creation of a generic temporal ionization rate profile, by averaging the solar minimum and solar maximum profiles given in the aforementioned work, leading to a profile independent of the solar cycle. This resulting empirical profile, although admittedly simple, is suitable since it serves as the baseline approach for the investigation of the ionization rate distribution effects on the ion density and, consequently, the atmospheric electrical conductivity distribution. Indeed, for in depth analysis more detailed profiles would be a better representation, regarding e.g. the annual variation, however, this goes beyond the scope of the paper at hand.
The mathematical formalism for the ionization rate, q, is listed in Table 2. The quantities q i and z i that depend on the geomagnetic latitude, and are used in the formulas, are plotted in Figure 2.

Recombination Rate
The ion-ion recombination rate, a, is derived by fitting the values of Bates (1982).
where [M] is the concentration of air molecules in units of 2.69 × 10 19 cm −3 . In the case of US Standard Atmosphere, [M] in units of m −3 can be calculated as: where P is the pressure, k is the Boltzmann constant, and T is the temperature.

Ionic Electrical Mobilities
The ion electrical mobilities μ ± are scaled along the altitude as (Zhou and Tinsley, 2007): where μ ± 0 is the electrical mobility of positive and negative ions at reference pressure, P 0 , and reference temperature, T 0 , and S is the Sutherland's constant. In the present work, μ ± 0 are adopted by Bricard (1965), and are listed in Table 1.

Diffusion Coefficients
The positive and negative small ion diffusion coefficient, D ± ion , is given by the Einstein's relation: The diffusion coefficient for dust particles of radius r i , D p,i , is given by the Einstein's relation: where μ m,i is the mechanical mobility, given by the Cunningham-Knudsen-Weber-Millikan equation (e.g., Tammet, 1995), η is the air dynamic viscosity, and C c is the slip correction factor. The dynamic viscosity, η, is scaled along the altitude according to Sutherland's model (e.g., Zhou and Tinsley, 2007, and references therein): where η 0 is the viscosity at reference temperature T 0 . The slip correction factor is adopted by Davies (1945) as: where l is the mean free path in air that is given by the expression (Jennings, 1988):

Particle Velocity
The magnitude of the velocity, v p,i , of a particle with radius r i at any time point is calculated by solving the Newton's second law of motion: Altitude range (km) q (pairs/(m 3 s)) Scale height (km) where m p,i is the particle mass, F ⃗ el is the electric force that acts upon the particle, F ⃗ grav is the gravity force, and F ⃗ drag is the drag force.
In 1D case, and assuming the positive direction to be upwards, Eq. 14 becomes: with m p,i being the particle mass: and 9 p being the mass density of a dust particle. When the electric force, q p,i E, is positive (negative charge in a downward pointing electric field, or positive charge in an upward pointing electric field) and greater than gravity, m p,i , the particle moves upwards (positive velocity or acceleration) and the drag force is negative pointing downwards in the same direction as gravity ( Figure 3A). When the electric force is positive but less than gravity, the particle moves downwards (negative velocity or acceleration) and the drag force points upwards ( Figure 3B). Finally, in the case of negative electric force (negative charge in an upward pointing electric field, or positive charge in a downward pointing electric field), the particle moves downwards, and the drag force points upwards ( Figure 3C).
The magnitude of the drag force is: where A p,i is the projected area of the assumed spherical dust particle, and C d is the drag coefficient, adopted by Clift and Gauvin (1971): with Re being the Reynolds number, defined as: It is noted that in both Eqs. 17, 19 the magnitude of the particle velocity v p,i is used. According to the definition of the Reynolds number and the drag force, the magnitude of the particle velocity with respect to the air is needed. Since in our analysis we neglect the effects of the wind, and therefore the wind speed is set to zero, the velocity of the particle with respect to an observer at the ground, v p,i , is the same as the velocity of the particle with respect to the air. The choice of Eq. 18 as an expression for the drag coefficient comes from the fact that it is valid for Reynolds numbers up to 10 5 , allowing the study of very large particles. Moreover, it is within 6% of experimental measurements (Clift et al., 2005). Additionally, it has been shown by Mallios et al. (2020) that, Eq. 18 compared to the drag coefficient derived by Proudman and Pearson (1957), which is the most accurate analytical expression for Re < 0.1, gives a relative error less than 2%. Finally, it is noted that since the slip correction factor, C c , has been introduced for the correction of the Stokes law (and therefore it is valid only in the Stokes region), it appears in Eq. 17 only in the case of Re < 0.1 (Mallios et al., 2020).
Finally, in order us to avoid very small time steps for the solution of Eq. 15, that will lead to large simulation times, we apply the following procedure for the calculation of the particle velocity, v p,i . The time step of the simulation is calculated based on the processes of Eqs. 1-3. Eq. 15 is discretized as: where F drag is calculated using the velocity at time point t, v t p,i . Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 709890 If F drag > q p,i E − m p,i g , it means that in the period of dt the particle has reached its terminal velocity. In this case the velocity is calculated by solving the equation: using the Brent's method (Press et al., 1992, p. 359).
In the case that F drag ≤ q p,i E − m p,i g , the terminal velocity has not been reached, and Eq. 20 is solved.

Ion Attachment Rates
For the quantification of ion attachment to dust particles, the particles are assumed perfect electrical conductors (Fuks, 1958, p. 81).

Continuum Regime
The formulation of the attachment rates that is adopted in this work, is valid in the continuum regime. This means that the ionic mean free path, l ± ion , is smaller than the particle radius, r i . Defining the particle diffusion Knudsen number, Kn ± D , as (Dahneke, 1983): where < v > ± is the average thermal velocity, and m ± is the mass of positive and negative ions respectively, the continuum regime is defined for Kn D < 1.
The ionic mass, m ± , can be correlated with the ion electrical mobility, μ ± 0 , at reference pressure, P 0 , and reference temperature, T 0 . The only available data set covering the mass range up to 2,122 amu has been derived by Kilpatric (1971), in measurements of molecular ion mobility in nitrogen, at a pressure of 1,013.25 hPa and temperature of 473 K. Meyerott et al. (1980) reproduced the results of Kilpatric (1971) as the most complete data for big cluster ions. Tammet (1995) based on the dataset of Kilpatric (1971) derived an ion mass-reduced mobility correlation at temperature equal to 273 K, and pressure equal to 1,013.25 hPa. Scaling the data set of Kilpatric (1971) to reference temperature, T 0 , and reference pressure, P 0 , the following ion mass-reduced mobility correlations are derived by numerical fitting the data points: where m ± are given in kg. It is noted that the mean absolute error of Eq. 23 is 6.8%, and of Eq. 24 is 6.3% (both less than 10%). Although these errors are larger than the mean absolute error of the expression derived by Tammet (1995), which does not exceed the value of 2%, they are much simpler to use and can be seen as a very good compromise between simplicity and accuracy.
Using the adopted values of μ ± 0 by Bricard (1965) in Eqs. 23 and 24, we get m + 4.486 × 10 -25 kg , and m − 2.124 × 10 -25 kg . Substituting these values to Eq. 22, it is found that, the continuum regime condition is satisfied for r i ≥ 0.2 μm up to altitude equal to 12 km, which is the typical limit between the upper troposphere-lower stratosphere. In the case of altitudes up to 20 km, where the lower edge of the stratosphere near the equator can be, the continuum regime condition is satisfied for r i ≥ 1.2 μm.

Physical Mechanisms
According to the field-diffusion theory, two major mechanisms act simultaneously for the ion attachment to particles. The first is the diffusion to the particle surface, and the second is the conduction in the electric field (Helsdon et al., 2002). The electric field is the superposition of the external large scale electric field, the induced dipole electric field, and the electric field due to charges on the particle (Lawless, 1996). In the presence of bipolar ion, particle charge approaches a steady-state level as charging time increases (Fjeld and McFarland, 1989). Gauntt et al. (1984) showed that neither the diffusion nor the field approximation alone evaluate correctly the steady state particle charge in the case of q e r i E ⃗ | | kT > 1, but the superposition of the two mechanisms, while not rigorous, yields reasonably accurate results (Gauntt et al., 1984). Hence, this method is suggested as a practical means for calculating charge acquired under bipolar charging conditions, and is adopted in the present work, since to the best of our knowledge up to now, no rigorous analytical solution of field diffusion theory equations for stationary or moving particles have been put forth.

Stationary Particles
In the case of stationary particles, the diffusion approximation is obtained by neglecting the effect of the external electric field on the particle and on the ion transport, while retaining the effect of the particle charge (Fjeld and McFarland, 1989). The total rate per unit volume of positive and negative ion loss due to their attachment to dust particles of radius r i , and number density N p,i is (Gunn, 1954): Φ s is called self potential and indicates the influence of the particle charge to the diffusion current that flows to the particle. When Φ s < 1, the thermal energy of the ions, kT/q e , is larger than the electric potential energy of the particle due to its charge. This means that ion diffusion dominates and the ion mobilities are the important factors that determine the type and the number of the ions that attach to the particle. When Φ s > 1, the surface electric field magnitude of the particle is larger than the ion thermal energy, and the ion drift in the presence of the particle surface electric field dominates. In the limit that Φ s ≪ 1, Eq. 25 becomes: The field approximation is obtained by neglecting diffusion, and calculating the charging rate from the transport of ions along the electric field lines to the surface of the particle (Fjeld and McFarland, 1989). The electrical forces due to the images of the ions with respect to the conducting particles have been neglected, because as has been shown by Marlow and Brock (1975), their effect is not important in the continuum regime. The Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 709890 total rate per unit volume of positive and negative ion loss due to their attachment to dust particles of radius r i , and number density N p,i is (Gunn, 1956): where q F 12πε 0 E ⃗ r 2 i is the limiting charge value for which the particle is considered to be totally positive or negatively charged.
Several approaches have been proposed in the past for the superposition of these two mechanisms, both in theoretical and model simulation basis. Several theoretical expressions that combine the mechanisms or that numerically solve the fielddiffusion problem have been proposed in the past (Klett, 1971;Smith and McDonald, 1975;Liu and Kapadia, 1978;Fjeld et al., 1983;Wang, 1983;Fjeld and McFarland, 1989;Lawless, 1996). Unfortunately, these expressions cannot be easily expanded to other cases such as moving particles, or different shapes, because they contain parts of approximate solutions that restrict the generalization. Moreover, some of them have convergence difficulties, or are valid in cases and conditions that are not commonly used (Lawless, 1996). On the other hand, a simple superposition of Eq. 27 with Eqs. 28, 29 (Chiu, 1978) or Eq. 25 with Eqs. 28, 29 (Helsdon et al., 2002), although not rigorous, offers the flexibility of generalization, because each individual charging mechanism has been studied extensively on a broad spectrum of applications.
The disadvantage of the superposition used by Chiu (1978) is that, in the case that q p,i ≥ q F , where the particle is not polarized, and therefore mainly the diffusion and the electric field created by the particle charge influence the attachment process, the expression does not give the Eq. 25. On the other hand, the superposition used by Helsdon et al. (2002) has the disadvantage of accounting for the electrical effect due to the particle net charge twice (both in diffusion and conduction mechanisms). In order to tackle these issues we propose the following superposition: Eqs. 30, 31 in the case that the particle is not polarized give the solution of the diffusion approximation, and in the case that the particle is polarized lead to the superposition used by Chiu (1978). Moreover, the effect of the electric field due to the particle net charge is not counted twice as in superposition by Helsdon et al. (2002). Long and Yao (2010) presented an evaluation of nine different particle charging models based on existing experimental results. One of the conclusions was that the field-diffusion combined model developed by Lawless (1996) should be the first choice relatively for numerical models of the particle dynamics. Since this model is shown to be the most preferable in terms of the theoretical consideration, accuracy and computational time, we use it as a reference for comparison with our proposed superposition (Eqs. 30, 31). The results are shown in Figure 4, where the attachment rates for positive ions, β + , are plotted as a function of particle charge, for different values of external electric field (the same conclusions are derived for the negative ions attachment rates). In Eqs. 25, 27-31 as attachment rate, β ± (with units 1/s), the expression in the right hand side that is multiplied with the appropriate ion number density, n ± , is defined (leading to an expression of the form dn ± dt β ± n ± ). Figure 4A shows the relative difference as a function of normalized particle charge, for different particle radii, in the presence of external electric field with magnitude 10 2 V/m. For r i 0.2 μm the mean relative difference in the studied range of the particle charge is 0.12 %. For r i 2 μm it is 1.08 %, for r i 20 μm it is 8.09 %, and for r i 200 μm it is 9.02 %.
The maximum particle charge is chosen to be 1,000 elementary charges, which is in the same order of magnitude with the charge derived from the charge density measured by Kamra (1972), and assuming particle density equal to 10 8 particles/m 3 , which is the annual mean of the lower limit of aerosol concentration over oceans (Tinsley and Zhou, 2006). The specific measurement was conducted over dust devils, which is an extreme phenomenon compared to the regular dust transport (in such conditions Nicoll et al. (2010) measured three orders of magnitude lower charge density) and the particle density is at the lower limit of particle densities observed in detached dust layers. For larger particle densities, keeping the charge density the same, the individual particle charge decreases. Therefore, the 1,000 elementary charges per particle can be used as an upper limit for the particle charge. Figure 4B shows the mean relative difference for a range of particle charges −1,000 to 1,000 elementary charges, as a function of particle radius, in the presence of external field with various magnitudes. The mean relative difference does not exceed the value of 50 % for a range of electric field magnitudes from 0 to 10 5 V/m. Therefore, the proposed superposition of Eqs. 30, 31 can be considered as physically Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 709890 reasonable, and although they are not strict analytical solutions of the ion transport equation that describes the ion attachment process, they can be used in numerical models to describe adequately the effects of the ion attachment to particles. The maximum electric field magnitude in Figure 4B is chosen to be in the order of magnitude of observed electric field values right before a lightning flash (e.g., Marshall et al., 2005;Stolzenburg et al., 2007). Regarding dust layers transported for days in the Earth's atmosphere, away from the emission source, no case of lightning discharge has been reported, indicating that, in these cases, the value of 10 5 V/m is unlikely to be reached or exceeded. Thus, it can be used as an upper limit of the electric field magnitude. On the other hand, in cases of dust storms or dust devils that occur on Earth (or other planets, e.g., Mars) and in the vicinity of the emission source, the value of 10 5 V/m (or the equivalent electric field threshold on that planet) can be reached. This can happen due to the intense electrification caused by the triboelectric effect, which is the dominant electrification mechanism during the emission of dust particles in the atmosphere, and the large scale charge separation caused by the strong updrafts (Krauss et al., 2003). A short review on the electrification of dust storms and on the resultant electric field strengths has been made by Riousset et al. (2020). When the electric field threshold is reached, or exceeded, a discharge occurs, leading to a decrease of the electric field magnitude. Therefore, 10 5 V/m acts again as an upper limit of the electric field magnitude.

Moving Particles
Eqs. 30, 31 can be easily transformed in the case of moving particles. According to Gunn (1954), the diffusion approximation for moving particles leads to the following expression for the total rate per unit volume of positive and negative ion loss: where the term A ± v represents the ventilation effect on the ion attachment process.
The coefficient F is an experimentally derived coefficient to match the theoretical results with actual measurements. It depends on the detailed air flow over the spherical particle, and can be described as a fitting function of the Reynolds number, Re, based on measurements by Kinzer and Gunn (1951) as: , Re ≤ 0.5 0.079Re 3 − 0.899Re 2 + 3.156Re − 1.289 , 0.5 < Re ≤ 5 2.218Re −0.169 , 5 < Re ≤ 423 5.165×10 −7 Re 2 − 2.828 × 10 −4 Re + 0.882 , Re > 423 The field approximation for moving particles has been studied by Whipple and Chalmers (1944), based on the theory of charge acquisition by particles falling through a region containing ions proposed by Wilson (1929). The total rate per unit volume of positive and negative ion loss is: FIGURE 4 | Comparison between Eq. 30 and formulation by Lawless (1996): (A) Relative difference as a function of normalized particle charge, for different particle radii, in the presence of external electric field with magnitude 10 2 V/m, and (B) mean relative difference for a range of particle charges −1,000 to 1,000 elementary charges, as a function of particle radius, in the presence of external field with various magnitudes.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 709890 where ⃗ vp,i↑↑E ⃗ means that the velocity vector and the electric field vector are in the same direction (the particle moves in the same direction as the electric field lines), while ⃗ vp,i↑↓E ⃗ means that the velocity vector and the electric field vector are in opposite directions (the particle moves in the opposite direction from the electric field lines).
Although Eqs 36, 37 have been derived under the stream-line flow assumption, it has been discussed by Whipple and Chalmers (1944) that the provided formalism is appropriate even in turbulent environments. As a matter of fact, Gott (1936) gives experimental evidence to show that turbulence does not affect the general validity of Wilson's theory.
A superposition similar to Eqs 30, 31 leads to the following expression for the combined diffusion-field theory approximation in the case of moving spherical particles:

Dust Particle Number Density Distribution
Several aircraft campaigns over recent years have attempted to measure the properties of Saharan mineral dust. Ryder et al. (2018) and Weinzierl et al. (2017) provide lists of the major aircraft campaigns with a focus on in-situ dust measurements. In the current work, we use the mean size distribution between altitudes of 1-6 km at STP, measured by Ryder et al. (2013b) during the Fennec 2011 aircraft campaign at Mauritania and Mali, since they represent measurements very close to the dust particle emission sources ( Figure 5A).

Size Bins Apportionment
The PSD discretization is of high importance for the correctness of the results. The number of bins, along with the size of each bin, must ensure that the initial measured distribution is properly introduced to the model without significant loss of information, deformation, or non convergence of the final results. The adopted size distribution is discretized using 19 bins, listed in Table 3. Figure 5B provides the representation of the discretized number density distribution from Fennec, using the aforementioned binning, along with 5, 10, 20 equilateral bins, respectively. Figure 5C shows a comparison between the measured distribution, with the reconstructed ones from the different types of binning. It becomes apparent that as the number of the equilateral bins increases, the initial distribution is better approximated. Nevertheless, a significantly large number of bins is required for the proper approximation of the whole range of the distribution (judging from the slow convergence at the small particles, more than 40). On the other hand, the adopted binning of the current work approximates very well a wider range of the initial distribution using a fairly small number of bins (19 bins). Furthermore, Figure 5D illustrates the relative difference between the reconstructed distributions and the initial one. As the number of bins increases and their length becomes much smaller than the mean radius of each bin, the relative difference with the original distribution decreases.
For the current work binning, the agreement becomes ∼1% for the narrow bins, while for the larger bins the relative difference increases to values comparable to other types of binning. Moreover, the relative difference is better distributed in a wider range of particle radii. Therefore, it is proven that the adopted binning leads to a better introduction of the initial distribution to the model, especially the bins that have length less than 1% of the mean radius (even numbered bins). In the Results section, we demonstrate that the electrical properties, calculated using the adopted binning, are the convergence limits of the values obtained by increasing the number of equilateral bins.

Spatial Particle Number Density Distribution
The dust particles are introduced in the simulation domain using two different spatial distributions. The first one is a rectangular distribution of the form ( Figure 5E): where h I i 3 km is the altitude of the distribution center, d I i 2,500 m is the distribution depth, and b z 3d z is a smoothing factor regarding the number of mesh points that will sample the boundary of the distribution (in our case 3). The two terms right before the bracket convert the number density N 0,i from STP to the model ambient temperature and pressure conditions at a given altitude, z.
The second is a Gaussian distribution, mathematically expressed as ( Figure 5F): where h II 5,000 m is the altitude of the distribution center, and d II i 250 m is the distribution depth, defined as the distance between the altitudes that the number density is reduced to half of its maximum value (full width at half maximum of a Gaussian distribution). The factor 12.6536 ensures that the integral of N II p,i along the altitude gives the same result as the integral of N I p,i along the altitude. This means that both distributions lead to the same surface particle density, and consequently to the same total number of particles.
The number density distribution N I p,i is similar with the one measured by Ryder et al. (2013a), while the distribution N II p,i denotes a study case where the particles are more concentrated, with higher number density, leading to different electrical properties.

Initial Fair Weather Electrical Conditions
At t 0 s all the dust particle bins are placed in the simulation domain, under fair weather conditions and are left to evolve under the influence of gravity, the electrical and drag forces. The initial fair weather electrical conditions are obtained by solving Eqs 1, 4 in the absence of dust particles, until steady state is reached. By doing that, the ambient fair weather ion densities, n ± fw , electrical conductivity, σ fw q e n + fw μ + + n − fw μ − , total charge, q tot, fw q e n + fw − n − fw , and electric field magnitude, E z, fw , are derived.
Since the ionization rate, q, which is the source of atmospheric ions, depends on the geomagnetic latitude, and therefore the latitude and longitude, three different locations have been chosen for the calculation of the fair weather electrical conditions, and the study of their spatial variation. These are: 1) Mauritania (MR, 24.316°, −10.962°), 2) Cape Verde (CV, 16°, −24°), and 3) Barbados (BB, 13.167°, −59.33°). The three locations have been chosen because they belong to the route of the transatlantic transport of Saharan dust from Africa to Caribbean (Weinzierl et al., 2017;van der Does et al., 2018). Figure 6 shows a comparison between the initial fair weather electrical conditions of the aforementioned locations. Figures  6A,B show that the ionization rate in Cape Verde and in Barbados is practically the same, with the relative difference being less than 1%, leading to the same fair weather electrical conditions.
On the other hand the ionization rate in Mauritania is higher, and can reach a maximum relative difference value of 40% compared with the ionization rate in Cape Verde. Regardless this significant difference in the ionization rates, the relative differences in the fair weather ion densities (Figures 6C,D) and the fair weather total charge densities ( Figures 6E,F) do not exceed the value of 20%. At this point we note that the small difference between the positive and negative ion densities allows only one field to be graphed in Figure 6C. In the case of the fair weather electric field, the relative difference is even smaller, and does not exceed the value of 10% (Figures 6G,H). Therefore, the ionization rate dependence, solely, on the geomagnetic latitude does not create significant variations of the initial fair weather electrical conditions along the transatlantic transport of the dust particles. This means that, for the given simplistic ionization rate formulation it does not make a difference on the choice of the reference location between Mauritania and Barbados. In the simulation results presented in the current work, the fair weather electrical conditions of Mauritania are used as initial electrical conditions, because this location is closer to the measured dust particle size distribution that is used in the model. At this point we want to emphasize that, in reality there are additional spatial and temporal variation of the ionization rate that for the sake of simplicity have been omitted, which will lead to larger fair weather conditions differences along the transatlantic passage than the ones presented here, and their impact on the electrical evolution of the dust particles are planned to be studied in future work using a higher dimension model. Figure 7 shows the time dynamics of the electrical properties of Bin 18 (particles with radius, r 100 μm), as well as the ambient conductivity and the large scale electric field magnitude. The ambient conductivity, σ, decreases, because the atmospheric ions attach to the total number of particles regardless their size, leading to a reduction of the ionic densities ( Figure 7A). The conductivity reduction within the dust layer results in an increase of the columnar resistance, and a decrease in the conduction current at the ground (Baumgaertner et al., 2014). Since the conductivity at the ground does not change from the fair weather value, the electric field magnitude, depicted in Figure 7B, decreases to reflect the reduction in the conduction current (Daskalopoulou et al., 2021). On the other hand, as the system evolves towards the steady state, the conduction current J z σE z tends to become constant along the altitude, and therefore when the conductivity decreases the electric field magnitude increases ( Figure 7B).

Particle Charging and Electrical Force
At time t 4 2,668 s, more than 70% of the Bin 18 particles have vanished ( Figure 7C). As these particles fall to the ground their charge density increases since more ions attach to them ( Figure 7D). Especially in the region below the main dust layer (below 1 km), where the ion density increases to the fair weather values, the falling dust particles encounter a larger pool of ions available for attachment than the one within the layer. This consequently leads to an increase of the particle charge ( Figure 7E). The dust particles are charged negatively, since more negative ions attach to them than positive ones, due to their higher mobilities. To be more precise, the ion attachment process depends both on the ion mobilities and the ion densities (Eqs 30,31,38,39). If the mobility ratio, μ − μ + is larger than the ion density ratio, n + n − , the ion mobilities influence the attachment process, otherwise, the ion densities do that. In our case, the ion mobility ratio is equal to 1.357, while the fair weather ion density ratio ranges between the values 1.0014 and 1, depending on the altitude, and thus the ion mobilities characterize, initially, the type of ions that attach to particles. As the negative ion density decreases faster than the positive ion density, the ion density rate increases up to the value where it becomes equal to the mobility ratio. Further increase of the ion density ratio, leads the positive ions to attach faster to the particle than the negative, resulting a reduction of the particle negative net charge, or even the charging of the particle with positive charge. In the case presented in Figure 7, the ion density ratio does not become larger than the mobility ratio, and all the particles of all bins have negative net charge. Regardless of the charge and the large scale electric field magnitude, the electrical force that is acted upon the particles is negligible compared to gravity. As it can be seen in Figure 7F, the ratio between the electrical force and the gravity, R F,18 is less than 2 × 10 -6 . Relative difference between ionization rates, (C) fair weather ion densities, (D) relative difference between fair weather ion densities, (E) fair weather total charge densities, (F) relative difference between fair weather total charge densities, (G) fair weather electric field magnitudes, (H) relative difference between fair weather electric field magnitudes. For a better representation of the bin particles dynamics, and so as to efficiently display the electrical properties of the other bins as well, we define the following average quantities. The average particle charge of a bin, <q p,i >, at every time step is: where z max,i , and z min,i , are the maximum and the minimum altitudes of the particle distribution that belongs to Bin i, respectively. The average electrical force magnitude of a bin, <F el,i >, at every time step is: The average velocity magnitude of a bin, <v p,i >, at every time step is: The average altitude of a bin, <z p,i >, at every time step is: The usefulness of the defined average quantities can be described as follows. At t 0 s, all the particle bins are placed at the same altitude in the simulation domain, and they form the dust layer.
Since each bin consists of particles with different sizes, it moves with different velocity. As time progresses, each bin starts being detached from the dust layer and moves towards the ground. At each time step depending on the velocity, each bin is located at different average altitude, < z p,i >, in the presence of different electric field magnitude, E z . During their motion, the particle bins are not deformed, because the magnitude of the forces acted on them is not sufficient to cause any deformation (e.g., Figure 7C). This results that when the bins are at the same average altitude, the same average electric field is present, and the only quantities that are different are the average charge, the average electrical force magnitude, and the average velocity. By plotting these properties as a function of the average altitude, we can compare them under the same conditions, and their differences are only due to the particle bin size, and the bin number density. Figures 8A,B show the time evolution of the conductivity and electric field magnitude ratios with respect to their fair weather values, at different altitudes. As the conductivity decreases with respect to the fair weather value, the electric field magnitude increases, with two exceptions. The first one is at z 6 km, which is close to the upper boundary of the layer. As layer boundary we define the region where the dust particle number density gradually decreases to values at least two orders of magnitude less than the values inside the layer (in the case of the layer upper boundary), or the region where the dust particle number density gradually increases to the maximum values inside the layer (in the case of the layer lower boundary). In this region as the conductivity increases to the initial fair weather values, the electric field decreases. Similar behavior has been derived by Zhou and Tinsley (2007), and by Harrison et al. (2020) in layer clouds. This has to do with the reduction of the current density due to the increase of the columnar resistance caused by the conductivity reduction inside the layer (Rycroft et al., 2008). As the system approaches the steady state, the current density along the altitude tends to become constant (which is the steady state condition inside conductors (see e.g. Landau and Lifshitz, 1963, p. 86)). Since the potential difference between the boundaries of the simulation domain is constant, the increase of the columnar resistance caused by the conductivity reduction in the layer, leads to a decrease of the current density. In the regions where the conductivity remains the same to the initial fair weather values, decrease of the current density results in a decrease of the electric field. Consequently, in the upper boundary of the dust layer, as the conductivity increases to fair weather values, the electric field decreases with respect to its fair weather value to reflect the decrease of the current density.
The second is at the ground (z 0 km), where initially the electric field decreases, but as time progresses and several bins arrive at the ground it increases. The initial decrease of the electric field, right after the placement of the dust layer, seems to contradict measurements and observations that show significant increase when dust is present (e.g. Katz et al., 2018;Zhang and Zhou, 2020).This can be explained, first of all, by the fact that the 1d model is not appropriate for the study of effects upon the arrival of the dust layer, since it assumes an infinite symmetrical horizontal extent of the layer. The particles are placed instantaneously above a location, and are left to move vertically according to the forces acting upon them. Effects caused by the horizontal transport of the layer that can create horizontal electrical currents that eventually close to the ground enhancing the fair weather electric field (due to the electrical structure of the layer or the horizontal wind) require 2d or 3d modeling for a more detailed study. Especially the presence of horizontal winds, that are responsible for the dust particle transport, can alter the fair weather electric field by modifying the electrical structure of the atmosphere, as has been shown by Anisimov et al. (2001). As a matter of fact, the differences between the measurements, presented by Katz et al. (2018), during day and night at Wise Observatory and Mt. Hermon are attributed, among other factors, to the effect of the winds at the lower layer.
While the dust particles are aloft, and in the absence of any charge separation within, the electric field at the ground decreases as a consequence of the increase of the columnar resistance as has been shown by Harrison et al. (2020) and Daskalopoulou et al. (2021). As dust particles approach the ground and the electrical conductivity decreases, because the dust particle number density increases, the electric field increases. The increase is not large because the calculation is terminated when 99.9% of particles with radius 5 μm vanish from the domain (hence smaller particles remain more than 5 days in the atmosphere, having terminal velocities in the order of a few cm per hour, and their settling is so slow that we choose not to include it in the current study, for the sake of saving simulation time). Small particles (that are the majority of the dust particle population) are still in the domain and above the ground. In the case of Katz et al. (2018) and Zhang and Zhou (2020) the measuring instruments (field mills) are well within the dust layer that reaches the ground. The increase of the electric field inside the layer and when it is already present has been reported by the measurements and agrees with our findings. Differences in the electric field polarities can exist as a consequence of multiple factors that are neglected in the current model, such as wind or other electrification mechanisms, which are planned to be studied in future work. Finally, we note that, although the conductivity decreases by a factor up to 0.85, the electric field enhancement factor is no greater than 1.08. Figure 8C shows the relation between the number of elementary charges per particle and the particle size, at different altitudes. As the particle size decreases, particle charge decreases also as it is highly depended on particle size. Additionally, smaller particles have greater number densities than the large ones, and the attached ions tend to be divided on more particles leading to a decrease of the attached number of ions per particle (particle competition described by Gunn (1954)). Moreover, the charge of the larger particles increases as they travel within the atmosphere, as the attachment process depends on the particle velocity and the generated flow pattern, described FIGURE 8 | Characteristics of particles distributed along the altitude according to N I p,i : (A) Conductivity reduction with respect to the fair weather value at different altitudes over time, (B) electric field magnitude modification with respect to the fair weather value at different altitudes over time, (C) average number of elementary charges per particle as a function of particle radius for different average bin distribution altitudes, (D) ratio between average electrical force magnitude and gravity magnitude as a function of particle radius for different average bin distribution altitudes, (E) average particle velocity as a function of particle radius for different average bin distribution altitudes, (F) residual time of particles in the atmosphere until 99.9% of their concentration is lost.
by the Reynolds number. The larger the particle velocity the larger the difference between the particle charge at different altitudes. Figure 8D illustrates the dependence of the ratio between the electrical force and gravity, on the particle size. As the particle size decreases both electrical force and gravity decrease (the first one because the particle charge decreases, and the second due to the size reduction). The gravitational force reduction is larger than the particle charge reduction, and the force ratio increases, but does not overcome a value of 0.01. Additionally, as large particles settle, the force ratio increases due to the subsequent increase of the particle charge with velocity, which leads to an increase of the electrical force.
The velocity of the particles is a result mainly of the gravity and the drag forces. It does not charge significantly over the altitude and decreases as the particle size decreases ( Figure 8E). The residual time in the atmosphere, until 99.9% of the particle concentration of each size falls to the ground, spans in the range between 1 h, for particles with radius 100 μm, up to ∼ 120 h (5 days) with particles of radius equal to 5 μm ( Figure 8F) The later results are consistent with the calculations presented by Mallios et al. (2020Mallios et al. ( , 2021b.

Dependence on the Layer Depth
In this case the electrical properties of the distribution N II p,i are examined. Figures 8D-9A show the time evolution of the conductivity (Figures 9A,B) and electric field magnitude ( Figures 9C,D) ratios with respect to their fair weather values at different altitudes. The number density in this case is almost one order of magnitude higher than the N I p,i , and the conductivity reduction inside the layer is reduced by a factor of ∼0.25, almost 3 times lower than the previous case. The electric field magnitude is increased similarly by a factor of 3.5. In the region outside the layer, the conductivity reduction/electric field enhancement at each altitude depends on the number density of each bin located at that altitude at any given time point (Figures 9B,D).
The higher number density in the distribution leads to a more significant particle competition for ions. This leads to lower charge per particle than in the previous case when the particles are inside the layer. This can be seen clearly in Figure 9E where the number of elementary charges for particles that are located at average altitude equal to 4.9 km (in the center of the layer), is lower than in the previous case (dashed blue line in Figure 9E). As the particles keep falling the velocity influences the attachment process and increases the particle charge up until an average altitude of 3.2 km, where the particle charge remains constant for each size, regardless of the altitude. Similar remarks hold for the case of the ratio between the electrical force and the gravity ( Figure 9F). The velocity does not vary significantly over the altitude ( Figure 9G), and since the electrical force is still insignificant compared to gravity, the particle velocity is the same as the previous case. The difference in the spatial distribution leads to a difference in the residual time that individual particles remain in the atmosphere by 10% ( Figure 9H).

Dependence on the Particle Number Density
In order to demonstrate the charging mechanism dependence on the particle number density, we increase N I p,i by a factor of 10, leading the number density to become comparable to the fair FIGURE 9 | Characteristics of particles distributed along the altitude according to N II p,i : (A) Conductivity reduction with respect to the fair weather value at different altitudes over time, (B) magnification of panel (A) in the region z < 5 km, (C) electric field magnitude modification with respect to the fair weather value at different altitudes over time, (D) magnification of panel (C) in the region z < 5 km, (E) average number of elementary charges per particle as a function of particle radius for different average bin distribution altitudes, (F) ratio between average electrical force magnitude and gravity magnitude as a function of particle radius for different average bin distribution altitudes, (G) average particle velocity as a function of particle radius for different average bin distribution altitudes, (H) residual time of particles in the atmosphere until 99.9% of their concentration is lost, and relative difference with the results of N I p,i . In panels (E)-(G) with dashed line of the same color and marker are the same results in the case of N I p,i for comparison.
Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 709890 weather ion density ( Figure 10A). This increase is not unrealistic, as it represents the maximum measured values by Ryder et al. (2013b) during the Fennec 2011 campaign. Such an increase of the ion density leads to a conductivity reduction by a factor up to ∼0.3 ( Figure 10B) inside the layer, and an electric field magnitude enhancement by a factor of 1.7 ( Figure 10C). The particle competition in this case reduces the average particle charge by a factor up to 0.1-0.2 in the case of particles with radius larger than 50 μm, and up to 0.4 for the rest sizes ( Figure 10D), compared to the case of Figure 8. Similar reduction is observed in the case of the electrical force ( Figure 10E). Another interesting point that is revealed by Figure 10E, that for some bins the electrical force changes polarity and becomes positive. This means that the particles are charged positively. In Figure 10F, where a snapshot of the charge density dynamics is depicted, it is clear that the charge polarity in the case of particles with radius larger than 50 μm changes to positive values (Bin 14 -Bin 18 ). This happens for the following reason. Since the ion attachment rate is proportional to the particle number density (Eqs 30,31,38,39), the negative ions attach more to bins with smaller radii (that have several order of magnitude more particles than bins with larger radii). In the case of the distribution with N I p,i , as the particles move to lower altitudes, since the particle density is much less than the ambient ion density, the ions replenish from the ionization process and there are always enough negative ions to attach to particles (the ratio of ion densities is less than the mobility ratio, see Section 3.1). In the case 10N I p,i case the particle density is comparable to the ambient ion density, and therefore the negative ions do not have time to replenish. Given that most negative ions attach to small particles, the ratio of ion densities becomes larger than the mobilitiy ratio for larger particles, and they can become positively charged. Consequently, ion attachment alone is capable to charge positively large particles and negatively small ones, as was previously stated only for the triboelectric effect (e.g., Zhao et al., 2002;Lacks et al., 2008;Shinbrot and Herrmann, 2008;Merrison, 2012). Similar conclusions were derived by Yair and Levin (1989) based on the ion attachment to aerosols formulation by Hoppel and Frick (1986). They studied polydisperse spherical conducting aerosols with radii up to 0.5 μm and found that for small aerosol concentrations, particles were charged mostly negatively, but that an increase in aerosol concentrations lead to differential charging, resulting in opposite charging of the small (mostly negative) and the large (mostly positive) particles. This effect was enhanced for size distributions which had a significant component of large particles. Their results showed dependency only on the particle concentration, as opposed to the current work findings, that show dependency on the particle concentrations with respect to the atmospheric ion concentrations. The latter distinction is attributed to the difference in the scale of the problem under examination. In the present study, the particle sizes are larger than the ionic mean free path and therefore the ionic diffusion is expected to influence the attachment process. In the case of Yair and Levin (1989) the aerosol sizes are comparable, or smaller, than the ionic mean free path and the processes of three body trapping and image capture are dominant. In general, it becomes clear that regardless the scale of the problem, the ion attachment to FIGURE 10 | Comparison between the case of N I p,i and the case of 10N I p,i : (A) Total particle number density compared to the ambient fair weather ion density, (B) conductivity reduction with respect to the fair weather value at different altitudes over time for 10N I p,i , (C) electric field magnitude modification with respect to the fair weather value at different altitudes over time for 10N I p,i , (D) average particle charge ratio between N I p,i and 10N I p,i , for several bins, (E) average electrical force ratio between N I p,i and 10N I p,i , for several bins, (F) snapshot of particle charge density distribution for some bins in the case of 10N Frontiers in Earth Science | www.frontiersin.org August 2021 | Volume 9 | Article 709890 particles can lead to different polarities of the particle charges in case of different particle sizes. Figure 11 shows the case of the densely packed N II p,i particle distribution, enhanced by a factor of 10. In the specific case, the particle number density becomes larger than the fair weather ion densities ( Figure 11A). The conductivity inside the layer reduces by a factor up to ∼0.03 ( Figure 11B), while the electric field is enhanced by a factor up to ∼20 with respect to the fair weather values ( Figure 11C) The change of charge polarity in the case of large particles is more apparent in this case Figures 11D-F.
From the results of Figures 10, 11 we deduce that, the increase of the dust particle density favors the increase of the ambient large scale electric field, but does not favor the charging process of the dust particles, because their net charge decreases. Overall, the increase of the dust particle density acts against the effects of the electrical force, because the reduction of the particle charge is more prominent than the increase of the electric field, and the electrical force ultimately decreases. This indicates that the ion attachment process (acting on the enhancement of the ambient electric field) combined with a second charging process, such as the triboelectric effect (acting on the enhancement of the particle charge), might create a sufficient electrical force that in turn influences the settling of the dust particles. We plan on investigating this further on future work. Figure 12 illustrates the results regarding the comparison of different types of binning. Figure 12A shows the number of elementary charges per particle as a function of the particle radius. As the number of equilateral bins increases, the results converge to the results calculated by the current work binning. This becomes more clear in Figures 12C,E, where initial size range is split in sub ranges 10-150 μm, and 0.1-10 μm, respectively. The same conclusions are derived in the case of the ratio between the electrical force magnitude and the gravity magnitude (Figures 2, 12B,D,F).

CONCLUSION
Ion attachment to dust particles is one of the major charging mechanisms during their gravitational settling within the atmospheric column. It can, therefore, lead to an enhancement of the ambient fair weather electric field, due to the ion density and, consequently, the atmospheric electrical conductivity reduction. This is the first consistent attempt, to our knowledge, to analytically represent the phenomena by integrating a measured size distribution of spherical dust particles, in the 1D model that incorporates the calculation of new attachment rates of non stationary spherical particles.
We have shown that ion attachment is capable of charging dust particles, with radius between 1 −100 μm, with a net charge that ranges between 1 −1,000 elementary charges depending on the dust particle number density. The particle velocity enhances the acquired particle net charge, but converges to a steady state value as the particle is traveling in the fair weather atmosphere away from the initial dust FIGURE 11 | Comparison between the case of N II p, i and the case of 10N II p,i : (A) Total particle number density compared to the ambient fair weather ion density, (B) conductivity reduction with respect to the fair weather value at different altitudes over time for 10N II p,i , (C) electric field magnitude modification with respect to the fair weather value at different altitudes over time for 10N II p,i , (D) average particle charge ratio between N II p,i and 10N II p,i , for several bins, (E) average electrical force ratio between N II p,i and 10N II p,i , for several bins, (F) magnification of panel (E) for average altitude range between 4 and 5 km.
layer. The influence of the total number density on the particle charging process appears stronger than its influence on the electric field enhancement. If the number density is very close or becomes even larger than the ambient ion density, large particles can be charged positively while within the initial dust layer. This is attributed to the lack of negative ions, which tend to attach to small particles, that have number densities several orders of magnitude greater than the larger ones. A similar charging behavior was found in atmospheric aerosols with sizes smaller than 0.5 μm in radius, by Yair and Levin (1989). The dust layer depth is expected to influence both the net particle charge and the ambient electric field enhancement, but overall, the electrical force applied to the particles remains the same. The magnitude of the electrical force is found to be, at best, two orders of magnitude lesser than gravity, which leads to no significant influence to the particle dynamics, for particles with radius in the range 1−100 μm. Nonetheless, other processes must be combined with the ion attachment, such as the triboelectric effect or the updrafts, for a better quantification of the electrical force possible influence on the particle gravitational settling. In parallel to this work, sophisticated profiling measurements with novel atmospheric electricity sensors are launched in planned experiments within the D-TECT project, in order to further evaluate our model. These steps will be concluded in future work. FIGURE 12 | Study on the proper discretization of the size distribution: (A) Number of elementary charges per particle as a function of particle radius for different types of binnings, (B) ratio between the electrical force magnitude and the gravity magnitude as a function of particle radius for different types of binnings, (C,D) magnification of panels (A) and (B) for particle radius in the range of 10 − 150 μm, (E,F) magnification of panels (A) and (B) for particle radius in the range of 0.1 − 10 μm.