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.
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.
: accumulation mode,
: coarse mode,
: 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 Cd 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.
where Cd is the drag coefficient, Ap 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, Vp the particle volume and g the acceleration of gravity.
where the Reynolds number Re is defined as follows:
where and S=110.4 K.
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.
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.
To slightly generalize matters, let us rewrite Eq. (5) as
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 and δ as
The terminal Reynolds number of the particle is equal to
Let us introduce
where R is the Reynolds number of the particle if it were to fall at speed : 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:
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 Cd from (Eq. 7). Independent of the formulation of f(Re), the relative deviation δ of v∞ from 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:
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).
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
This expression of δ(R) yields an error relative to the exact solution <1 % up to Re=10 and <2.5 % up to Re=103 (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=103, 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 µ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).
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:
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):
From the Knudsen number, slip correction factors can be designed to account for the non-continuous effects, turning Eq. (1) into
we obtain the following:
The terminal Reynolds number of the particle is equal to
Let us introduce:
With these notations, , and Eq. (24) becomes
Equation (27) is the same fixed-point equation as Eq. (13) but with parameter (defined from Eq. 26) instead of R (from Eq. 12). Therefore, its solution is , 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 Cc 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:
Calculate v∞ as
Figure 4 shows the impact of the slip correction term and the drag correction term on . Two regimes are clearly separated in this figure, with 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 is less than 5 %. This zone in which the Stokes equation is directly applicable covers the diameter range at ground-level pressure and the range at 200 hPa. Some authors (e.g., Mallios et al., 2020) argue that the slip correction factor Cc 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), Cc 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∞.
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 starting from and successively cut the solution interval [Rem;ReM] into halves until the relative error on speed becomes smaller than ε: ; then .
Fixed-point method. Find the solution δ to Eq. (27) starting from δ0=0 and iterate it with until the relative error on speed becomes smaller than ε: . Then, .
AerSett method. Directly calculate .
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 108 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 we have , so that for all three methods we speed up the calculation by assuming that when . As can be seen in Fig. 4 (see the 0.99 iso-line of ), 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 . 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 () to 7.8 () for the bisection method, and from 1.5 () to 7.4 () 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.
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∞:
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 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.
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).
No data sets were used in this article.
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.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This paper was edited by Samuel Remy and reviewed by two anonymous referees.
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. a
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. a
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. a, b, c, d, e
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. a, b, c, d, e, f
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. a
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. a
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. a, b, c, d
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. a, b, c
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. a
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. a, b, c, d, e
Seinfeld, J. H. and Pandis, S. N.: Atmospheric chemistry and physics: From air pollution to climate change, Wiley-Interscience, ISBN 0471178152, 1997. a
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. a
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. a
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. a, b, c
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. a, b, c
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. a, b, c, d
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. a