Curve fitting algorithm for multimodal particle size distributions – a theoretical basis

Here we detail an open-source curve fitting algorithm for multimodal particle size distributions (MPSDs) and evaluate it against a ten-year dataset of ambient particle size distribution (PSD) measurements collected at Storm Peak Laboratory, a remote mountainous research site. This algorithm is grounded in traditional aerosol statistics and assumes measured particle distributions are the sum of several lognormal PSDs. It is designed to be free of any predefined mode templates or mode number constraints. For a MPSD measurement, the total number concentration (N_i), geometric standard deviation (σ_g), and geometric mean diameter (${\stackrel{\mathrm{‾}}{D}}_{\text{pg}}$) of each mode is estimated using a Levenberg–Marquardt nonlinear least-squares algorithm. These fitted modes are then iteratively subtracted from the measured PSD until convergence and/or accuracy thresholds are met. Rigorous evaluation of ambient aerosol data reveals a tri-modal distribution is a poor assumption for Storm Peak Laboratory, particularly during new particle formation events. Four or more modes were necessary for 55.7 % of data associated with new particle formation. Furthermore, the algorithm was used to characterize complex laboratory PSDs where size selected ammonium sulfate aerosol was coated in oxidized biogenic secondary organic matter. The multimodal algorithm with a tri-modal setting exhibited strong agreement with existing MPSD fitting models but with significantly improved computational time. In summary, this algorithm provides an effective method to analyze PSD datasets for in situ laboratory and ambient measurements. To improve accessibility of this algorithm to the broader aerosol research community, we also include supplemental functions to format datasets from common mobility particle size spectrometers.

1 Introduction

Atmospheric aerosols directly interact with Earth's climate by absorbing or scattering solar radiation, modulating the global radiative balance (Myhre et al., 2013). Furthermore, aerosol are critical to cloud formation by serving as cloud condensation nuclei (CCN) or ice nucleating particles (INPs) (Lohmann and Feichter, 2005). Aerosol-cloud interactions or indirect climate effects are commonly summarized as modifications to cloud properties such as cloud albedo or lifetime (Albrecht, 1989; Lohmann, 2007; Twomey, 1974; Twomey et al., 1984). Combined, these direct and indirect effects constitute the largest uncertainty in estimates of global radiation budgets (Myhre et al., 2013). Beyond climatic effects, atmospheric aerosol (or colloquially particulate matter (PM) in the scope of air quality) can be a significant air pollutant with adverse health effects. Aerosols have been demonstrated to negatively affect cardiopulmonary health, with fine particulate matter (PM_2.5) constituting an elevated risk (Pope and Dockery, 2006). The magnitude of health effects due to atmospheric aerosol is an area of active research and the following properties are typically the foci of studies: size, concentration, and composition (Fuzzi et al., 2015). Indeed, these properties have also shown to disproportionally influence the climatic effects of aerosol described above (Myhre et al., 2013).

Properties critical to studies measuring both the climatic and health impacts of atmospheric aerosol are size, concentration and composition. As such, measurements of the particle size distributions (PSDs) of aerosol are fundamental to understanding atmospheric aerosol. A common trait of PSDs is a multimodal distribution, indicative of a variety of atmospheric processes such as primary emission (e.g., mechanical generation such as dust or sea spray) or secondary formation (e.g., new particle formation from gas-phase). Analysis and parameterization of these modes is of particular interest to infer aerosol formation processes and composition.

Multimodal parameterizations of PSDs are determined using manual, semi-automated, or automated methods. Manual or semi-automated methods require some form of user control and interpretation to provide initial guesses of parameters (Mäkelä et al., 2000). While this manual control can be advantageous for particularly complex PSDs, it becomes a barrier when analyzing long-term multi-year datasets. Automated algorithms are available in both commercial (DistFit™, Chimera Technologies) and non-commercial form (Hussein et al., 2005; Taylor et al., 2014). The algorithm designed by (Hussein et al., 2005) has seen frequent implementation (also known as DO-FIT) in analyzing long-term datasets from ambient research sites (Franco et al., 2022; Herrmann et al., 2015; Ondráček et al., 2009). Among the automated algorithms, simultaneous fitting or solving for all modes is performed and typically limited to three or less modes in the sub-micrometer range (i.e. nucleation, Aitken, and accumulation). Furthermore, a recent evaluation of the DO-FIT model for multiple ambient settings determined high correlation with observed PSDs; however, tended to overestimate total number concentration by about 9 % (Zhu and Wang, 2024). While the methodology for the non-commercial methods is thoroughly described to allow replication, the technical skill necessary to replicate the algorithm may present a barrier for use by the broader aerosol research community.

This article presents an open-source, automated algorithm for parameterizing and fitting particle size distributions (PSDs) based on their theoretical descriptions, with an emphasis on ease of use. In contrast to existing methods where multiple aerosol modes are fitted simultaneously, this algorithm differs by using an iterative subtractive approach for fitting PSDs. To our knowledge, such an automated iterative subtraction technique has not previously been explicitly implemented in the aerosol sciences, although it has proven successful in characterizing geological sediment size distributions (Weltje and Prins, 2007). Another advantage of the technique presented is freedom from predefined mode templates or assumptions of maximum number of modes.

The algorithm was validated and developed with complex laboratory measurements containing six or more modes where size selected ammonium sulfate particles were coated in oxidized biogenic secondary organic matter. To evaluate performance for ambient aerosol distributions, the algorithm was used to characterize ten-years of 10 min averaged particle size distribution measurements from Storm Peak Laboratory (SPL), a remote mountainous research observatory in Steamboat Springs, Colorado (Hallar et al., 2025). Performance metrics for the multimodal classification of this dataset are described. Furthermore, the validity of a tri-modal distribution for SPL was tested, particularly during new particle formation (NPF) events (Gordon et al., 2017) which frequently occur at SPL (Hallar et al., 2011; Hirshorn et al., 2022). Lastly, we compared the multimodal algorithm with the existing DO-FIT model for ambient PSDs. In summary, the presented algorithm offers a robust, automated, and assumption-free approach for resolving complex multimodal aerosol particle size distributions across diverse experimental and observational datasets.

2 Methods

2.1 Laboratory Aerosol Generation and Measurement

Complex aerosol distributions consisting of a mixture of self-nucleated SOA from oxidized gas-phase precursors and monodisperse sulfate salts were used to develop this algorithm. Briefly, ammonium sulfate (AS, 99 %, A4915; Sigma-Aldrich) and bisulfate (ABS, 99.99 %, 455 849; Sigma-Aldrich) were atomized (Model 3076; TSI Inc., Shoreview, MN 3077) at a flow rate of 0.9 L min⁻¹ and dried using a Nafion™ dryer (Model MD-700-6S-3, PermaPure, Lakewood, NJ). Dried particles with an electrical mobility diameter of 300 nm were selected using a differential mobility analyzer (DMA, Model 3081A; TSI Inc., Shoreview, MN) with a 10:1 sheath-to-sample flow ratio. Particles were then sampled into a potential aerosol mass oxidation flow reactor (PAM-OFR, Aerodyne Research Inc., Billerica, MA) where gas-phase limonene (97 %, 183 164; Sigma-Aldrich) or α-pinene (99 %, 274 399; Sigma-Aldrich) was oxidized by $\mathrm{OH}/{\mathrm{O}}_{\mathrm{3}}$ producing both organic coatings on the sulfate particles and self-nucleated SOA. Following the PAM-OFR, O₃ was removed by an ozone denuder and aerosol were further dried by a desiccant dryer and Nafion™ dryer in series (Model MD-700-12S-3, PermaPure, Lakewood, NJ).

Particle size distributions (PSDs) were measured using a scanning electrical mobility sizer (SEMS, Model 2002; Brechtel Manufacturing Inc., Hayward, CA) also operating at a minimum 10 : 1 sheath-to-sample flow ratio. Sampling intervals were varied depending on aerosol generation stage to balance resolution and temporal variations; self-nucleation testing of SOA using a PAM oxidation flow reactor (30 s), monodisperse sulfate seed sampling (60 s), or coating evaluation (120 s).

For each experiment, 362 scans were obtained on average with coating evaluation scans representing a majority. A total of 16 554 SEMS scans were obtained and used to develop the curve fitting algorithm (Sect. 2.4–2.7).

2.2 Ambient Aerosol Measurement

Ambient particle size distribution measurements at SPL from 14 October 2010–24 October 2020 were used in evaluating the algorithm. PSD measurements were obtained using a scanning mobility particle sizer (SMPS, Model 3936; TSI Inc., Shoreview, MN 3077) and condensation particle counter (CPC, Model 3010; TSI Inc., Shoreview, MN 3077) at a time resolution of 5 min. To mitigate instrument and ambient noise, the data was averaged to 10 min intervals for algorithm development and evaluation on ambient data. For comparison purposes with the DO-FIT model, the high-resolution data (∼100 bins) was reduced to 32 bins and kept at a 5 min interval. Data was quality controlled to the standards set for level 1 data of the European Monitoring and Evaluation Programme (EBAS) database including multiple charge and diffusion corrections. Both visual and statistical NPF classifications for this time period were retrieved from a previous study (Hirshorn et al., 2022). Detailed specifications of the aerosol instrumentation, inlet, and site details of SPL are beyond the scope of this paper and are described elsewhere (Hallar et al., 2011, 2025; Petersen et al., 2019).

2.3 Theory

Atmospheric aerosol size distributions are commonly described as the sum of multiple lognormal PSDs. Assuming each mode is lognormally distributed and the PSD can be fully described by the summation of j lognormal distributions, a continuous multimodal PSD (MPSD) can be fully defined by Eq. (1) (Seinfeld and Pandis, 2016)

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }{n}_{\mathrm{N}}^{\circ }\left(\log D_{\mathrm{p}}\right)={\displaystyle \frac{\mathrm{d}N}{\mathrm{d}\log D_{\mathrm{p}}}}\\ {\displaystyle }& {\displaystyle }={\sum }_{i=\mathrm{1}}^j{\displaystyle \frac{N_i}{(\mathrm{2}\mathit{\pi }{)}^{\mathrm{1}/\mathrm{2}}\log {\mathit{\sigma }}_{\mathrm{g},i}}}\exp \left(-{\displaystyle \frac{(\log D_{\mathrm{p}}-\log {\stackrel{\mathrm{‾}}{D}}_{\text{pg},i}{)}^{\mathrm{2}}}{\mathrm{2}{\log }^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{g},i}}}\right),\\ {\displaystyle }& {\displaystyle }[n\cdot {\mathrm{cm}}^{-\mathrm{3}}]\end{array}\end{array}$$

where N_i is the total number concentration, σ_g,i is the geometric standard deviation, ${\stackrel{\mathrm{‾}}{D}}_{\text{pg},i}$ is the geometric mean diameter of the ith mode. Here, log is shorthand for the base 10 logarithm. A full description of the derivation of this equation and its parameters can be found elsewhere (Hinds and Zhu, 2022; Seinfeld and Pandis, 2016). Notation used throughout this article will be consistent with the notation used by Seinfeld and Pandis, 2016, with noted minor modifications made for indices such as j instead of n to represent each lognormal distribution (see Table 3).

2.4 Pre-Processing

Depending on the aerosol process of interest, the desired temporal resolution may differ. For example, a lower temporal resolution (e.g. 30 min averaging) is sufficient to characterize aerosol produced by stable laboratory generation. In this application, time averaging the PSD is beneficial in reducing variance caused by instrument noise. Conversely, the analysis of new particle formation requires a higher resolution to observe processes that can occur at or below an instrument measurement frequency (e.g. 2 min SMPS scan) and no time averaging is more appropriate. To accommodate a range of desired temporal resolutions, interval averaging is performed in minutes. Given a time resolved lognormal PSD with m binned diameter ranges ($D_{\mathrm{p}},D_{\mathrm{p}}+\mathrm{d}{D}_{\mathrm{p}}$) and t samples, the time averaged lognormal PSD, $\langle {\mathit{n}}_{\mathrm{N}}^{\circ }\left(\log D_{\mathrm{p}}\right)\rangle$ (Eq. 2), becomes a 1×m matrix

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }\langle {\mathit{n}}_{\mathrm{N}}^{\circ }\left(\log D_{\mathrm{p}}\right)\rangle =\langle {\mathit{n}}_{\mathrm{N},\mathrm{1}}^{\circ }\left(\log D_{\mathrm{p},\mathrm{1}}\right)\rangle ,\langle {\mathit{n}}_{\mathrm{N},\mathrm{2}}^{\circ }\left(\log D_{\mathrm{p},\mathrm{2}}\right)\rangle ,\mathrm{\cdots },\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{t}}{\sum }_{l=\mathrm{1}}^t{n}_{\mathrm{N},m,l}^{\circ }\left(\log D_{\mathrm{p},m}\right),[n\cdot {\mathrm{cm}}^{-\mathrm{3}}]\end{array}\end{array}$$

where $n_{\mathrm{N},m,l}^{\circ }\left(\log D_{\mathrm{p},m}\right)$ denotes the log-normalized number concentration for l time and m bin. Similarly, we obtain the time averaged differential particle number concentration 〈n_N(D_p)dD_p〉 by first multiplying $n_{\mathrm{N}}^{\circ }\left(\log D_{\mathrm{p}}\right)$ by dlog D_p and similarly applying Eq. (2).

2.5 Peak Parameterization

Multi-modal peak identification was performed using the findpeaks function from the Practical Numerical Mathematical Functions (pracma) R package (Borchers, 2011). For mode identification, the recommended minimum distance between peaks must be equal to or greater than 5 bins and are sorted by peak magnitude (local maximum number concentration) for each retrieval. To account for noise, findpeaks was applied to a smoothed cubic spline of the residuals for more robust peak identification. For each mode, lower D_p,min and upper D_p,max bounds of the identified peak mode diameter D_p,i are retrieved. For the purposes of initializing the curve fitting algorithm (see Sect. 2.6) it was assumed that $D_{\text{pg},i}\approx D_{\mathrm{p},i}$. The mode total number concentration N_i from the lower and upper diameters was calculated by the summation in Eq. (3)

$$\begin{array}{}\text{(3)}& N_i=\sum _{D_{\mathrm{p},\text{min}}}^{D_{\mathrm{p},\text{max}}}\langle {\mathit{n}}_{\mathrm{N}}\left(D_{\mathrm{p}}\right)\mathrm{d}{D}_{\mathrm{p}}\rangle ,[n\cdot {\mathrm{cm}}^{-\mathrm{3}}]\end{array}$$

For instances in which the predicted total concentration of an adjacent mode was overestimated, N_i may become negative at which case the subsequent peak is selected for fitting until the maximum number of iterations is reached. The mode specific geometric standard deviations σ_g,i used to initialize fitting were determined by Eq. (4)

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathit{\sigma }}_{\mathrm{g},i}=\\ & \sqrt{{\displaystyle \frac{\mathrm{1}}{N_i-\mathrm{1}}}\sum _{D_{\mathrm{p},\text{min}}}^{D_{\mathrm{p},\text{max}}}\langle {\mathit{n}}_{\mathrm{N}}^{\circ }\left(\log D_{\mathrm{p}}\right)\rangle \cdot (\log D_{\mathrm{p}}-\log D_{\mathrm{p},i}{)}^{\mathrm{2}}}\end{array}\end{array}$$

Similarly to the procedure used in calculating N_i, controls to prevent negative values from appearing within Eq. (4) were implemented where negative concentrations of $\langle {\mathit{n}}_{\mathrm{N}}^{\circ }\left(\log D_{\mathrm{p}}\right)\rangle$ are temporarily set to 0. If calculations of σ_g,i failed, the subsequent peak is tried for fitting until the maximum number of iterations is reached.

2.6 Curve Fitting

Using the parameters described above, each mode was fitted to a lognormal size distribution (Eq. 5) using the nlsLM function, a Levenberg–Marquardt nonlinear least-squares algorithm (LM-NLS) from the (minpack.lm) R package (Elzhov et al., 2022). A thorough description of LM-NLS is provided by (Moré, 1978). For completeness, the functional notation of a predicted lognormal PSD using LM-NLS is given by Eq. (5). Discretized values of D_p are modifiable depending on the size range of interest but must have the same resolution of the reported resolution in the data file (typically 10⁻² nm) to match actual and predicted binned concentrations. The LM-NLS algorithm outperformed standard non-linear least squares functions (e.g., Gauss-Newton Algorithm or GNA) for this application and avoided minimization problems arising from poor starting conditions, something commonly encountered using generic GNA-NLS methods. Poor mode fits are automatically rejected if any of the lognormal parameters $({\stackrel{\mathrm{‾}}{D}}_{\text{pg},i},N_i,{\mathit{\sigma }}_{\text{gi}})$ exceed or meet the significance threshold p=0.05.

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }{\widehat{\mathit{n}}}_{\mathrm{N},i}^{\circ }\left(\log D_{\mathrm{p}}\right)=f(D_{\mathrm{p}},{\stackrel{\mathrm{‾}}{D}}_{\text{pg},i},N_i,{\mathit{\sigma }}_{\text{gi}})\\ {\displaystyle }& {\displaystyle }=\mathrm{2.303}{\displaystyle \frac{N_i}{(\mathrm{2}\mathit{\pi }{)}^{\mathrm{1}/\mathrm{2}}\ln {\mathit{\sigma }}_{\mathrm{g},i}}}\exp \left(-{\displaystyle \frac{(\ln D_{\mathrm{p}}-\ln {\stackrel{\mathrm{‾}}{D}}_{\text{pg},i}{)}^{\mathrm{2}}}{\mathrm{2}{\ln }^{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{g},i}}}\right),\\ {\displaystyle }& {\displaystyle }[n\cdot {\mathrm{cm}}^{-\mathrm{3}}]\end{array}\end{array}$$

Each fitted mode ${\widehat{\mathit{n}}}_{\mathrm{N},i}^{\circ }\left(\log D_{\mathrm{p}}\right)$ is then iteratively subtracted from the previous residual ${\mathit{r}}_{i-\mathrm{1}}^{\circ }\left(\log D_{\mathrm{p}}\right)$. For the i=1 iteration, we define the initial residual distribution ${\mathit{r}}_{\mathrm{0}}^{\circ }$ as ${\mathit{r}}_{\mathrm{0}}^{\circ }\left(\log D_{\mathrm{p}}\right)=\langle {\mathit{n}}_{\mathrm{N}}^{\circ }\left(\log D_{\mathrm{p}}\right)\rangle$. The residual PSD after the assignment of j modes is then defined by Eq. (6)

$$\begin{array}{}\text{(6)}& {\mathit{r}}_j^{\circ }\left(\log D_{\mathrm{p}}\right)=\langle {\mathit{n}}_{\mathrm{N}}^{\circ }\left(\log D_{\mathrm{p}}\right)\rangle -\sum _{i=\mathrm{1}}^j{\widehat{\mathit{n}}}_{\mathrm{N},i}^{\circ }\left(\log D_{\mathrm{p}}\right).\end{array}$$

This is performed until any of the following conditions are met: fraction of variance unexplained (FVU) of ${\mathit{r}}_j^{\circ }$ is less than or equal to the FVU tolerance (Eq. 7), the maximum number of iterations has been reached, or the maximum number of modes allowed is reached. The second condition is utilized to prevent an “infinite loop” from occurring while the third is used to prevent overfitting. In cases where the estimated peak is much higher than the identified peak (1.5 times larger), the residuals are set to 0 to identify the next peak. While this will bias the FVU within the iterative loop, it is recalculated when evaluating the final model. Note the FVU can equivalently be defined as the square of the correlation coefficient of measured to predicted particle concentration subtracted from 1.

$$\begin{array}{}\text{(7)}& \text{FVU}={\displaystyle \frac{\text{var}\left[{\mathit{r}}_j^{\circ }\right(\log D_{\mathrm{p}}\left)\right]}{\text{var}\left[\langle {\mathit{n}}_{\mathrm{N}}^{\circ }\left(\log D_{\mathrm{p}}\right)\rangle \right]}}=\mathrm{1}-R^{\mathrm{2}}=\mathrm{0.01}\end{array}$$

Upon completion, the theoretical continuous MPSD (Eq. 1) can be approximated as the following summation of j fitted modes.

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\displaystyle }{n}_{\mathrm{N}}^{\circ }\left(\log D_{\mathrm{p}}\right)& {\displaystyle }\cong {\widehat{\mathit{n}}}_{\mathrm{N},j}^{\circ }\left(\log D_{\mathrm{p}}\right)\\ {\displaystyle }& {\displaystyle }={\sum }_{i=\mathrm{1}}^j{\widehat{\mathit{n}}}_{\mathrm{N},i}^{\circ }\left(\log D_{\mathrm{p}}\right),[n\cdot {\mathrm{cm}}^{-\mathrm{3}}]\end{array}\end{array}$$

2.7 Evaluation

Determining the accuracy or “goodness-of-fit” of the sum of j modes predicted by the LM-NLS algorithm is based on residual error analysis. The statistical metrics used to evaluate the model are the min-max normalized RMSE (NRMSE, see Eq. 9) and FVU of the fitted MPSD and actual data (Eq. 7), where k is an index of summation for matched m binned concentrations

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }\text{NRMSE}={\displaystyle \frac{\mathrm{1}}{\text{max}(\langle {\mathit{n}}_{\mathrm{N}}^{\circ }(\log D_{\mathrm{p}})\rangle )-\text{min}(\langle {\mathit{n}}_{\mathrm{N}}^{\circ }(\log D_{\mathrm{p}})\rangle )}}\\ {\displaystyle }& {\displaystyle }\times \sqrt{{\displaystyle \frac{\mathrm{1}}{m}}{\sum }_{k=\mathrm{1}}^m(\langle {\mathit{n}}_{\mathrm{N},k}^{\circ }(\log D_{\mathrm{p},k})\rangle -{\mathit{r}}_{i,k}^{\circ }(\log D_{\mathrm{p},k}){)}^{\mathrm{2}}}\end{array}\end{array}$$

For size-selected measurements, normalizing to the concentration range reduces the bias high concentration peaks (particularly relevant for monodisperse particle selection or nucleation particles) contribute to the RMSE. A successful MPSD fitting is considered when NRMSE<0.05 (5 %) and FVU<0.05 (5 %) by default, but are modifiable thresholds depending on the desired accuracy for a specific application. For the ambient SPL measurements presented, FVU and NRMSE thresholds were changed from the default 5 %–10 % to account for higher ambient noise. Examples and description of algorithm outputs are included in the Supplement.

3 Results

3.1 Laboratory Measurements

Measurements of the generated PSDs predictably exhibited a multimodal structure due to multiply charged particles, singly charged aggregates (i.e. doublets), and/or self-nucleation of SOA (see Figs. 1 and 2). Both multiply charged and singly charged aggregates were identified as expected as DMA-selected particles were passed through a charge neutralizer before entering the SEMS. The most prominent mode was the size selected inorganic sulfate mode with a mean particle diameter (D_p) of 300 nm (∼310 nm when coated in an oxidized organic layer). A nucleation and Aitken mode were also frequently observed corresponding to nucleation and subsequent agglomeration. For a modified experiment where inorganic sulfate particles with a mean particle diameter (D_p) of 150 nm were selected, singly charged aggregates up to singly-charged quadruplets were effectively characterized along with doubly and triply charged singlets (see Fig. 1).

Figure 1 Example of a successful multimodal curve fitting for a SEMS measurement of DMA-selected AS (D_p=150 nm) particles. Lognormal parameters for identified modes are listed in Table 1. Shaded gray region represents a NRMSE of 0.02 with a corresponding total particle uncertainty (RMSE) of 2537 $\mathrm{d}N/\mathrm{d}\log D_{\mathrm{p}}$ [$n\cdot {\mathrm{cm}}^{-\mathrm{3}}$].

Figure 2 Example of a successful curve fitting of a MPSD for size selected AS (D_p=300 nm) and α-pinene derived SOA aerosol. Lognormal parameters for identified modes are listed in Table 1. Shaded gray region represents a NRMSE of 0.03 with a corresponding total particle uncertainty (RMSE) of 2132 $\mathrm{d}N/\mathrm{d}\log D_{\mathrm{p}}$ [$n\cdot {\mathrm{cm}}^{-\mathrm{3}}$].

3.2 Ambient Measurements – SPL

The algorithm successfully characterized 59.4 % of the 208 060 10 min averaged scans at the accuracy thresholds of 0.1 (10 %) FVU and NRMSE. Under higher particle concentration conditions (e.g., $\ge {\mathrm{10}}^{\mathrm{5}}\mathrm{d}N/\mathrm{d}\log D_{\mathrm{p}}$), the percentage of successful fittings increases to 68.8 %. For these accuracy thresholds, the mean measured to predicted particle ratio was 0.991 and overall underpredicted concentrations by 0.899 %. This corresponded to a coefficient of determination (R²) of 0.999 when a simple linear regression was applied. The algorithm's ability to estimate total number concentration on a scan-by-scan basis is illustrated in Fig. 3. When all fittings are considered, including those below the thresholds defined above, the coefficient of determination (R²) of measured to predicted particle concentrations is 0.847.

Figure 3 Measured concentration vs. predicted concentration dN/dlogDp for fittings above the accuracy thresholds of 0.1 (10 %) FVU and NRMSE. Black dots represent prediction ratios for individual 10 min averaged scans at SPL from 14 January 2010–24 October 2020. The red line is the simple linear regression of the predicted and measured particle concentration with coefficient of determination and degrees of freedom (df) listed in the title.

Both the mean and median number of modes necessary to explain 90 % of the variance for all SPL measurements was 4±2 modes, where the error is standard deviations rounded to the nearest mode. Examples of measurement fittings where more than 4 modes were necessary for both low and high particle concentrations are given by Figs. 4 and 5 respectively. In the case of low particle concentrations, five modes were necessary to explain 92 % of particle concentrations. Six modes were needed for the high concentration case (e.g., visually identified NPF event with strong nucleation mode) where 99.57 % of all particles were predicted. Extending beyond these two examples, 37.8 % of all successful fittings required 5 or more modes.

Figure 4 Example of a successful curve fitting of a MPSD at SPL (22 January 2022 11:18:40 UTC) for non-NPF conditions. Lognormal parameters for identified modes are listed in Table 1. Shaded gray region represents a NRMSE of 0.07 with a corresponding total particle uncertainty (RMSE) of 108 $\mathrm{d}N/\mathrm{d}\log D_{\mathrm{p}}$ [$n\cdot {\mathrm{cm}}^{-\mathrm{3}}$].

Figure 5 Example of a successful curve fitting of a MPSD at SPL (30 April 2022 18:08:40 UTC) for the beginning of a NPF event. Lognormal parameters for identified modes are listed in Table 1 and modes 2–6 are not labeled due to low concentration compared to mode 1. Shaded gray region represents a NRMSE of 0.02 with a corresponding total particle uncertainty (RMSE) of 704 $\mathrm{d}N/\mathrm{d}\log D_{\mathrm{p}}$ [$n\cdot {\mathrm{cm}}^{-\mathrm{3}}$].

Table 1 A summary of the lognormal parameters specific to the modes in Figs. 1–4.

3.3 Evaluation of the Tri-Modal Assumption for NPF at SPL

For the time period spanning from 15 January 2010–22 October 2020, NPF events were visually classified on 372 d according to the methods described in (Hirshorn et al., 2022). The remainder were classed as non-NPF and consisted of either undefined or non-event days. Multimodal analysis was performed using the algorithm presented with fittings deemed successful if they met the following accuracy conditions: 90 % or more measured particles predicted and a NRMSE of 0.10 or lower attained. The relative frequency of modes necessary to fit each 10 min scan for both NPF and non-NPF event days are summarized in Fig. 6. For visually classified NPF days, 4+ modes are required to explain 55.7 % of scans. This requirement of higher order modes becomes more pronounced if using the statistical Gaussian NPF classification detailed by (Hirshorn et al., 2022). This method provides specific start and end times where active nucleation and growth are calculated. For this subset of data, both the mean and median number of modes necessary to explain 90 % of particle concentration increases to 5±2 modes where the error is standard deviations rounded to the nearest mode. Furthermore, 73.2 % of scans associated to statistically identified NPF required four or more modes. For all non-NPF PSDs, only 27.7 % of fits were tri-modal (Fig. 6). Regardless of the NPF classification method utilized, a tri-modal assumption fails to explain more than half of NPF events at SPL.

Figure 6 Relative frequency of modes (%) necessary to successfully fit a PSD (NRMSE≤10 % and FVU≤10 %) on days where NPF was visually classified (top) or days where NPF was confirmed to not have occurred (bottom). Annotations indicate the percentage of data associated with 4 or more modes.

3.4 Comparison to the DO-FIT model

To robustly compare the multimodal algorithm against an established fitting technique, we performed a comparison with the DO-FIT model which has been extensively used in characterizing PSD measurements (Hussein et al., 2005; Zhu and Wang, 2024). Both multimodal and DO-FIT were allowed to fit a maximum of three modes. In the case of DO-FIT, mode identification was based on pre-existing parameter schemes based on ambient measurements (Hussein et al., 2005). Both models were run (or rerun in the case of multimodal) on a 5 min time resolution dataset of SPL data (15 January 2010–22 October 2020) at a lower bin resolution scheme (32 bins) to reduce computational cost. In total, 411 775 scans were examined.

First, we note that the multimodal algorithm was significantly faster than the DO-FIT model, requiring 7.5 and 7784.13 single-threaded CPU hours, respectively. Model fitting failure, defined here as either automatic self-termination or inability to fit measured data, was comparable between methods with multimodal and DO-FIT accounting for 9.96 % and 5.86 %, respectively. Overall fitting accuracy metrics are summarized in Table 2, where values containing mean values were trimmed such that extreme values in the upper and lower 5 % were removed to account for computational outliers.

Table 2 Model comparison metrics. All metrics are reported in units of $\mathrm{d}N/\mathrm{d}\log D_{\mathrm{p}}$ [$n\cdot {\mathrm{cm}}^{-\mathrm{3}}$], except for NRMSE and R² which are dimensionless.

Table 3 List of symbols.

To summarize model fitting in lognormal parameter spaces, the probability distribution of each model P(N,D_pg) was estimated with a normalized kernel density estimate (KDE) using the Modern Applied Statistics with S R package (Venables et al., 2002). The resulting estimates are shown in Fig. 7. In both cases, the accumulation mode is most distinct, with agreement in location and density between both. Primary differences were observed in the nucleation and Aitken regions, where multimodal was most likely to assign modes between 20 and 70 nm for N between 10² and 10⁴ cm⁻³.

Figure 7 KDE estimate of the probability distribution of each lognormal fitting algorithm (as labelled) as a function of N and D_pg. Saturated regions of each panel indicate regions where the model is most likely to assign modes for the ambient PSDs collected at SPL from 15 January 2010–22 October 2020.

Residuals for each model were compared across the full reduced dataset to examine the relationship between model prediction error and aerosol regimes (i.e. within D_p and N space). Data was pre-binned by existing diameters and further grouped by log-spaced bins of N. Given over or underprediction in aerosol models can result in predicted particle concentrations spanning several orders of magnitude, we represent the magnitude and direction of predictive error using the log bias ratio (Eq. 10):

$$\begin{array}{}\text{(10)}& \text{Log Bias Ratio}={\log }_{\mathrm{10}}\left({\displaystyle \frac{\widehat{y}}{y}}\right)\end{array}$$

In Fig. 8, the mean of the log bias ratio across all scans is shown with computational errors removed (i.e. extreme non-physical values ±5 %). For both models, underprediction of the SPL dataset was most frequent. Overprediction was observed only for high aerosol loadings. At low concentrations (i.e. N less than 10² cm⁻³) DO-FIT outperformed multimodal where automated peak identification is difficult. DO-FIT also outperformed multimodal for both sub 20 nm bins ($N<{\mathrm{10}}^{\mathrm{3}}{\mathrm{cm}}^{-\mathrm{3}}$) and data between 100 and 200 nm irrespective of N. For all other regions, multimodal generally outperformed DO-FIT.

Figure 8 Contour map of the magnitude in prediction error of the multimodal (left) and DO-FIT (right) PSD fitting algorithms. Red indicates overprediction whereas blue indicates underprediction, with the magnitude of error given in log scale (i.e. a value of −2 indicates the model underpredicts in this PSD region by 100×). Gray fill indicates regions where fitting failed.

4 Discussion

4.1 Tri-Modal Assumption

Two to three modes are generally sufficient to characterize ambient PSDs for a variety of environments (Franco et al., 2022; Herrmann et al., 2015; Hussein et al., 2005; Mäkelä et al., 2000; Zhu and Wang, 2024). As such, existing commercial algorithms such as DistFit™ (Chimera Technologies) are designed to use a maximum of three modes when fully automated. Similarly, existing non-commercial algorithms assume (or merge smaller modes) to achieve at most a trimodal distribution (Hussein et al., 2005; Mäkelä et al., 2000; Zhu and Wang, 2024). For the algorithm developed by (Hussein et al., 2005), modes are assigned by iterating lognormal parameters (N_i,σ_g, and ${\stackrel{\mathrm{‾}}{D}}_{\text{pg}}$) in a predefined range. In contrast, the algorithm we developed does not have any assumptions on mode parameters and is only limited by user-defined arguments for maximum allowed modes and accuracy tolerances. After iteratively subtracting modes in order of decreasing concentration, the algorithm will continue to assign modes unless desired tolerance is achieved or maximum modes are reached. This flexibility proved invaluable when evaluating both the laboratory measurements at Purdue University and ambient measurements at SPL. Existing algorithms with limited mode numbers would not have been able to characterize complex laboratory measurements (e.g. DMA size-selected) such as those presented here where a minimum of six modes is present. Furthermore, the assumption of a tri-modal distribution was inadequate for PSD data at SPL from 2010–2020, particularly during NPF events. For statistically classified NPF, ∼75 % of PSDs required four or more modes to explain 90 % of all particle variance. When considering all fittings at the same accuracy, 4±2 modes were required.

4.2 Accuracy

Evaluation of the algorithm against SPL data revealed that the algorithm is highly effective in fitting ambient PSDs. When considering multimodal fittings without any accuracy thresholds for ten years of SPL PSD data, the model exhibited high predictive skill (R²=0.847) for total predicted vs. measured particle concentration. When applying more selective accuracy thresholds to filter poor/failed fittings, the algorithm's predictive skill increased at the expense of a lower number of available fittings. Using the accuracy thresholds NRMSE≤10 % and FVU≤10 %, the algorithm slightly underestimated total particle concentration at SPL by 0.9 % on average. Zhu and Wang, 2024, determined the multi-lognormal particle distribution using 3 or less modes overpredicted concentrations by 9 % when examining hourly data from eight diverse research sites across Europe with a high correlation of observed and fitted results (R²>0.75). If we consider the results of Zhu and Wang, 2024 for Puy de Dôme only (which is most similar to SPL among the research sites they examined), a bi- or tri-modal distribution was accurate in fitting nucleation mode particles but performed very poorly otherwise (R²<0.12). In comparison to the DO-FIT model for a reduced dataset on SPL, we show when similarly limited to a tri-modal assumption, the multimodal algorithm achieved higher correlation with measured data (R²=0.722 vs. R²=0.616). Our results show using higher order modes for a similarly mountainous forested region was highly accurate.

4.3 Limitations

For the laboratory lognormal PSDs used to develop this algorithm, scans with total number concentrations exceeding ${\mathrm{10}}^{\mathrm{10}}\mathrm{d}N/\mathrm{d}\log D_{\mathrm{p}}$ or with unstable sheath flows (high standard deviation from set flow) were removed. Inclusion of these scans resulted in large overpredictions in concentration for all modes. An inherent limitation of any curve fitting procedure is the requirement that a local maximum is identified. For measured nucleation modes, the lower bound is typically the peak diameter with the lower tail outside the effective measurement range of instruments. Thus, the model will fail to curve-fit nucleation modes with a mode at or below the lower end of the measurement range. Low concentrations additionally preclude peak identification, preventing curve fitting (Fig. 8). Failure to identify edge cases of multiply charged particles or aggregates (ex. highest charged singlet or largest aggregate particle) were also encountered if a FVU tolerance of 0.01 (1 %) or higher was selected. Failure of edge cases is best summarized in Fig. 8. In the example presented in Fig. 2, a 0.1 % tolerance was used so all multiply charged particles and agglomerates were accounted for, emphasizing the strength of user-adjustable tolerance thresholds.

Significance testing of each parameter of the fitted lognormal PSD retrieved from the LM-NLS algorithm was successful (p≤0.05) with model runs containing parameters above this threshold automatically removed. We strongly caution however that MPSDs with ∼7 or more fitted modes will result in residual modes with low degrees of freedom e.g., significance testing is not an appropriate measure for goodness-of-fit. However, even with the increased statistical power for the sum of j modes, we found significance testing using either parametric or non-parametric testing to be a poor metric in evaluating the performance of the multimodal curve fitting algorithm in fitting a MPSD and recommend using the combination of the FVU and NRMSE.

Figure 9 Flowchart summarizing the data formatting process to apply the methods described here to a variety of datasets.

4.4 Applications

This multimodal curve fitting algorithm is open-source and available to the aerosol research community as a documented multi-argument R function. Both an actively maintained version and original copy (as submitted) are publicly accessible in addition to an example file describing the process in formatting datasets for analysis and applying the multimodal R function (see Supplement). This algorithm is not specific to any instrument or file format but must contain binned log-normalized concentration data (ideally corrected for any particle losses) with numeric names corresponding to the diameter range (e.g., 10.16 nm) listed in the data file. An assumption of this algorithm is that the data to be analyzed has been corrected for instrument errors and has been quality controlled. It is assumed the data can be pre-processed according to the procedure described in Sect. 2.4. To assist in this step, four supplemental functions are provided to pre-process TSI or Brechtel measurement data (instrument software version dependent). Furthermore, both the NASA-Ames (level 1 and 2) and netCDF standardized formats are accepted (see Fig. 9). Other instrument files or formats not listed can additionally be included upon request.

5 Conclusions

This multimodal curve fitting algorithm provides a method to estimate and parameterize MPSDs using familiar theory foundational to aerosol statistics. As MPSDs are assumed to be the summation of j modes, the algorithm iteratively subtracts modes in order of maximum concentration until user-defined accuracy and convergence tolerances are achieved. The algorithm is not limited by predefined ranges of lognormal parameters or number of modes and was designed for analyzing both laboratory and ambient PSD measurements. In particular, the freedom from predefined mode templates was invaluable for characterizing ten-years of ambient PSDs measured at SPL. The commonly assumed tri-modal distribution for ambient aerosol particles was adequate for only 27.7 % of non-NPF data analyzed where 48.4 % of data required 4 or more modes. During periods of statistically defined NPF this assumption became increasingly unreliable, where approximately 75 % of PSDs required four or more modes to explain 90 % of particle variance with a mean/median of 5±2 modes. A tri-modal limited run of the multimodal algorithm showed similar agreement with the DO-FIT algorithm but achieved higher performance metrics (Table 2). We note that these improved metrics can be attributed to the purely mathematical nature of multimodal whereas DO-FIT is arguably more grounded by decades of MPSD observations. We believe this research tool can be utilized by the aerosol research community to analyze aerosol size distributions and more critically, improve our understanding of the size-specific effects aerosol have on air quality and climate.