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

Sylvain Mailler, Laurent Menut, Arineh Cholakian, Romain Pennel

► To cite this version:

Sylvain Mailler, Laurent Menut, Arineh Cholakian, Romain Pennel. AerSett v1.0: a simple and straightforward model for the settling speed of big spherical atmospheric aerosols. Geoscientific Model Development, 2023, 16, pp.1119 - 1127. 10.5194/gmd-16-1119-2023. hal-04107228

HAL Id: hal-04107228

https://enpc.hal.science/hal-04107228

Submitted on 26 May 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Sylvain Mailler¹,², Laurent Menut¹, Arineh Cholakian¹, and Romain Pennel¹

¹LMD/IPSL, École Polytechnique, Institut Polytechnique de Paris, ENS, PSL Research University, Sorbonne Université, CNRS, Palaiseau, France

²École des Ponts ParisTech, Marne-la-Vallée, France

**Correspondence:** Sylvain Mailler (sylvain.mailler@lmd.ipsl.fr)

Received: 4 October 2022 – Discussion started: 21 October 2022
Revised: 13 January 2023 – Accepted: 2 February 2023 – Published: 15 February 2023

**Abstract.** 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:

- 0.1 µm < *D* < 2.5 µm: accumulation mode,

- 2.5 µm < *D* < 20 µm: coarse mode,

- 20 µm < *D* < ∞: 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 long-wave 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_\infty$ 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_\infty$ of a spherical particle with diameter $D$ is given by

$$ \frac{1}{2} C_d A_p \rho_a v_\infty^2 = (\rho_p - \rho_a) V_p g, \quad (1) $$

where $C_d$ is the drag coefficient, $A_p$ the projected area of the particle normal to the flow, $\rho_a$ the density of air, $\rho_p$ the density of the settling particle, $v_\infty$ 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

$$ C_d = \frac{12}{Re}, \quad (2) $$

where the Reynolds number $Re$ is defined as follows:

$$ Re = \frac{\rho_a D v_\infty}{2\mu}. \quad (3) $$
For all the calculations below, we have used the following empirical expression for dynamic viscosity $\mu$ (NOAA-/NASA/USAF, 1976):

$$ \mu = \frac{\beta T^{\frac{3}{2}}}{T + S}, \qquad (4) $$

where $\beta = 1.458 \times 10^{-6} \text{ kg s}^{-1} \text{ m}^{-1} \text{ K}^{-\frac{1}{2}}$ and $S = 110.4 \text{ 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^5$:

$$ C_d = \frac{12}{Re} \left(1 + 0.2415 Re^{0.687}\right) + \frac{0.42}{1 + \frac{19.019}{Re^{1.16}}}, \qquad (5) $$

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_\infty$ and $C_d$. For small Reynolds numbers ($Re < 0.1$), Eq. (5) is reduced to Eq. (2) and we obtain:

$$ v_{\infty}^{\text{Stokes}} = \frac{D^2 (\rho_{\text{p}} - \rho_{\text{a}}) g}{18 \mu}. \qquad (6) $$

If the Reynolds number exceeds 0.1, Eq. (2) does not hold, and $v_\infty$ 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 \mu\text{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 $\mu\text{m}$, reaching -30% when $D \approx 100 \mu\text{m}$ and -90% when $D \approx 1000 \mu\text{m}$ (Fig. 1c). While particles with diameter $D \approx 1000 \mu\text{m}$ are not a concern for chemistry-transport modeling, those with $D \approx 100 \mu\text{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 $\mu\text{m}$ or beyond (see Sect. 6). This is why we will expose a way to estimate $v_\infty$ from the physical parameters of the problem in a straightforward way.

### 3 Expressing $v_\infty$ from the parameters of the problem

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

$$ C_d = \frac{12}{Re} (1 + f (Re)), \qquad (7) $$

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

Injecting Eq. (7) into Eq. (1), with $A_p = \frac{\pi D^2}{4}$ and $V_p = \frac{\pi D^3}{6}$, we obtain:

$$ v_{\infty} = \frac{D^2 (\rho_{\text{p}} - \rho_{\text{a}}) g}{18 \mu} \times \frac{1}{(1 + f (Re))}, \qquad (8) $$

$$ = v_{\infty}^{\text{Stokes}} \times \frac{1}{(1 + f (Re))}. \qquad (9) $$

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

$$ v_{\infty} = v_{\infty}^{\text{Stokes}} (1 + \delta). \qquad (10) $$

The terminal Reynolds number of the particle is equal to

$$ Re = \frac{\rho_{\text{a}} D v_{\infty}^{\text{Stokes}} (1 + \delta)}{2\mu} = (1 + \delta) \frac{\rho_{\text{a}} D^3 (\rho_{\text{p}} - \rho_{\text{a}}) g}{36\mu^2}. \qquad (11) $$

Let us introduce

$$ R = \frac{\rho_{\text{a}} D^3 (\rho_{\text{p}} - \rho_{\text{a}}) g}{36\mu^2}, \qquad (12) $$

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

$$ 1 + \delta = \frac{1}{1 + f((1 + \delta) R)}. \qquad (13) $$

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{12}{Re}$ (Eq. 7). Independent of the formulation of $f(Re)$, the relative deviation $\delta$ of $v_\infty$ from $v_\infty^{\text{Stokes}}$ is the solution of the fixed-point Eq. (13). Therefore, $\delta$ is a function of $R$, and only of $R$, which we will note $\delta(R)$. In other words, the settling speed of particles can be expressed as follows:

$$ v_{\infty} = \frac{D^2 (\rho_{\text{p}} - \rho_{\text{a}}) g}{18 \mu} \times (1 + \delta(R)), \quad \text{with} \qquad (14) $$

$$ R = \frac{\rho_{\text{a}} D^3 (\rho_{\text{p}} - \rho_{\text{a}}) g}{36\mu^2}. \qquad (15) $$

This shows that the dependence of $v_\infty$ on parameters $D$, $\rho_{\text{p}}$, $\rho_{\text{a}}$, $g$, $\mu$ and on function $f(Re)$ has a very specific form, and that, if one wants to tabulate the solutions, the values of $\delta(R)$ could be obtained once and for all solving Eq. (13) for the entire relevant range of the possible values of $R$ (where the
<figure><img src="image_3.png" /><figcaption>Figure 1. (a) Settling speed of a spherical particle where $\rho = 2650 \text{ kg m}^{-3}$ in air with $P = 101325 \text{ Pa}$ and $T = 293.15 \text{ 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_\infty$ from $v_\infty^{\text{Stokes}}$.</figcaption></figure>

virtual Reynolds number $R$ is defined by Eq. 12), instead of tabulating $v_\infty$ for all the possible combinations of $D$, $\rho_p$, $\rho_a$, $g$ and $\mu$. Instead of performing such a tabulation, in the next section, we will show that it is possible to find an analytical expression approximating $\delta(R)$.

## 4 The $\delta(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)$:

$$f(Re) = 0.2415Re^{0.687} + \frac{Re}{12} \times \frac{0.42}{1 + \frac{19019}{Re^{1.16}}} \quad (16)$$

Equation (13) can then be solved iteratively, as suggested in van Boxel (1998). The resulting function $\delta(R)$ is shown in Fig. 2a. Due to the sigmoid shape of $\delta(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

$$\delta(R) \approx - \left(1 + e^{-0.4335(\ln R - 0.8921)}\right)^{-1.905}, \quad (17)$$

$$\approx - \left(1 + \left(\frac{R}{2.440}\right)^{-0.4335}\right)^{-1.905}. \quad (18)$$

This expression of $\delta(R)$ yields an error relative to the exact solution < 1% up to $Re = 10$ and < 2.5% up to $Re = 10^3$ (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 $\delta(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^3$, 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 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 \text{ µm}$, and less than 2% for all particles with diameter less than $D < 1000 \text{ µ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 \simeq 400 \text{ µm}$ as in van der Does et al., 2018).

## 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 \ll 1$, where $Kn$ is the Knudsen number of the falling particle:

$$Kn = \frac{2\lambda}{D}, \quad (19)$$

where $\lambda$ 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 $\mu$ and density
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) $\delta(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_\infty$ relative to the exact solution. (c) Percent difference of the settling speed $v_\infty$ calculated with the Clift and Gauvin (1971) parameterization relative to the one calculated with the Cheng (2009) parameterization.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a) $v_\infty$ from the Clift-Gauvin expression (Eq. 16) as a function of atmospheric pressure $P$ and particle diameter $D$ for an atmospheric particle with density $\rho = 2650 \text{ kg m}^{-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).</figcaption></figure>

$\rho_a$ (Jennings, 1988):

$$ \lambda = \sqrt{\frac{\pi}{8}} \times \frac{\frac{0.4987445}{\sqrt{P\rho_a}}}{\sqrt{P\rho_a}}. \qquad (20) $$

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

$$ \frac{1}{2} \frac{C_d}{C_c} A_p \rho_a v_\infty^2 = (\rho_p - \rho_a) V_p g, \qquad (21) $$

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

$$ C_c = 1 + Kn \left( 1.257 + 0.4 \exp \left( -\frac{1.1}{Kn} \right) \right). \qquad (22) $$

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

$$ \tilde{v}_\infty^{\text{Stokes}} = C_c v_\infty^{\text{Stokes}}, \qquad (23) $$

we obtain the following:

$$ v_\infty = \tilde{v}_\infty^{\text{Stokes}} \times \frac{1}{(1 + f(Re))}. \qquad (24) $$

The terminal Reynolds number of the particle is equal to

$$ Re = \frac{\rho_a D \tilde{v}_\infty^{\text{Stokes}} (1 + \delta)}{2\mu} = (1 + \delta) \frac{C_c \rho_a D^3 (\rho_p - \rho_a) g}{36\mu^2} \qquad (25) $$

Let us introduce:

$$ \tilde{R} = \frac{\rho_a D \tilde{v}_\infty^{\text{Stokes}}}{2\mu}. \qquad (26) $$

With these notations, $Re = (1 + \delta) \tilde{R}$, and Eq. (24) becomes

$$ (1 + \delta) = \frac{1}{(1 + f((1 + \delta) \tilde{R}))}. \qquad (27) $$

Equation (27) is the same fixed-point equation as Eq. (13) but with parameter $\tilde{R}$ (defined from Eq. 26) instead of $R$ (from Eq. 12). Therefore, its solution is $\delta(\tilde{R})$, where the $\delta$ 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_\infty$ is a fortunate consequence of the slip correction factor being a function of the Knudsen number only, with no dependence on particle speed.

To summarize the previous development, the method we have designed to calculate $v_\infty$ 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 $\tilde{v}_\infty^{\text{Stokes}}$ from Eq. (23) and $\tilde{R}$ from Eq. (26).

3. Calculate $v_\infty$ as

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

Figure 4 shows the impact of the slip correction term and the drag correction term on $\tilde{v}_\infty$. Two regimes are clearly separated in this figure, with $\frac{v_\infty}{v_\infty^{\text{Stokes}}} > 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_\infty$ from $v_\infty^{\text{Stokes}}$ is less than 5%. This zone in which the Stokes equation is directly applicable covers the 3 µm < $D$ < 35 µm diameter range at ground-level pressure and the 10 µm < $D$ < 50 µ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_\infty$.

## 6 Implementation and computational efficiency

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

- **Bisection method.** Find $Re$ as the solution to $Re = \frac{R}{1+f(Re)}$ starting from $Re \in [0; \tilde{R}]$ and successively cut the solution interval [$Re_M; Re_M$] into halves until the relative error on speed becomes smaller than $\varepsilon$: $\frac{Re_M - Re_M}{Re_M} < \varepsilon$; then $v_\infty = \tilde{v}_\infty^{\text{Stokes}} \frac{Re_\varepsilon}{R}$.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Ratio $\frac{v_\infty}{v_\infty^{\text{Stokes}}}$ as a function of atmospheric pressure $P$ and particle diameter $D$ for a spherical particle with $\rho_p = 2650 \text{ kg m}^{-3}$ (color shades). Black contours are iso-lines for the slip correction factor $C_c = \frac{\tilde{v}_\infty^{\text{Stokes}}}{v_\infty^{\text{Stokes}}}$, and white lines are iso-lines of $1 + \delta(\tilde{R}) = \frac{v_\infty}{\tilde{v}_\infty^{\text{Stokes}}}$. The red line corresponds to $Re = 0.1$, where $Re > 0.1$ above this line and $Re < 0.1$ below.</figcaption></figure>

- **Fixed-point method.** Find the solution $\delta$ to Eq. (27) starting from $\delta_0 = 0$ and iterate it with $\delta_{i+1} = \frac{1}{1+f(\tilde{R} \times (1+\delta_i))} - 1$ until the relative error on speed becomes smaller than $\varepsilon$: $\frac{|\delta_{i+1} - \delta_i|}{1+\delta_{i+1}} < \varepsilon$. Then, $v_\infty = \tilde{v}_\infty^{\text{Stokes}} (1+\delta)$.

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

The bisection and fixed-point methods depend on a tolerance parameter $\varepsilon$, which we have set to $\varepsilon = 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_\infty$ for $10^8$ 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 $\tilde{R} < 0.0116$ we have $1 + \delta(\tilde{R}) > 0.99$, so that for all three methods we speed up the calculation by assuming that $v_\infty = \tilde{v}_\infty^{\text{Stokes}}$ when $\tilde{R} < 0.0116$. As can be seen in Fig. 4 (see the 0.99 iso-line of $1 + \delta(\tilde{R})$), this means that for all particles with a diameter $D < 20 \mu\text{m}$ (the precise threshold value depends on pressure and temperature), no drag correction factor is applied.
<figure><img src="image_7.png" /><figcaption>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.</figcaption></figure>

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 $\tilde{R} < 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 (10 µm < D < 100 µm) to 7.8 (100 µm < D < 1000 µm) for the bisection method, and from 1.5 (10 µm < D < 100 µm) to 7.4 (100 µm < D < 1000 µ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.

## 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_∞$:

$$ \tilde{v}_{\infty}^{\text{Stokes}} = C_c \frac{D^2 (\rho_p - \rho_a) g}{18 \mu}, \quad (28) $$

$$ \tilde{R} = \frac{\rho_a D \tilde{v}_{\infty}^{\text{Stokes}}}{2\mu}, \quad (29) $$

$$ v_{\infty} = \tilde{v}_{\infty}^{\text{Stokes}} \times \left[ 1 - \left( 1 + \left( \frac{\tilde{R}}{2.440} \right)^{-0.4335} \right)^{-1.905} \right]. \quad (30) $$

The slip correction factor $C_c$ can be obtained by the Davies (1945) formula (Eq. 22). To reduce computation time, for $\tilde{R} < 0.116$, Eq. (30) can be replaced by just $v_∞ = \tilde{v}_{\infty}^{\text{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 \simeq 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_∞^{\text{Stokes}}$ reaches about 25% for a particle with diameter $D \simeq 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 $\delta$, 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.

**Code availability.** All the simulations and figures have been performed with Python scripts and Fortran code available according to the GNU General Public License v3.0 at https://doi.org/10.5281/zenodo.7535171 (Mailler, 2023).

All the available versions of the AerSett Fortran module are available according to the GNU General Public License v3.0 from https://doi.org/10.5281/zenodo.7535114 (Mailler et al., 2023).

**Data availability.** No data sets were used in this article.

**Author contributions.** All authors have contributed to the article by stirring the ideas, writing and correcting the paper, and providing bibliographical reference. SM has developed the Python scripts used to produce the plots.

**Competing interests.** The contact author has declared that none of the authors has any competing interests.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Acknowledgements.** The authors acknowledge the developers of the Python modules "Thermo" (Bell and Contributors, 2016-2021) and "SciPy" (Virtanen et al., 2020), used for our numerical treatment of the problem.

**Review statement.** This paper was edited by Samuel Remy and reviewed by two anonymous referees.

## References

Adebiyi, A. A. and Kok, J. F.: Climate models miss most of the coarse dust in the atmosphere, Sci. Adv., https://doi.org/10.1126/sciadv.aaz9507, 2020.

Betzer, P. R., Carder, K. L., Duce, R. A., Merrill, J. T., Tindale, R. W., Uematsu, M., Costello, D. K., Young, R. W., Feely, R. A., Breland, J. A., Bernstein, R. E., and Greco, A. M.: Long-range transport of giant mineral aerosol particles, Nature, 336, 568–571, https://doi.org/10.1038/336568a0, 1988.

Bell, C. and Contributors: Thermo: Chemical properties component of Chemical Engineering Design Library (ChEDL), GitHub [code], https://github.com/CalebBell/thermo (last access: 10 February 2023), 2016–2021.

Cheng, N.-S.: Comparison of formulas for drag coefficient and settling velocity of spherical particles, Powder Technol., 189, 395–398, https://doi.org/10.1016/j.powtec.2008.07.006, 2009.

Clift, R. and Gauvin, W. H.: Motion of entrained particles in gas streams, Can. J. Chem. Eng., 49, 439–448, https://doi.org/10.1002/cjce.5450490403, 1971.

Cunningham, E.: On the velocity of steady fall of spherical particles through fluid medium, P. Roy. Soc. Lond. A, 83, 357–365, https://doi.org/10.1098/rspa.1910.0024, 1910.

Davies, C. N.: Definitive equations for the fluid resistance of spheres, P. Phys. Soc., 57, 259–270, https://doi.org/10.1088/0959-5309/57/4/301, 1945.

Drakaki, E., Amiridis, V., Tsekeri, A., Gkikas, A., Proestakis, E., Mallios, S., Solomos, S., Spyrou, C., Marinou, E., Ryder, C. L., Bouris, D., and Katsafados, P.: Modeling coarse and giant desert dust particles, Atmos. Chem. Phys., 22, 12727–12748, https://doi.org/10.5194/acp-22-12727-2022, 2022.

Goossens, W. R.: Review of the empirical correlations for the drag coefficient of rigid spheres, Powder Technol., 352, 350–359, https://doi.org/10.1016/j.powtec.2019.04.075, 2019.

Jennings, S.: The mean free path in air, J. Aerosol Sci., 19, 159–166, https://doi.org/10.1016/0021-8502(88)90219-4, 1988.

Khvorostyanov, V. I. and Curry, J. A.: Fall Velocities of Hydrometeors in the Atmosphere: Refinements to a Continuous Analytical Power Law, J. Atmos. Sci., 62, 4343–4357, https://doi.org/10.1175/JAS3622.1, 2005.

Lapple, C. E. and Sheperd, C. B.: Calculation of particle trajectories, Ind. Eng. Chem., 32, 605–617, https://doi.org/10.1021/ie50365a007, 1940.

LeGrand, S. L., Polashenski, C., Letcher, T. W., Creighton, G. A., Peckham, S. E., and Cetola, J. D.: The AFWA dust emission scheme for the GOCART aerosol model in WRF-Chem v3.8.1, Geosci. Model Dev., 12, 131–166, https://doi.org/10.5194/gmd-12-131-2019, 2019.

Mailler, S.: FreeFall: v1: v1 (Version v1), Zenodo [code], https://doi.org/10.5281/zenodo.7535171, 2023.

Mailler, S., Menut, L., Cholakian, A., and Pennel, R.: AerSett (v1.0), Zenodo [code], https://doi.org/10.5281/zenodo.7535115, 2023.

Mallios, S. A., Drakaki, E., and Amiridis, V.: Effects of dust particle sphericity and orientation on their gravitational settling in the earth's atmosphere, J. Aerosol Sci., 150, 105634, https://doi.org/10.1016/j.jaerosci.2020.105634, 2020.

NOAA/NASA/USAF: U.S Standard Atmosphere 1976, Tech. Rep. NASA-TM-X-74335, NOAA-S/T-76-1562, https://www.ngdc.noaa.gov/stp/space-weather/online-publications/miscellaneous/
us-standard-atmosphere-1976/us-standard-atmosphere\_st76-1562_noaa.pdf (last access: 10 February 2023), 1976.

Rémy, S., Kipling, Z., Flemming, J., Boucher, O., Nabat, P., Michou, M., Bozzo, A., Ades, M., Huijnen, V., Benedetti, A., Engelen, R., Peuch, V.-H., and Morcrette, J.-J.: Description and evaluation of the tropospheric aerosol scheme in the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecasting System (IFS-AER, cycle 45R1), Geosci. Model Dev., 12, 4627–4659, https://doi.org/10.5194/gmd-12-4627-2019, 2019.

Ryder, C. L., Highwood, E. J., Walser, A., Seibert, P., Philipp, A., and Weinzierl, B.: Coarse and giant particles are ubiquitous in Saharan dust export regions and are radiatively significant over the Sahara, Atmos. Chem. Phys., 19, 15353–15376, https://doi.org/10.5194/acp-19-15353-2019, 2019.

Seinfeld, J. H. and Pandis, S. N.: Atmospheric chemistry and physics: From air pollution to climate change, Wiley-Interscience, ISBN 0471178152, 1997.

Shu, Q., Murphy, B., Pleim, J. E., Schwede, D., Henderson, B. H., Pye, H. O. T., Appel, K. W., Khan, T. R., and Perlinger, J. A.: Particle dry deposition algorithms in CMAQ version 5.3: characterization of critical parameters and land use dependence using DePoBoxTool version 1.0, Geosci. Model Dev. Discuss. [preprint], https://doi.org/10.5194/gmd-2021-129, 2021.

Sič, B., El Amraoui, L., Marécal, V., Josse, B., Arteta, J., Guth, J., Joly, M., and Hamer, P. D.: Modelling of primary aerosols in the chemical transport model MOCAGE: development and evaluation of aerosol physical parameterizations, Geosci. Model Dev., 8, 381–408, https://doi.org/10.5194/gmd-8-381-2015, 2015.

Stokes, G.: On the Effect of the Internal Friction of Fluids on the Motion of Pendulums, vol. 9, part II of Transactions of the Cambridge Philosophical Society, https://www3.nd.edu/~powers/ame.60635/stokes1851.pdf (last access: 10 February 2023), 1851.

van Boxel, J.: Numerical model for the fall speed of raindrops in a rainfall simulator, Tech. Rep. 1998/1, I. C. E. special report, https://hdl.handle.net/11245/1.523366 (last access: 10 February 2023), 1998.

van der Does, M., Knippertz, P., Zschenderlein, P., Giles Harrison, R., and Stuut, J.-B. W.: The mysterious long-range transport of giant mineral dust particles, Sci. Adv., 4, eaau2768, https://doi.org/10.1126/sciadv.aau2768, 2018.

Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, İ., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and SciPy 1.0 Contributors: SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nat. Methods, 17, 261–272, https://doi.org/10.1038/s41592-019-0686-2, 2020.

Zhang, L., Gong, S., Padro, J., and Barrie, L.: A size-segregated particle dry deposition scheme for an atmospheric aerosol module, Atmos. Environ., 35, 549–560, 2001.