AerSett v1.0: a simple and straightforward model for the settling speed of big spherical atmospheric aerosols

This study introduces AerSett v1.0 (AERosol SETTling version 1.0), a model giving the settling speed of big spherical aerosols in the atmosphere without going through an iterative equation resolution. We prove that, for all spherical atmospheric aerosols with diameter D up to 1000 µm, this direct and explicit method including the drag coefficient formulation of Clift and Gauvin (1971) and the Davies (1945) slip correction factor gives results within 2 % of the exact solution obtained from the numerical resolution of a non-linear fixed-point equation. This error is acceptable considering the uncertainties on the drag coefficient formulations themselves. For D<100 µm, the error is below 0.5 %. We provide a Fortran implementation of this simple and straightforward model, hoping that more chemistry–transport models (CTMs) and general circulation models will be able to take into account large-particle drag correction to the settling speed of big spherical aerosol particles in the atmosphere, without performing an iterative and time-consuming calculation.

1 Introduction

One of the main purposes of chemistry–transport models (CTMs) is to simulate the aerosol concentration in the atmosphere as accurately as possible. The settling velocity of aerosols is a key driver of their dry removal from the atmosphere (Zhang et al., 2001). With dry removal being the only sink for atmospheric aerosol under dry conditions, any error on representing the settling velocity of atmospheric aerosol will have a direct impact on their modeled concentrations. Fortunately, for particles with diameter D<10 µm, which are the most significant for health effects, the Stokes law (Stokes, 1851) along with slip correction factors (Cunningham, 1910; Davies, 1945) gives a straightforward and accurate way to calculate the settling speed of aerosol particles.

However, dust particles exist in the atmosphere with different sizes from D≃0.1 µm to D>100 µm (Ryder et al., 2019). These authors classify mineral dust particles according to three modes as follows:

• $\mathrm{0.1}\mathrm{\mu }\mathrm{m}<D<\mathrm{2.5}\mathrm{\mu }\mathrm{m}$: accumulation mode,

• $\mathrm{2.5}\mathrm{\mu }\mathrm{m}<D<\mathrm{20}\mathrm{\mu }\mathrm{m}$: coarse mode,

• $\mathrm{20}\mathrm{\mu }\mathrm{m}<D<\mathrm{\infty }$: giant mode.

The contribution of the giant mode is substantial, at least over the Sahara: Ryder et al. (2019) show that not considering giant dust particles over the Sahara results in an underestimation of mass concentration by 40 % as well as extinction by as much as 18 % for shortwave radiation and 26 % for longwave radiation. Dust particles with diameter D up to 100 µm are present not only above the Sahara (Ryder et al., 2019) but also far away from emission sources: Betzer et al. (1988) have observed dust particles with D>75 µm in the atmosphere more than 10 000 km away from their source. More recently, van der Does et al. (2018) have observed dust particles with diameter up to 450 µm over the Atlantic Ocean, more than 2400 km away from the West African coast. Modeling coarse and giant dust particles is still a very challenging task. For example, Drakaki et al. (2022) show that the WRFV4.2.1 model with a version of the GOCART–AFWA dust scheme (see LeGrand et al., 2019), modified to include coarse and giant dust particles, underestimates the lifetime of coarse and giant dust particles in the atmosphere. They show that their simulation results are closer to observation when they include an artificial reduction of particles' settling velocities by 60 % to 80 % (depending on the diameter). This reduction is a way to account for underrepresented mechanisms such as non-sphericity of particles (Mallios et al., 2020), or their electric charges, which have been discussed as possible factors explaining a longer atmospheric lifetime of coarse dust particles (Adebiyi and Kok, 2020). Several observational and modeling studies have addressed the question of coarse and giant dust particles in the atmosphere: doing an exhaustive bibliographical overview on this question falls beyond the scope of the present study. For a more complete bibliography, the reader is referred to van der Does et al. (2018), Ryder et al. (2019) and Drakaki et al. (2022).

The settling speed of giant particles deviates substantially from the Stokes law, an effect that can be taken into account using mathematical formulations known as large-particle drag corrections. Usually, these large-particle drag corrections are performed by using empirical formulations of the drag coefficient C_d as a function of the Reynolds number Re (typically the one provided by Clift and Gauvin, 1971), and numerically solving an equation to obtain an estimate of the settling speed v_∞ as a function of the characteristics of the particle and of ambient air. This method is robust and permits studies like Drakaki et al. (2022) while taking this effect into account, but it induces large calculation costs. Since the important impact of giant dust particles on the dust concentration and optical effect has been demonstrated (e.g., Ryder et al., 2019), there is an emerging need to solve the problems that hinder the representation of giant dust particles in CTMs and general circulation models. Designing a robust and efficient method to calculate the settling speed of giant dust particles is a step in this direction. Until now, the gravitational settling speed in most chemistry–transport models is calculated with a plain Stokes formulation and a slip correction factor for the smallest particles, as in Zhang et al. (2001) but without a large-particle drag correction (e.g., Sič et al., 2015; Rémy et al., 2019; Shu et al., 2021). An exception to this is the recent development exposed by Drakaki et al. (2022) in the GOCART–AFWA dust scheme of WRFV4.2.1. In that study, the Clift and Gauvin (1971) drag coefficient correction is taken into account by a bisection method, performed at each time step, in each model cell and for each model size bin to calculate the settling speed as a function of the particle properties and the atmospheric conditions.

Since it has been highlighted by Drakaki et al. (2022) that taking large-particle drag correction into account is important for representing the settling of giant dust particles, the goal of this article is to give a simple, robust and computationally efficient expression to calculate the settling speed of spherical atmospheric aerosol, including large-particle drag correction. As discussed thoroughly in Goossens (2019), several parameterizations exist for the drag coefficient, each of them fitting the reference data of Lapple and Sheperd (1940) only in part. These parameterizations give drag coefficients which differ between them and from measurements by a few percents. Among these parameterizations, the Clift and Gauvin (1971) and Cheng (2009) formulations seem to perform better according to the objective scores presented in Table 2 of Goossens (2019). In the present study, we base our calculations on the Clift and Gauvin (1971) drag coefficient formulation and use Cheng (2009) as a benchmark of the uncertainty that can exist between several state-of-the-art drag coefficient formulations.

In Sect. 2, we review the basic equations that govern the settling speed of spherical aerosol particles in the atmosphere (leaving aside slip correction); in Sect. 3, we lay the basis for a new methodology for finding the settling speed of the particle as a function of the particle and surrounding gas parameters. In Sect. 4, we apply this method to the case of the drag coefficient parameterization of Clift and Gauvin (1971) and give an approximated expression of the large-particle drag correction factor. Sections 2, 3 and 4 focus on the large-particle drag correction, ignoring the slip correction factor. In Sect. 5, we extend the results of the previous sections by showing how to combine the slip correction factor with the large-particle drag correction, and we present the equations that define the AerSett (AERosol SETTling) method, including both effects. In Sect. 6, we evaluate the computational benefit of this method compared to other numerical resolution strategies. Finally, we conclude in Sect. 7 by summarizing the method we propose for the calculation of the settling speed of large spherical aerosol particles as well as future prospects for improving the representation of the settling speed of aerosol particles in chemistry–transport models.

2 Basic equations

We will follow the notations of Mallios et al. (2020). The settling speed v_∞ of a spherical particle with diameter D is given by

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{\mathrm{d}}{A}_{\mathrm{p}}{\mathit{\rho }}_{\mathrm{a}}{v}_{\mathrm{\infty }}^{\mathrm{2}}=\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}}\right)V_{\mathrm{p}}g,\end{array}$$

where C_d is the drag coefficient, A_p the projected area of the particle normal to the flow, ρ_a the density of air, ρ_p the density of the settling particle, v_∞ the settling speed of the particle, V_p the particle volume and g the acceleration of gravity.

In the case of spherical particles and for extremely small Reynolds numbers, Stokes (1851) has shown that

$$\begin{array}{}\text{(2)}& C_{\mathrm{d}}={\displaystyle \frac{\mathrm{12}}{\mathit{Re}}},\end{array}$$

where the Reynolds number Re is defined as follows:

$$\begin{array}{}\text{(3)}& \mathit{Re}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}D{v}_{\mathrm{\infty }}}{\mathrm{2}\mathit{\mu }}}.\end{array}$$

For all the calculations below, we have used the following empirical expression for dynamic viscosity μ (NOAA/NASA/USAF, 1976):

$$\begin{array}{}\text{(4)}& \mathit{\mu }={\displaystyle \frac{\mathit{\beta }{T}^{\frac{\mathrm{3}}{\mathrm{2}}}}{T+S}},\end{array}$$

where $\mathit{\beta }=\mathrm{1.458}\times {\mathrm{10}}^{-\mathrm{6}}\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\frac{\mathrm{1}}{\mathrm{2}}}$ and S=110.4 K.

While Eq. (2) only holds for extremely small Reynolds number (Re<0.1), Clift and Gauvin (1971) has given an empirical expression of C_d as a function of Re, which holds up to Re=10⁵:

$$\begin{array}{}\text{(5)}& C_{\mathrm{d}}={\displaystyle \frac{\mathrm{12}}{\mathit{Re}}}\left(\mathrm{1}+\mathrm{0.2415}{\mathit{Re}}^{\mathrm{0.687}}\right)+{\displaystyle \frac{\mathrm{0.42}}{\mathrm{1}+\frac{\mathrm{19}\mathrm{019}}{{\mathit{Re}}^{\mathrm{1.16}}}}},\end{array}$$

where exponents 0.687 and 1.16 and coefficients 0.2415, 0.42 and 19 019 are empirical values that permit Eq. (5) to fit experimental data. Hereinafter, we will refer to Eq. (5) as the Clift–Gauvin formula.

Equations (1) and (5) are two equations for two unknowns, v_∞ and C_d. For small Reynolds numbers (Re<0.1), Eq. (5) is reduced to Eq. (2) and we obtain:

$$\begin{array}{}\text{(6)}& v_{\mathrm{\infty }}^{\mathrm{Stokes}}={\displaystyle \frac{D^{\mathrm{2}}\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}}\right)g}{\mathrm{18}\mathit{\mu }}}.\end{array}$$

If the Reynolds number exceeds 0.1, Eq. (2) does not hold, and v_∞ does not have a general analytical expression. The solution to Eq. (2) can be obtained by an iterative fixed-point method as suggested in van Boxel (1998), or by a bisection method as in Drakaki et al. (2022). The results of this numerical calculation are shown in Fig. 1. Figure 1a shows that the Stokes equation (Eq. 6) gives excellent results for D<20 µm, but it also shows that deviations from it due to the departure of the drag coefficient from Eq. (2) gradually arise when D exceeds 20 µm, reaching −30 % when D≃100 µm and −90 % when D≃1000 µm (Fig. 1c). While particles with diameter D≃1000 µm are not a concern for chemistry–transport modeling, those with D≃100 µm are an emerging concern due to recent observations of particles with such diameters far away from their source (van der Does et al., 2018).

Solving Eq. (1) with an iterative method as suggested in van Boxel (1998) demands several iterations when the diameter of the particle gets close to 100 µm or beyond (see Sect. 6). This is why we will expose a way to estimate v_∞ from the physical parameters of the problem in a straightforward way.

Figure 1 (a) Settling speed of a spherical particle where $\mathit{\rho }=\mathrm{2650}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ in air with P=101 325 Pa and T=293.15 K as a function of diameter where C_d from Eqs. (5) (blue) and (6) (orange). (b) Reynolds number from Eqs. (5) (blue) and (6) (orange). (c) Relative difference of v_∞ from $v_{\mathrm{\infty }}^{\mathrm{Stokes}}$.

3 Expressing v_∞ from the parameters of the problem

To slightly generalize matters, let us rewrite Eq. (5) as

$$\begin{array}{}\text{(7)}& C_{\mathrm{d}}={\displaystyle \frac{\mathrm{12}}{\mathit{Re}}}\left(\mathrm{1}+f\left(\mathit{Re}\right)\right),\end{array}$$

where f is a function characterizing the deviation of C_d from its Stokes (1851) expression.

Injecting Eq. (7) into Eq. (1), with $A_{\mathrm{p}}=\frac{\mathit{\pi }{D}^{\mathrm{2}}}{\mathrm{4}}$ and $V_{\mathrm{p}}=\frac{\mathit{\pi }{D}^{\mathrm{3}}}{\mathrm{6}}$, we obtain:

$$\begin{array}{}\text{(8)}& {\displaystyle }{v}_{\mathrm{\infty }}& {\displaystyle }={\displaystyle \frac{D^{\mathrm{2}}\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}}\right)g}{\mathrm{18}\mathit{\mu }}}\times {\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}+f\left(\mathit{Re}\right)\right)}},\text{(9)}& {\displaystyle }& {\displaystyle }=v_{\mathrm{\infty }}^{\mathrm{Stokes}}\times {\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}+f\left(\mathit{Re}\right)\right)}}.\end{array}$$

Of course, Eq. (9) does not give an explicit formulation of v_∞ since Re depends on v_∞ (Eq. 3). However, we are looking for a way to take advantage of Eq. (9) to obtain an explicit estimate of v_∞ so that we introduce the deviation of v_∞ from $v_{\mathrm{\infty }}^{\mathrm{Stokes}}$ and δ as

$$\begin{array}{}\text{(10)}& v_{\mathrm{\infty }}=v_{\mathrm{\infty }}^{\mathrm{Stokes}}\left(\mathrm{1}+\mathit{\delta }\right).\end{array}$$

The terminal Reynolds number of the particle is equal to

$$\begin{array}{}\text{(11)}& \mathit{Re}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}D{v}_{\mathrm{\infty }}^{\mathrm{Stokes}}\left(\mathrm{1}+\mathit{\delta }\right)}{\mathrm{2}\mathit{\mu }}}=\left(\mathrm{1}+\mathit{\delta }\right){\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}{D}^{\mathrm{3}}\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}}\right)g}{\mathrm{36}{\mathit{\mu }}^{\mathrm{2}}}}.\end{array}$$

Let us introduce

$$\begin{array}{}\text{(12)}& R={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}{D}^{\mathrm{3}}\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}}\right)g}{\mathrm{36}{\mathit{\mu }}^{\mathrm{2}}}},\end{array}$$

where R is the Reynolds number of the particle if it were to fall at speed $v_{\mathrm{\infty }}^{\mathrm{Stokes}}$: we will call it the virtual Reynolds number. The Reynolds number of the particle falling at its corrected settling speed v_∞ is R(1+δ). With these modifications, Eq. (9) can be rewritten as follows:

$$\begin{array}{}\text{(13)}& \mathrm{1}+\mathit{\delta }={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+f\left(\left(\mathrm{1}+\mathit{\delta }\right)R\right)}}.\end{array}$$

In this equation, R is a non-dimensional number that depends on the characteristics of the problem (Eq. 12), and f is the relative deviation of the drag coefficient C_d from $\frac{\mathrm{12}}{\mathit{Re}}$ (Eq. 7). Independent of the formulation of f(Re), the relative deviation δ of v_∞ from $v_{\mathrm{\infty }}^{\mathrm{Stokes}}$ is the solution of the fixed-point Eq. (13). Therefore, δ is a function of R, and only of R, which we will note δ(R). In other words, the settling speed of particles can be expressed as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }{v}_{\mathrm{\infty }}& {\displaystyle }={\displaystyle \frac{D^{\mathrm{2}}\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}}\right)g}{\mathrm{18}\mathit{\mu }}}\times \left(\mathrm{1}+\mathit{\delta }\left(R\right)\right),\mathrm{with}\text{(15)}& {\displaystyle }R& {\displaystyle }={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}{D}^{\mathrm{3}}\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}}\right)g}{\mathrm{36}{\mathit{\mu }}^{\mathrm{2}}}}.\end{array}$$

This shows that the dependance of v_∞ on parameters D, ρ_p, ρ_a, g, μ and on function f(Re) has a very specific form, and that, if one wants to tabulate the solutions, the values of δ(R) could be obtained once and for all solving Eq. (13) for the entire relevant range of the possible values of R (where the virtual Reynolds number R is defined by Eq. 12), instead of tabulating v_∞ for all the possible combinations of D, ρ_p, ρ_a, g and μ. Instead of performing such a tabulation, in the next section, we will show that it is possible to find an analytical expression approximating δ(R).

4 The δ(R) function in the Clift–Gauvin case

Now, we proceed supposing that the expression of C_d as a function of Re is that of Clift and Gauvin (1971) (Eq. 5), yielding the following expression for f(Re):

$$\begin{array}{}\text{(16)}& f\left(\mathit{Re}\right)=\mathrm{0.2415}{\mathit{Re}}^{\mathrm{0.687}}+{\displaystyle \frac{\mathit{Re}}{\mathrm{12}}}\times {\displaystyle \frac{\mathrm{0.42}}{\mathrm{1}+\frac{\mathrm{19}\mathrm{019}}{{\mathit{Re}}^{\mathrm{1.16}}}}}.\end{array}$$

Equation (13) can then be solved iteratively, as suggested in van Boxel (1998). The resulting function δ(R) is shown in Fig. 2a. Due to the sigmoid shape of δ(R), fitting it with a logistic function of ln R is tempting and gives relatively good results. However, a generalized logistic function gives an even better agreement with the exact solution (Fig. 2a). The equation obtained with this fit is

$$\begin{array}{}\text{(17)}& {\displaystyle }\mathit{\delta }\left(R\right)& {\displaystyle }\simeq -{\left(\mathrm{1}+{\mathrm{e}}^{-\mathrm{0.4335}\left(\ln R-\mathrm{0.8921}\right)}\right)}^{-\mathrm{1.905}},\text{(18)}& {\displaystyle }& {\displaystyle }\simeq -{\left(\mathrm{1}+{\left({\displaystyle \frac{R}{\mathrm{2.440}}}\right)}^{-\mathrm{0.4335}}\right)}^{\mathrm{-}\mathrm{1.905}}.\end{array}$$

This expression of δ(R) yields an error relative to the exact solution <1 % up to Re=10 and <2.5 % up to Re=10³ (Fig. 2b). Considering that Goossens (2019) indicate that the Clift–Gauvin empirical formulation of the drag coefficient (Eq. 5) is true within 7 % of the reference drag coefficient values of Lapple and Sheperd (1940), approximating δ(R) by the generalized logistic function given in Eq. (18) is accurate enough not to degrade the evaluation of the settling speed, at least until Re=10³, which is well beyond the typical range of the Reynolds number for atmospheric aerosol in free fall (Fig. 1b). Figure 2c indeed shows that the error committed by applying the fit formula (Eq. 18) instead of actually solving Eq. (13) is smaller than the discrepancy between the Clift and Gauvin (1971) formulation and the Cheng (2009) parameterization. With the Clift and Gauvin (1971) and Cheng (2009) drag formulations being the two best-performing formulations according to the objective scores presented in Goossens (2019) (their Table 2), this confirms that the error introduced by using Eq. (18) instead of the exact solution of Eq. (13) is smaller than the uncertainty of state-of-the-art drag coefficient formulations.

Figure 2 (a) δ(R) (solution of Eq. 13) as a function of R when f is defined from the Clift–Gauvin expression (Eq. 16). The red line is a fit of the solution by a generalized logistic function (Eq. 18). (b) Percent error of the fitted expression of v_∞ relative to the exact solution. (c) Percent difference of the settling speed v_∞ calculated with the Clift and Gauvin (1971) parameterization relative to the one calculated with the Cheng (2009) parameterization.

Figure 3a shows that, for all realistic particle sizes and all realistic atmospheric pressures in the troposphere, the Reynolds number is below 100, and Fig. 3b shows that the error induced by using Eq. (18) to evaluate the solution of Eq. (13) is less than 0.5 % for all particles with D<100 µm, and less than 2 % for all particles with diameter less than D<1000 µm: this shows that the domain of applicability of Eq. (18) largely covers the size range of giant dust particles (for which the typical diameter is below or around 100 µm, with exceptional observations of particles with D≃400 µm as in van der Does et al., 2018).

Figure 3 (a) v_∞ from the Clift–Gauvin expression (Eq. 16) as a function of atmospheric pressure P and particle diameter D for an atmospheric particle with density $\mathit{\rho }=\mathrm{2650}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. The dependance of temperature on pressure has been taken from the US standard atmosphere (NOAA/NASA/USAF, 1976). Black contours are iso-Re contours. (b) Percent error committed by using Eq. (18) instead of solving Eq. (13).

5 Inclusion of the slip correction factor

So far, for the sake of simplicity, we have assumed that continuum fluid mechanics apply to our falling particles. As explained in, e.g., Seinfeld and Pandis (1997) (Chap. 8), this assumption holds only if Kn≪1, where Kn is the Knudsen number of the falling particle:

$$\begin{array}{}\text{(19)}& \mathit{Kn}={\displaystyle \frac{\mathrm{2}\mathit{\lambda }}{D}},\end{array}$$

where λ is the mean free path of molecules in air. The mean free path in air is given by the following empirical equation as a function of pressure P, dynamic viscosity μ and density ρ_a (Jennings, 1988):

$$\begin{array}{}\text{(20)}& \mathit{\lambda }=\sqrt{{\displaystyle \frac{\mathit{\pi }}{\mathrm{8}}}}\times {\displaystyle \frac{\frac{\mathit{\mu }}{\mathrm{0.4987445}}}{\sqrt{P{\mathit{\rho }}_{\mathrm{a}}}}}.\end{array}$$

From the Knudsen number, slip correction factors can be designed to account for the non-continuous effects, turning Eq. (1) into

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{C_{\mathrm{d}}}{C_{\mathrm{c}}}}{A}_{\mathrm{p}}{\mathit{\rho }}_{\mathrm{a}}{v}_{\mathrm{\infty }}^{\mathrm{2}}=\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}}\right)V_{\mathrm{p}}g,\end{array}$$

where the slip correction factor C_c (also known as the Cunningham (1910) correction factor) is usually estimated from the Davies (1945) expression:

$$\begin{array}{}\text{(22)}& C_{\mathrm{c}}=\mathrm{1}+\mathit{Kn}\left(\mathrm{1.257}+\mathrm{0.4}\exp \left(-{\displaystyle \frac{\mathrm{1.1}}{\mathit{Kn}}}\right)\right).\end{array}$$

From Eqs. (7) and (21), and if we introduce the Stokes terminal velocity including the slip correction term, ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}$, as

$$\begin{array}{}\text{(23)}& {\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}=C_{\mathrm{c}}{v}_{\mathrm{\infty }}^{\mathrm{Stokes}},\end{array}$$

we obtain the following:

$$\begin{array}{}\text{(24)}& v_{\mathrm{\infty }}={\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}\times {\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}+f\left(\mathit{Re}\right)\right)}}.\end{array}$$

The terminal Reynolds number of the particle is equal to

$$\begin{array}{}\text{(25)}& \mathit{Re}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}D{\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}\left(\mathrm{1}+\mathit{\delta }\right)}{\mathrm{2}\mathit{\mu }}}=\left(\mathrm{1}+\mathit{\delta }\right){\displaystyle \frac{C_{\mathrm{c}}{\mathit{\rho }}_{\mathrm{a}}{D}^{\mathrm{3}}\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}}\right)g}{\mathrm{36}{\mathit{\mu }}^{\mathrm{2}}}}\end{array}$$

Let us introduce:

$$\begin{array}{}\text{(26)}& \stackrel{\mathrm{̃}}{R}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}D{\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}}{\mathrm{2}\mathit{\mu }}}.\end{array}$$

With these notations, $\mathit{Re}=\left(\mathrm{1}+\mathit{\delta }\right)\stackrel{\mathrm{̃}}{R}$, and Eq. (24) becomes

$$\begin{array}{}\text{(27)}& \left(\mathrm{1}+\mathit{\delta }\right)={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}+f\left(\left(\mathrm{1}+\mathit{\delta }\right)\stackrel{\mathrm{̃}}{R}\right)\right)}}.\end{array}$$

Equation (27) is the same fixed-point equation as Eq. (13) but with parameter $\stackrel{\mathrm{̃}}{R}$ (defined from Eq. 26) instead of R (from Eq. 12). Therefore, its solution is $\mathit{\delta }\left(\stackrel{\mathrm{̃}}{R}\right)$, where the δ function is the same as represented in Fig. 2a, and it can correspondingly be approximated by Eq. (18), with the same error term as represented in Fig. 2b. The fact that the introduction of the slip correction factor C_c changes the mathematical method only very slightly to obtain the expression of the settling speed v_∞ is a fortunate consequence of the slip correction factor being a function of the Knudsen number only, with no dependance on particle speed.

To summarize the previous development, the method we have designed to calculate v_∞ including the slip correction term and the drag correction term takes the following steps:

1. Calculate Kn from Eq. (19) and C_c from Eq. (22).

2. Calculate ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}$ from Eq. (23) and $\stackrel{\mathrm{̃}}{R}$ from Eq. (26).

3. Calculate v_∞ as

   $v_{\mathrm{\infty }}={\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}\times \left[\mathrm{1}-{\left(\mathrm{1}+{\left(\frac{\stackrel{\mathrm{̃}}{R}}{\mathrm{2.440}}\right)}^{-\mathrm{0.4335}}\right)}^{\mathrm{-}\mathrm{1.905}}\right]$.

Figure 4 shows the impact of the slip correction term and the drag correction term on ${\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}$. Two regimes are clearly separated in this figure, with $\frac{v_{\mathrm{\infty }}}{v_{\mathrm{\infty }}^{\mathrm{Stokes}}}>\mathrm{1}$ for the smaller particles, for which the slip correction term dominates, and the larger particles for which the drag correction term dominates. Between these two regimes, a relatively large zone exists in which the departure of v_∞ from $v_{\mathrm{\infty }}^{\mathrm{Stokes}}$ is less than 5 %. This zone in which the Stokes equation is directly applicable covers the $\mathrm{3}\mathrm{\mu }\mathrm{m}<D<\mathrm{35}\mathrm{\mu }\mathrm{m}$ diameter range at ground-level pressure and the $\mathrm{10}\mathrm{\mu }\mathrm{m}<D<\mathrm{50}\mathrm{\mu }\mathrm{m}$ range at 200 hPa. Some authors (e.g., Mallios et al., 2020) argue that the slip correction factor C_c should only be applied in the Stokes regime (Re<0.1). However, in Fig. 4, it has been applied regardless of the Reynolds number. Analysis of Fig. 4 shows that the choice of applying the slip correction factor for Re>0.1 or not has little consequence, at least in the troposphere, because in all the portions of the pressure–diameter diagram where Re>0.1 (above the red line in Fig. 4), C_c is comprised between 1 and 1.02, so whether or not this factor which is extremely close to 1 is applied has no practical consequence on the simulated value of v_∞.

Figure 4 Ratio $\frac{v_{\mathrm{\infty }}}{v_{\mathrm{\infty }}^{\mathrm{Stokes}}}$ as a function of atmospheric pressure P and particle diameter D for a spherical particle with ${\mathit{\rho }}_{\mathrm{p}}=\mathrm{2650}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ (color shades). Black contours are iso-lines for the slip correction factor $C_{\mathrm{c}}=\frac{{\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}}{v_{\mathrm{\infty }}^{\mathrm{Stokes}}}$, and white lines are iso-lines of $\mathrm{1}+\mathit{\delta }\left(\stackrel{\mathrm{̃}}{R}\right)=\frac{v_{\mathrm{\infty }}}{{\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}}$. The red line corresponds to Re=0.1, where Re>0.1 above this line and Re<0.1 below.

6 Implementation and computational efficiency

Apart from its simplicity, another possible advantage of the straightforward evaluation of v_∞ using the three steps defined above is its computational efficiency compared to the currently available methods:

• Bisection method. Find Re as the solution to $\mathit{Re}=\frac{\stackrel{\mathrm{̃}}{R}}{\mathrm{1}+f\left(\mathit{Re}\right)}$ starting from $\mathit{Re}\in \left[\mathrm{0};\stackrel{\mathrm{̃}}{R}\right]$ and successively cut the solution interval [Re_m;Re_M] into halves until the relative error on speed becomes smaller than ε: $\frac{{\mathit{Re}}_{\mathrm{M}}-{\mathit{Re}}_{\mathrm{m}}}{{\mathit{Re}}_{\mathrm{m}}}<\mathit{\epsilon }$; then $v_{\mathrm{\infty }}={\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}\frac{\mathit{Re}}{\stackrel{\mathrm{̃}}{R}}$.

• Fixed-point method. Find the solution δ to Eq. (27) starting from δ₀=0 and iterate it with ${\mathit{\delta }}_{i+\mathrm{1}}=\frac{\mathrm{1}}{\mathrm{1}+f\left(\stackrel{\mathrm{̃}}{R}\times \left(\mathrm{1}+{\mathit{\delta }}_i\right)\right)}-\mathrm{1}$ until the relative error on speed becomes smaller than ε: $\frac{\left|{\mathit{\delta }}_{i+\mathrm{1}}-{\mathit{\delta }}_i\right|}{\mathrm{1}+{\mathit{\delta }}_{i+\mathrm{1}}}<\mathit{\epsilon }$. Then, $v_{\mathrm{\infty }}={\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}\left(\mathrm{1}+\mathit{\delta }\right)$.

• AerSett method. Directly calculate $v_{\mathrm{\infty }}={\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}\times \left[\mathrm{1}-{\left(\mathrm{1}+{\left(\frac{\stackrel{\mathrm{̃}}{R}}{\mathrm{2.440}}\right)}^{-\mathrm{0.4335}}\right)}^{\mathrm{-}\mathrm{1.905}}\right]$.

The bisection and fixed-point methods depend on a tolerance parameter ε, which we have set to ε=0.02, corresponding to the maximal error of 2 % observed using the AerSett method (Fig. 3b).

A Fortran code has been designed to estimate the calculation cost for these three methods (Fig. 5). This has been done by performing the calculation of v_∞ for 10⁸ random values of D, at a random altitude z in the atmosphere between 0 and 12 000 m. The pressure and temperature values P(z) and T(z) are estimated from the US Standard Atmospheric profile (NOAA/NASA/USAF, 1976).

To optimize calculation speed, we have observed that for $\stackrel{\mathrm{̃}}{R}<\mathrm{0.0116}$ we have $\mathrm{1}+\mathit{\delta }\left(\stackrel{\mathrm{̃}}{R}\right)>\mathrm{0.99}$, so that for all three methods we speed up the calculation by assuming that $v_{\mathrm{\infty }}={\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}$ when $\stackrel{\mathrm{̃}}{R}<\mathrm{0.0116}$. As can be seen in Fig. 4 (see the 0.99 iso-line of $\mathrm{1}+\mathit{\delta }\left(\stackrel{\mathrm{̃}}{R}\right)$), this means that for all particles with a diameter D<20 µm (the precise threshold value depends on pressure and temperature), no drag correction factor is applied.

Figure 5 shows that, for D<10 µm, the computation time (around 12 ns per call) is the same for the three methods, which is expected since the drag correction factor is not calculated for these diameters, resulting in $\stackrel{\mathrm{̃}}{R}<\mathrm{0.0116}$. In average over the entire interval, the bisection method and the fixed-point method are about 4 times slower than the AerSett method. For D>100 µm, the bisection method is 6.6 times slower than AerSett, and the fixed-point method is 8 times slower. It is also worth noting that in Fig. 5, the fixed-point method is faster than the bisection method for small values of D but slower for large values of D. The average number of iterations to obtain convergence increases from 4.5 ($\mathrm{10}\mathrm{\mu }\mathrm{m}<D<\mathrm{100}\mathrm{\mu }\mathrm{m}$) to 7.8 ($\mathrm{100}\mathrm{\mu }\mathrm{m}<D<\mathrm{1000}\mathrm{\mu }\mathrm{m}$) for the bisection method, and from 1.5 ($\mathrm{10}\mathrm{\mu }\mathrm{m}<D<\mathrm{100}\mathrm{\mu }\mathrm{m}$) to 7.4 ($\mathrm{100}\mathrm{\mu }\mathrm{m}<D<\mathrm{1000}\mathrm{\mu }\mathrm{m}$) for the fixed-point method. This sharper increase in the number of iterations to obtain convergence in the fixed-point method explains why bisection becomes more efficient for large values of D.

Figure 5 Execution time in nanoseconds (ns) for the three numerical methods described in Sect. 6 as a function of D. The estimate has been performed for 4 ranges of D (0.1–1, 1–10, 10–100, 100–1000 µm). In each of these intervals, the calculation of v_∞ as a function of D, ρ_p, ρ_a, λ and μ has been performed 10⁸ times with each method, for 10⁸ values of D selected randomly within the range. The execution time in nanoseconds corresponds to the total execution time (for 10⁸ calls to the function) divided by the number of calls (10⁸) so that it represents an estimate of the execution time for one call to the function. The test has been performed on a laptop with an Intel Core i7-1165G7 CPU.

7 Conclusions

As a conclusion, we have found that the following method is suitable to evaluate the settling speed of spherical aerosol particles in the atmosphere, v_∞:

$$\begin{array}{}\text{(28)}& {\displaystyle }{\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}=C_{\mathrm{c}}{\displaystyle \frac{D^{\mathrm{2}}\left({\mathit{\rho }}_{\mathrm{p}}-{\mathit{\rho }}_{\mathrm{a}}\right)g}{\mathrm{18}\mathit{\mu }}},\text{(29)}& {\displaystyle }\stackrel{\mathrm{̃}}{R}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}D{\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}}{\mathrm{2}\mathit{\mu }}},\text{(30)}& {\displaystyle }{v}_{\mathrm{\infty }}={\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}\times \left[\mathrm{1}-{\left(\mathrm{1}+{\left({\displaystyle \frac{\stackrel{\mathrm{̃}}{R}}{\mathrm{2.440}}}\right)}^{-\mathrm{0.4335}}\right)}^{\mathrm{-}\mathrm{1.905}}\right].\end{array}$$

The slip correction factor C_c can be obtained by the Davies (1945) formula (Eq. 22). To reduce computation time, for $\stackrel{\mathrm{̃}}{R}<\mathrm{0.116}$, Eq. (30) can be replaced by just $v_{\mathrm{\infty }}={\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}$, changing the result by less than 1 %.

Equations (28), (29) and (30) constitute the AerSett model v1.0 (AerSett for AERosol SETTling). The error induced by applying this model compared to an iterative calculation of v_∞ is less than 0.5 % for particles with diameter D<100 µm and less than 2 % for particles with D<1000 µm (Fig. 3b). Particles with larger diameters fall so rapidly that they are not relevant as atmospheric aerosol: other parameterizations exist for the falling hydrometeors (Khvorostyanov and Curry, 2005), taking into account the shape of snow flakes, the deformation of raindrops due to their speed, etc. We have shown that the error due to using Eq. (30) is smaller than the uncertainty that exists between different state-of-the-art formulations of the drag coefficient, showing that this error is not a problem for modeling.

Equation (30) takes into account both the slip correction factor (for small particles) and the large-particle drag correction factor. Figure 4 shows that for particles smaller than D≃10 µm, the slip correction factor dominates, while for particles larger than 20 µm, the large-particle drag correction factor dominates. The reduction in settling speed relative to $v_{\mathrm{\infty }}^{\mathrm{Stokes}}$ reaches about 25 % for a particle with diameter D≃100 µm, a typical diameter for giant dust particles. So, if chemistry–transport models are to represent the observed giant dust particles (which is still a challenge), large-particle drag correction for the gravitational settling speed needs to be taken into account, and we think that Eq. (30), valid for all spherical particles with D<1000 µm and at least from the surface to p=200 hPa (Figs. 3b and 4), is the simplest available method to do it. While designed specifically to take into account large-particle drag correction, Eq. (30) includes the slip correction term as well, which is critical for the finer atmospheric aerosol with D<10 µm. Using Eq. (30) yields an error smaller than 2 % relative to the exact combined effect of the Davies (1945) slip correction factor and the Clift and Gauvin (1971) large-particle drag correction factor. Therefore, it is possible to use this formulation systematically for spherical particles in chemistry–transport models. An implementation of Eqs. (28)–(30) in a Fortran module along with the necessary thermodynamical calculations is available for download and use (see “Code availability” section). We have shown (Sect. 6, and in particular Fig. 5) that the use of this Fortran routine allows us to gain a factor of 4 (in general) to 8 (for particles with D>100 µm) relative to bisection or fixed-point methods.

To go further in understanding the settling speed of giant dust particles and be able to represent them in chemistry–transport models, it is necessary to extend simple models such as AerSett to the case of non-spherical particles, and give simple and straightforward estimates of the drag correction factor δ, not only as a function of particle density and diameter, as in the present study, but also including other factors such as particle eccentricity which have been shown to have a strong impact on the settling speed of dust particles (Mallios et al., 2020, and references therein). To solve the persistent mystery of the processes allowing giant dust particles to stay airborne over long distances, new findings on physical processes such as the electric charges of the particles and their effect on settling velocities are still needed.