Producing aerosol size distributions consistent with optical particle counter measurements using space-based measurements of aerosol extinction coefficient

Stratospheric aerosol has been observed by several long-lived observational systems. These include the University of Wyoming series of balloon-borne optical particle counters (OPCs) (1971–2020) and the Stratospheric Aerosol and Gas Experiment (SAGE) series of instruments, particularly SAGE II (1984–2005). Inferences of aerosol surface area density (SAD) and volume density are straightforward using data from OPCs. Conversely, many numerical methods to infer size distributions and SAD have been applied to SAGE II observations, but all are limited by the restricted number of independent wavelengths of the SAGE optical measurements. We have developed a new method that uses OPC observations to constrain SAGE II inferences of aerosol properties. We start by noting that, whatever the details of the underlying size distribution, the SAGE II-measured aerosol extinction coefficient ratio (525–1020 nm) must reflect the shape of the underlying aerosol size distribution for particles that dominate the extinction coefficient values (roughly radii from 0.1 to 0.5 µm). Since this extinction ratio can be easily calculated from OPC measurements, we use the OPC size distribution measurements, across a broad range of aerosol levels from background to highly volcanic, to compute the associated 525–1020 nm extinction coefficient ratios for each measurement. We then sort the OPC measurements by these ratios (across a range of roughly 1–6) into discrete ratio bins and derive mean bimodal log-normal size distributions for each bin using a particle swarm optimization. These fits can be applied to SAGE II observations without the need for further retrieval calculations, effectively producing an OPC-like product consisting of the six bimodal parameters for all SAGE II observations. This method successfully captures the median behavior of the OPC inferences of bulk parameters like aerosol surface area and volume densities, although we also observe a significant altitude dependence particularly in the lower stratosphere. In addition, there are occasional deviations of SAD from the fit behavior by as large as a factor of 10 for individual OPC measurements of SAD, primarily due to variations in small-radii particle number density (roughly smaller than about 0.15 µm). The presence of such particles is effectively invisible to extinction coefficient measurements such as those by SAGE II.

1 Introduction

Space-based measurements of stratospheric aerosol extinction coefficient have been made continuously since the 1978 launch of the Stratospheric Aerosol Measurement (SAM II) instrument aboard the Nimbus-7 spacecraft. For practically the same period of time, these measurements have been used to infer underlying properties of the aerosol focused on the aerosol size distribution (ASD) and properties that impact chemistry and climate such as aerosol surface area density (SAD) and total aerosol mass or volume density (VD). Among the earliest efforts to infer ASD using space-based measurements was use of the 450 and 1000 nm aerosol extinction coefficient measurements by the original Stratospheric Aerosol and Gas Experiment instrument (SAGE, 1979–1981) to fit a single-mode log-normal fit of fixed width (1.6) and inferring the mode radius and number density (Yue and Deepak, 1983). Since the stratospheric aerosol size distribution tends roughly to resemble a single-mode log-normal (SLN) distribution (Pinnick et al., 1976), though other mathematical forms exist, the SLN remains a common starting point for many ASD algorithms based on space-based observations of aerosol extinction coefficient (e.g., Russell et al., 1996; Arfeuille et al., 2013; Nyaku et al., 2020; Knepp et al., 2023; Wrana et al., 2021) often making use of the long-lived SAGE II mission (1984–2005). SAGE II measurements, with aerosol extinction coefficient measurements at four wavelengths (385, 452, 525, and 1020 nm), remain an object of considerable scientific attention given that they include observations of part of the recovery from the 1982 eruption of El Chichón; the 1991 eruption of Mt. Pinatubo and its recovery; and a relatively volcanically quiet, low-aerosol-loading period from 2000 to the end of its mission in 2005. As a result, this record remains a core part of the long-term stratospheric aerosol record (Kovilakam et al., 2023) and still plays a significant role in the study of the impact of volcanic activity on climate and the processes that lead to ozone destruction (Rieger et al., 2020; Revell et al., 2017; Pauling et al., 2023).

After 40 years, it might be expected that there would be a generally accepted approach (or approaches) from which robust determinations of ASD from the SAGE II measurements are routinely inferred, and this topic would be of limited further effort; however, this is not the case. One reason for the proliferation of diverse methods for solving for aerosol properties is that, while the retrieval methods are almost always mathematically straightforward, all retrieval methods effectively constrain the ASD solution space such that not all mathematically (as opposed to physically) plausible solutions for ASD are allowed. This is necessary because the measurements contain insufficient information to uniquely identify the ASD across the span of radii that are relevant to chemistry and climate. For instance, Thomason et al. (2008) showed that SAGE II aerosol extinction measurements have an explicit minimum SAD that is consistent with observations but an upper limit that is unbound (effectively infinite) by allowing very large numbers of very small and very ineffective scattering particles. Retrieval algorithm constraints often take the form of a fixed mathematical form for the aerosol size distribution (e.g., the SLN) or constrains the way aerosol number density varies as a function of particle radius (e.g., smoothness). Most a priori constraints are not at face value unreasonable, but, when applied, they fundamentally affect the outcomes for ASD and bulk property retrievals in ways that are not easy to account for in an uncertainty estimate and can vary substantially from constraint to constraint. Unsurprisingly, these factors point out the lack of robustness in SAGE II-based ASD and other aerosol property retrievals for which no solution based solely on the measurements is possible.

Given the limited number of SAGE II measurement wavelengths and the correlation between, particularly the short-wavelength, channels (385, 452, 525 nm), inferring an ASD more complex than an SLN (with three free parameters) such as a bimodal size distribution (with six) is impossible. However, in situ measurements of stratospheric aerosol size distributions are often more complex than an SLN. For instance, the University of Wyoming optical particle counter (WOPC) measurements are often fit best with a bimodal log-normal aerosol size distribution (Deshler et al., 2003, 2019). Both modes do not necessarily contribute significantly to a computed aerosol extinction coefficient at SAGE II wavelengths. In many cases, the computed extinction at SAGE II wavelengths does not differ significantly between a fitted SLN distribution and a bimodal distribution. Both modes, however, are often important in the estimation of aerosol bulk properties like SAD, which can be dependent on small particles, primarily missed by optical extinction measurements (Deshler et al., 2003). While we will make extensive use of WOPC data to infer aerosol bulk properties from SAGE II retrievals, this is not primarily a validation or intercomparison of the measurements of these two instruments which have been presented elsewhere (Hervig and Deshler, 2002; Deshler et al., 2003, 2019; Kovilakam and Deshler, 2015). To a greater degree, we are trying to determine whether it is possible to infer the magnitudes and variability observed in WOPC-derived key parameters like SAD from SAGE II measurements. In that regard, we are treating the WOPC measurements as a test bed for SAGE II retrievals.

Therefore, in this paper, we discuss the SAGE II measurements and demonstrate some of the limiting factors for ASD inferences. We then estimate aerosol extinction coefficient at SAGE II wavelengths and SAD using the WOPC data alone and show how WOPC SAD varies with computed aerosol extinction and its wavelength dependence. The degree to which this relationship is well-behaved directly addresses how well the SAGE II measurements can be used to infer ASD or SAD consistent with WOPC values. We demonstrate that, while the median behavior of WOPC observations can be replicated, there remain substantial SAD and VD positive outliers, primarily in the lower stratosphere, that are larger than the median value by factors as large as 10. While it applies only to the median behavior, we produce a WOPC-based bimodal log-normal ASD, varying with aerosol extinction coefficient wavelength dependence, that potentially allows a bimodal aerosol size distribution to be assigned to any SAGE II multi-wavelength stratospheric aerosol extinction coefficient measurement set. While the analysis and its outcomes are strictly only relevant to WOPC/SAGE II comparisons in the midlatitudes of the Northern Hemisphere, the outcomes reflect a fundamental limitation on what is possible for aerosol property estimates from SAGE-like measurements. Work is ongoing now to extend this analysis to SAGE III/ISS.

2 Some issues related to estimating size distribution using SAGE II data

The SAGE II aerosol extinction coefficient ensemble consists of measurements at four wavelengths (385, 452, 525, and 1020 nm) that usually extend from the upper troposphere to 40 km. Assuming Mie scattering, the aerosol extinction coefficient, k_λ, can be mathematically expressed as

$$\begin{array}{}\text{(1)}& {\displaystyle }{k}_{\mathit{\lambda }}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathit{\pi }{r}^{\mathrm{2}}{Q}_{\mathit{\lambda }}(r,m_{\mathit{\lambda }}){\displaystyle \frac{\mathrm{d}n\left(r\right)}{\mathrm{d}r}}\mathrm{d}r\end{array}$$

or

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }{k}_{\mathit{\lambda }}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{3}}{\mathrm{4}r}}{Q}_{\mathit{\lambda }}(r,m_{\mathit{\lambda }}){\displaystyle \frac{\mathrm{d}V\left(r\right)}{\mathrm{d}r}}\mathrm{d}r,\end{array}$$

where Q_λ is the Mie kernel for particles of radius r and with an index of refraction m_λ, and the aerosol size distribution is $\mathrm{d}n\left(r\right)/\mathrm{d}r$ in number per unit radius and $\mathrm{d}V\left(r\right)/\mathrm{d}r$ in volume of aerosol per unit radius. Figure 1a shows the values of Q_λ(r,m_λ) using a refractive index typical of stratospheric conditions and sulfuric acid–water aerosol. Figure 1b shows the per-unit-volume kernels weighted by the measurement wavelength or $\mathrm{3}{Q}_{\mathit{\lambda }}(r,m_{\mathit{\lambda }})/\mathrm{4}r$. From these figures, it is clear that the SAGE II measurement ensemble does not contain significant information for particles much less than 0.1 µm, even for the shortest wavelength measurement. As a result, estimates of total number density and similar parameters dependent on low-order moments of the aerosol size distribution (e.g., number density) are not well constrained by the measurements and, in fact, depend on how a retrieval process fills this information gap. Thus, while SAGE II measurements can almost always be used to find a unique log-normal (or similar low-free-parameter) aerosol size distribution that reproduces the measurements, there is no guarantee that all high-value geophysical parameters like SAD will be adequately calculated.

Figure 1 (a) Mie extinction kernels and (b) kernels scaled to per-unit-aerosol volume times wavelength for the SAGE II channels assuming spherical water–sulfuric acid droplets at 220 K and a composition of 75 % H₂SO₄ and 25 % H₂O. The real refractive indices used were 1.432, 1.432, and 1.421 for 453, 525, and 1020 nm, respectively, with zero for all imaginary parts.

There are further complications in performing SAGE II ASD retrievals. The measurements at 385 nm are not considered reliable, except at relatively high extinction coefficient values ($>{\mathrm{10}}^{-\mathrm{3}}{\mathrm{km}}^{-\mathrm{1}}$) (Thomason et al., 2018), and are not recommended for general use, reducing available measurements for ASD retrievals to only three wavelengths. In addition, as Fig. 1b shows, the two remaining short-wavelength measurements (452 and 525 nm) have significant overlap in their extinction kernels and thus provide limited unique information between them, particularly in light of their associated uncertainties. This can be demonstrated with a relatively simple exercise. First, we compute an Ångström coefficient using SAGE II extinction coefficient data at 525 and 1020 nm and then use this value to extrapolate to aerosol extinction coefficient at 452 nm using one of the former measurements as the base. While the mean difference between the estimated and measured aerosol extinction coefficient is primarily a measure of how well the extrapolation works, its relative rms (root-mean-square) mean difference is a measure of how much unique information exists in the 452 nm measurement. Crudely, if the relative rms mean difference is greater than the 452 nm measurement uncertainty, then it is possible that there is some usable additional information contained in the measurement.

In Fig. 2, we show the outcome for April 1999 with measurement locations north of 20° N where observations at all three wavelengths exist and where extinction coefficient is less than 0.01 km⁻¹ (a crude cloud filter) and greater than ${\mathrm{10}}^{-\mathrm{5}}{\mathrm{km}}^{-\mathrm{1}}$, below which measurement quality decreases rapidly. In this figure, we show the mean ratio between the predicted extinction coefficient and the measured extinction coefficient at 452 nm as a function of altitude (dashed line, scale at the top). On average, the Ångström extrapolation does well between 13 and 24 km where the predicted value is within 5 % of the observed values. The departures increase to about +10 % at 10 km and −20 % at 30 km, primarily demonstrating the limitations in the interpolation method. This figure also shows that the relative rms mean difference (solid line, bottom scale) between estimated and measured aerosol extinction coefficient, a stand-in for inferred noise, is routinely about half the size of the reported measurement uncertainty (dotted line, bottom scale). In this context, the differences between reported and inferred uncertainties are likely somewhat exaggerated due to the correlation in measurement uncertainty, particularly among SAGE II's short-wavelength channels and, to a lesser extent, the fact that the reported uncertainties contain both systematic and precision elements (Damadeo et al., 2013). Nonetheless, the degree to which the 452 nm channel can be inferred from the values at 525 and 1020 nm strongly suggests that either there is limited variability in aerosol size distribution for particles which control the 525–1020 nm extinction coefficient ratio or, if significant variability does exist, the ability of the 452 nm channel to illuminate that variability is very low.

Figure 2 For the Northern Hemisphere (>20° N) in April 1999 at 452 nm: (1) ratio of estimated to measured aerosol extinction coefficient (dashed, top scale), (2) relative rms difference between estimated and measured aerosol extinction coefficient (solid, bottom scale), and (3) median relative aerosol extinction coefficient measurement uncertainty (dotted, bottom scale).

If the ability of the 452 nm channel to illuminate variability in the ASD is low, as we will show below, then fitting a meaningful low-parameter-size distribution, like an SLN, is problematic. This assessment is corroborated by past efforts to infer size distributions from these measurements in which single-mode log-normal fits with SAGE II data produce distributions that are rather narrow (e.g., Wang et al., 1989), as well as findings that it is possible to fit the extinction coefficient measurement spectra using a vanishingly narrow distribution (a delta function) (Thomason et al., 2008). We conclude that the SAGE II aerosol extinction measurement ensemble has at best two pieces of information that are most clearly represented in the overall magnitude of extinction at 525 and 1020 nm and their extinction coefficient ratio (or the extinction coefficient spectral slope). This is in basic agreement with assessments of the information content of the measurements (e.g., Thomason et al., 1997). As well, a similar conclusion was reached by Thomason and Poole (1993) using a different technique. With so little information contained in the measurements, essentially all SAGE II aerosol size distribution retrievals have little recourse except to be dependent on the retrieval method. In other words, the outcomes from the retrieval process are likely controlled by the assumptions made at the outset of the effort, and the robustness of the inferences is debatable.

While we are focused on SAGE II, it is worth considering whether other measurement types such as the limb-scatter technique employed by OSIRIS (Optical Spectrograph and Infra-Red Imaging System) can be used to infer aerosol size distributions in a similar way. It is clear from Fig. 1 that adding short-wavelength measurements (e.g., at 385 nm) would increase the information about the small particles present. Practically, however, robust measurements of stratospheric aerosol extinction coefficient at wavelengths much shorter than the 385 nm channel on SAGE II are difficult due to the effects of molecular scattering and absorption by ozone and other gases. Simply increasing the number of visible and near-infrared measurements may, through repetitive information, improve the resolution of size distribution retrievals in the 0.1–0.5 µm radius range; however, the degree to which this is true depends on the precision of the measurements and the details of the measurements' spectral location. Still, such an instrument (e.g., the current SAGE III/ISS instrument) may not radically improve the ability to infer aerosol size distributions which could represent the small particles (<0.1 µm radius). A possible long-term solution would be to add measurements where sulfuric acid aerosol strongly absorbs in the infrared. The sulfuric acid extinction kernels in the infrared are relatively flat across radii relevant to stratospheric aerosol and are, thus, a near, but not exact, measure of the total volume of aerosol present in the measurement volume (Thomason, 2012). Combined with visible and near-infrared measurements (where scattering dominates), infrared measurements could provide some constraint to what occurs at smaller particle sizes (Thomason, 2012; Boone et al., 2023), though most likely this would still require significant constraints to the possible solutions. Volume is inherently insensitive to small particles, in contrast to surface area.

A characteristic of almost all size distribution retrievals is that they tend to stand on their own, apart from information arising from in situ measurements of aerosol size distributions (e.g., by optical particle counters) beyond very general considerations such as that the aerosol size distribution is generally compatible with single or multimodal log-normal aerosol size distributions. Given the weak information content of the SAGE II measurements, we now consider the possibility of using in situ information much more explicitly. In the following sections, we will examine the ability to use the in situ measurements from the WOPC to assess the variability in size distributions as a function of the 525–1020 nm extinction coefficient ratio and attempt to infer WOPC-compatible aerosol size distributions from SAGE II extinction coefficient measurements. We will evaluate the success of this effort primarily by how well such inferences can reproduce WOPC-like SAD values.

3 University of Wyoming OPC measurements

The University of Wyoming (UW) in situ balloon-borne measurements of aerosol size distributions have been made continuously since 1971 (Deshler et al., 2003; Deshler, 2023). Vertical profiles of size-resolved cumulative aerosol concentration are provided along with unimodal/bimodal log-normal fits. The number density profiles are provided at full resolution and 0.5 km resolution, the size distribution fits at 0.5 km resolution, and these include calculated SAD and VD. The instrument originally used was developed by Rosen (1964) and utilizes the method of dark-field microscopy, focusing diffracted light from a particular angle onto a photomultiplier tube, which converts photons to voltages. The fundamental measurement of an OPC is the scattered intensity, or voltage, from an illuminated particle. Calibrations and the OPC response function then associate these voltages with a particle size, and the number of particles above a certain size is accumulated into size bins. Light scattered by aerosol particles was originally measured at a 25° forward angle in the UW project, with the Dust instrument measuring two to four sizes from 0.15–0.3 µm. All sizes are given as radius. This was changed to a 40° angle in 1991 to allow for size resolutions between 0.3 and 10.0 µm (Hofmann and Deshler, 1991; Deshler et al., 2003), i.e., the WOPC. The switch to a laser particle counter began in 2008, measuring side-scattering in a large solid angle centered on 90° (Ward et al., 2014), i.e., the WLPC (Wyoming laser particle counter). The WOPC provided 8–12 sizes, 0.19–2.0 or 10.0 µm, while the WLPC provided eight sizes, 0.09–4.5 µm. The measurements are made from the surface to ∼30 km. Included with most measurements is a second instrument to measure all particles >0.01 µm using a condensation nuclei counter which measures particles by growing the particles to optical detection by supersaturating the air stream with ethylene glycol vapor (Rosen and Hofmann, 1977; Campbell and Deshler, 2014). The data from the Dust counter from 1971–1988 are available at https://ndacc.larc.nasa.gov/ (last access: 22 June 2025). For flights after 1988, with the Dust, and for the WOPC (1989–2013) and WLPC (2008–2020) data, see Deshler (2023). This data record is now being extended with the LASP Optical Particle Counter (LOPC) from the Laboratory for Atmospheric and Space Physics in Boulder, Colorado (Kalnajas and Deshler, 2022), which provides >50 channels spanning 0.15–10 µm. These data are also available from Deshler (2023). The instrument with the most overlap with SAGE II is the WOPC, and it will be used as the reference OPC through the rest of this paper.

The size distribution measurements are fit with a unimodal or bimodal log-normal distribution, depending on the count of channels available and which shape produces the best fit. A log-normal size distribution consists of the total number concentration N_j, the median radius μ_j, and the distribution width σ_j for each mode j. The unimodal/bimodal log-normal size distribution is given by the following integrals where a is the integration variable:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }N(r>r_{\text{ch}})=\underset{r_{\text{ch}}}{\overset{\mathrm{\infty }}{\int }}\left[\sum _j{\displaystyle \frac{\mathrm{d}{n}_j}{\mathrm{d}\ln \left(a\right)}}\right]\mathrm{d}\ln \left(a\right)\\ \text{(3)}& {\displaystyle }& {\displaystyle }=\underset{r_{\text{ch}}}{\overset{\mathrm{\infty }}{\int }}\left[\sum _j{\displaystyle \frac{N_j}{\sqrt{\mathrm{2}\mathit{\pi }}\ln \left({\mathit{\sigma }}_j\right)}}\exp \left({\displaystyle \frac{-{\ln }^{\mathrm{2}}\left(a/{\mathit{\mu }}_j\right)}{\mathrm{2}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}}\right)\right]\mathrm{d}\ln \left(a\right).\end{array}$$

However, to better account for instrument counting efficiency, this equation has been modified to reflect the instrument's ability to count aerosols at the channel boundary (Deshler et al., 2019). The equation that is used to fit the measured concentrations is now

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{N}_{\text{ch}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left[\sum _j\mathrm{d}{n}_j/\mathrm{d}\ln \left(a\right)\right]\cdot {\text{CEF}}_{\text{ch}}\left(a\right)\cdot \mathrm{d}\ln \left(a\right)\end{array}$$

where CEF is the counting efficiency of the OPC instrument, which is modeled as a cumulative Gaussian distribution:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{\text{CEF}}_{\text{ch}}\left(r\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{1}+\text{erf}\left({\displaystyle \frac{r-\mathit{\mu }}{\sqrt{\mathrm{2}}\mathit{\sigma }}}\right)\right],\end{array}$$

where erf() is the error function. In this equation, μ is the size of the 50 % counting efficiency point, and it is the size reported in the WOPC data files. The other parameter, σ, is the rate at which the instrument counting efficiency approaches its limits of 0 and 1. The previous method of fitting assumed a perfect efficiency of the instrument to count all particles above the target radius and none below. Accounting for the realistic counting efficiency of the instrument (Deshler et al., 2019) has significantly improved the agreement between extinction coefficient computed using the resulting WOPC size distributions with those measured by SAGE II. Figure 3 shows an example of the OPC channel measurements during a period of high aerosol loading and the fitted bimodal log-normal distribution for the data using the new fitting method.

Figure 3 Example of WOPC measurements (black crosses). The black lines show differential size distributions for the two modes fit during the processing of WOPC data. The cumulative bimodal log-normal distribution is shown in red. Error bars are shown accounting for the Poisson error. The legend shows the date the measurements were made, the altitude (Alt), the total particle counts for the coarse and fine modes (No1$/$No2), the median radii (ro1$/$ro2), and the widths (so1$/$so2) in each mode.

4 Variability of SAD to extinction ratio as a function of extinction ratio

Using WOPC size distributions, it is straightforward to compute SAD and VD using analytic functions. It is also straightforward to compute aerosol extinction coefficient at any wavelength using these size distributions and Eq. (1). For these calculations, we assume that aerosol is composed of sulfuric acid–water that, in turn, defines the refractive index used to compute the Mie kernels. This is usually a very appropriate assumption, but significant exceptions can occur, particularly following injections of smoke or ash into the stratosphere and, to a lesser extent, by the presence of organic aerosol or other non-absorbing, non-sulfuric acid aerosol. While several large smoke events have occurred over the past decade, they are a relatively minor component of the SAGE II aerosol record (Thomason and Knepp, 2023) and will not be considered further in this discussion, though for some specific instances and for other instruments including the ongoing SAGE III/ISS mission, the composition cannot be as easily ignored. Using only WOPC data, we compute three sets of ratios (the SAD to 1020 nm extinction coefficient ratio (SADR), the VD to 1020 nm extinction coefficient ratio (VDR), and the 525–1020 nm extinction coefficient ratio (R)) in two altitude ranges (13–19.0 and 19.5–25 km) for the period where WOPC overlaps with the SAGE II mission. This period spans the heavily volcanic Mt. Pinatubo period as well as the fairly quiescent period between 1999 and 2005.

Figure 4a and b show the comparisons of SADR vs. R for the two altitude ranges. Organizing these plots by extinction ratio makes sense as the underlying size distribution, in the tail of the size distribution that dominates the extinction calculation, must be similar to all others that produce a roughly similar R. Superficially, the distribution of data on these curves is similar, though it is clear that the lower-altitude range shows much more scattering than the higher-altitude set. The medians of SADR and VDR for both altitudes are the same for R, near 1, but diverge and maintain a difference of almost a factor of 2 for much of the range of R values.

Figure 4 Panels (a) and (b) show SAD/1020 nm aerosol extinction coefficient (AEC) to $\mathrm{525}/\mathrm{1020}$ nm aerosol extinction ratio for the low- and high-altitude groups, respectively. The dashed lines show the 20th and 80th percentile values in bins of 0.2 ratio. The solid line shows the median values. Panels (c) and (d) are the same as (a) and (b) but show VD/1020 nm extinction coefficient to $\mathrm{525}/\mathrm{1020}$ nm extinction ratio.

Figure 5 The two median radii parameters of a bimodal size distribution are shown for 13–19 km (a, c) and 19.5–25 km (b, d) as a function of the $\mathrm{525}/\mathrm{1020}$ nm extinction ratio. The median (bold blue) and 20th and 80th percentile (dashed blue) values are shown for each extinction ratio bin. Results of the swarm optimization fit are shown in red.

Figure 6 The same as Fig. 5 but for the two mode widths.

Figure 7 The same as Figs. 5 and 6 but for the concentration ratio.

The values of SADR show a very non-linear conversion between SADR and R which varies from about 1500 for R around 1 to ∼50 000 for R around 6. While there are differences in the details, the actual conversion factors for extinction coefficient to SAD are not wildly different than those from Thomason et al. (2008). The distribution of scatter around the median SAD line is clearly not Gaussian and shows a significant positive (upward in the figures) skewness. The scatter and skewness of these points demonstrate how difficult it is to infer SAD from SAGE II measurements, even with an assist from in situ observations like the WOPC. The extremes in SADR ranges are well over an order of magnitude at some values of R. At a more restricted scale, the range of the 20th and 80th percentile levels of SADR is around a factor of 4 for R between 1 and 2.5 for the lower-altitude range and decreases to around a factor of 1.7 for R around 5 in the lower-altitude analysis. In the higher-altitude analysis, the 20th and 80th percentile SADR range is notably narrower and is a factor of between 1.5–2 for values of R between 1 and 6. These ranges are larger, particularly at low values of R, than those estimated in previous analyses (e.g., Thomason et al., 2008). The range in SADR at any given value of R is almost exclusively due to variations in small particle number concentrations that are poor scatterers and thus not reflected in any SAGE II-like measurement. Thus, while it may be possible to reproduce some of the WOPC SAD behavior from SAGE-like measurements, it is clearly impossible to reproduce all of the SAD variability observed by the WOPC.

The analysis of WOPC VDR vs. R is better behaved as VD is more dependent on large particles than the lower-order moment of the size distribution, SAD. In this case, the maximum range in both altitude ranges, shown in Fig. 4c and d, are mostly less than 10 at the lower-altitude range and less than 6 for the higher-altitude range. Similarly, the range between the 20th and 80th percentile curves is smaller for VDR than SADR. We show around a factor of 2 difference for values of R less than 2.6 in the lower-altitude analysis and a factor of 1.5 difference for R between 2.8 and 5. For the upper-altitude analysis, the percentile range fluctuates between a factor of 1.1 and 1.4 for R until R exceeds 6.8, where the range increases.

5 Finding average size distributions

While the large outliers in SADR and VDR, particularly in the lower-altitude range, cannot be captured using these extinction coefficient measurements, there is still some ability to capture median behavior; therefore, there is some utility to associating SAGE II-like aerosol extinction measurements with WOPC-compatible size distributions. Therefore, we pursue the development of representative bimodal log-normal size distributions for the WOPC as a function of the inferred 525–1020 nm aerosol extinction ratio. We used only WOPC measurements for which a bimodal log-normal distribution is most appropriate and exclude those for which only a unimodal size distribution is applicable. Data are analyzed in the two altitude ranges used above, 13–19 and 19.5–25 km, reflecting the observed differences in the distribution of inferred SAD / 1020 nm ratio for those altitude ranges. We further subdivide the data by inferred 525–1020 nm extinction coefficient ratio into bins of 0.2 width for extinction coefficient ratios from 1.0 to 8.2.

To simplify from six parameters, we scale all size distributions by the total particle number so that there are five free parameters to retrieve for bimodal distributions: the fraction of the data in the first (small) mode and the width and median radius of both the first and second modes. To retrieve these parameters, we employ a particle swarm optimization algorithm (Hu and Eberhart, 2002), where the many individual sets of WOPC data are referred to as “particles” in the algorithm's nomenclature, and the particles “fly” through parameter space, trying to minimize or maximize some function. This approach is unique in that it requires only the objective (or in our case minimization) function, and it is not dependent on gradients or derivatives of this function. This makes it fairly simple on the mathematical complexity scale of retrieval algorithms while, as we find, providing robust solutions. We define the objective function, OF, to be the sum of errors for each of the unknown parameters or

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }\text{OF}=r_{\mathrm{01}\mathit{\_}\text{err}}+s_{\mathrm{01}\mathit{\_}\text{err}}+r_{\mathrm{02}\mathit{\_}\text{err}}+s_{\mathrm{02}\mathit{\_}\text{err}}+f_{\text{err}}+R_{\text{err}}\cdot w,\end{array}$$

where $r_{\mathrm{01}\mathit{\_}\text{err}}/r_{\mathrm{02}\mathit{\_}\text{err}}$ is the two median radii errors, $s_{\mathrm{01}\mathit{\_}\text{err}}/s_{\mathrm{02}\mathit{\_}\text{err}}$ is the errors for the widths of the modes, f_err is the error for the ratio of the concentration of the first mode to the total concentration, R_err is the extinction coefficient ratio error, and w is a weight. The value of w is selected to prioritize the target extinction ratio (the bin center) among possible solutions since we value this outcome for this exercise. The error values are all positive as they are based on the absolute value of the difference between the particle's parameter value and the parameter median, divided by the parameter standard deviation, defined here as

$$\begin{array}{ll}{\displaystyle }\text{parameter error}=& {\displaystyle }\left|\right(\text{particle parameter value}\\ {\displaystyle }& {\displaystyle }-\text{parameter median}\left)\right|/\\ \text{(7)}& {\displaystyle }& {\displaystyle }\text{parameter standard deviation},\end{array}$$

where the median and standard deviation values for a parameter are determined from all the values within a particular extinction ratio bin, and the “particle parameter value” is a particle's value for that particular parameter in the parameter space. The extinction ratio error (R_err) is simply the current particle's calculated extinction ratio value minus the center of the target extinction ratio bin. In this method, many individual sets of WOPC data are effectively used to explore the parameter space through a series of iterations, t, to find a global minimum in the OF for the collection of particles in each target extinction ratio bin. The iterative process is given by

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }\text{OF}\left[t+\mathrm{1}\right]=\text{OF}\left[t\right]+v\left[t\right],\end{array}$$

where v[t] is a velocity parameter governed by an attraction to an individual particle's best value and the best value found among all of the particles in the bin. Where the best value is defined as the position in parameter space which results in the greatest minimization of the objective function, the velocity is given by

$$\begin{array}{}\text{(9)}& {\displaystyle }{\displaystyle }v\left[t+\mathrm{1}\right]=bv\left[t\right]+dv\left[t\right],\end{array}$$

where

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }dv\left[t\right]=w_{\mathrm{1}}c\left[p_{\text{best}}\left(t\right)-p\left(t\right)\right]+w_{\mathrm{2}}c\left[g_{\text{best}}\left(t\right)-p\left(t\right)\right].\end{array}$$

b is a velocity damping factor, w₁ and w₂ are weights, c is a uniform random deviate, p is the current parameter set, p_best is the parameter set yielding the lowest value for OF found for the particle among all previous iterations, and g_best (or global best) is the parameter set within p_best that yielded the minimum value of OF for any particle at any iteration. The variable b controls roughly how quickly the solution moves in its current direction, while the random perturbation created by the use of c influences how strongly the solution can change directions or explore the solution space to reduce the OF value for the best individual position and the best overall solution. The weights w₁ and w₂ effectively control the degree to which the solution search can explore the full space (w₁>w₂) or push more directly toward the current consensus best or swarm solution (w₂>w₁). In swarm optimization, this is referred to as weighting between exploration vs. exploitation. For this, we have chosen the weights to slightly emphasize the attraction of the particles toward the global best, prioritizing exploitation of the particles. The v[0] values are random perturbations of the initial variable values which are driven toward more instructive values by subsequent iterations. Obviously, there are a number of empirical knobs to turn if a solution is not found easily. This generally depends on the character of the data, the variability, and the noise. In practice, we found that varying the values of b, w₁, and w₂ does not strongly affect the ultimate solutions; though, sometimes it affects how rapidly the function approached them.

Figure 8 The ratios of SAD (a, b) and VD (c, d) to 1020 nm extinction coefficient as a function of the $\mathrm{525}/\mathrm{1020}$ nm extinction ratio for 13–19 km (a, c) and 19.5–25 km (b, d) are shown. The continuous blue lines show the median with the dashed lines showing the 20th and 80th percentile values for all WOPC data from Fig. 4. The red line shows the dependence of the ratios of SAD and VD to 1020 nm extinction coefficient as a function of the 525–1020 nm aerosol extinction coefficient ratio for the swarm optimized fits.

Figures 5–7 show the final fit values for the lower- and upper-altitude groups for each of the five parameters. We also include the 20th and 80th percentile values for individual fits provided by the WOPC dataset and the median value (as a function of extinction ratio), comprised of 1541 WOPC values for the lower-altitude group and 1515 values for the upper-altitude group. We find, in general, that the parameters found by the swarm optimization are close to the median values, suggesting that the solution space is well-behaved. Some large deviations in the 80th percentile curve occur primarily at higher extinction ratio values where extinction is also generally smaller and subject to higher measurement noise. Parameter value ranges are mostly fairly constrained for a given parameter, apart from the low extinction ratio bins for concentration ratio and the first-mode median radius where there is a wider range of values. Generally, the second-mode parameter values show significantly greater variance with many more outliers than the first mode. Some of the spread in these parameters may reflect geophysical processes, like volcanic events, so grouping data for analysis in ways that reflect the state of the atmosphere may reduce the spread in the derived quantities. This will be pursued in the future.

The influence of the second mode falls rapidly as the extinction-coefficient ratio grows. When the extinction ratio is close to 1, the first-mode fraction is ≈0.90 at the low-altitude range and ≈0.80 at the high-altitude range; it rises above 0.95 for ratios of ∼2 and is essentially unity for ratios >4. The decrease in the importance of the second mode with increasing extinction coefficient ratio (and decreasing aerosol levels in general) is not surprising since higher ratios indicate scarcity of larger particles in the distribution. In these cases, it is difficult for the WOPC measurements to separate the two modes as aerosol particles become increasingly small.

Figure 9 Surface area density (a) and volume density (b) at 45° N as a function of time, calculated using the GloSSAC 525 and 1020 nm extinction values and the corresponding SADR/VDR values calculated from the derived WOPC size distributions.

Figure 8 shows the values for SADR and VDR derived using the best fit parameters from the swarm optimization for both altitude ranges compared to the median as well as 20th and 80th percentile lines for SADR and VDR, as a function of 525–1020 nm extinction ratio, for all the WOPC data. The agreement between the median values for SADR and VDR with those from the swarm fits is reasonably good, given the spread in the observed values and keeping in mind that they are multi-decade log plots. The range in SADR is larger in a relative sense than for VDR. This is almost exclusively due to a greater sensitivity of SAD to the variations in the number density of small particles compared to VD. As noted before, particles smaller than about 0.15 µm are ineffective scatterers and thus do not significantly impact aerosol extinction coefficient at SAGE II wavelengths. Unsurprisingly then, it is clear that this approach does not capture the range of values shown in Fig. 4, where the deviations from the median values can be very large. It is not overstating the case that it is impossible that any approach based on SAGE-like measurements could infer such large variations. Generally, this outcome shows that the derived size distributions at least produce values for both parameters that are consistent with median WOPC values, and that we have successfully derived a process by which a bimodal log-normal aerosol size distribution, consistent with the WOPC, can be assigned to all SAGE II observations. Assigning uncertainty to the swarm optimization technique would require further optimizations at each extinction ratio over the range of uncertainties in that ratio, and this is beyond the scope of the work undertaken here.

The strong dependence of SADR and VDR on extinction ratio shown in Fig. 8 is in contrast to the relatively minor dependence of the log-normal parameters on extinction ratio shown in Figs. 5–7. This is reflective of how sensitive SAD and VD are to small changes in median radii and distribution widths. Both SAD and VD have a highly non-linear dependence on median radii and distribution widths, which enter into the SAD and VD calculations through power-law and exponential relationships. Thus, seemingly small differences in median radii and width lead to large differences in SAD and VD.

Figure 10 SAD median values, binned by AEC ratio, derived from three different methods. SAGE II v7.00 SAD (black) with its reported uncertainty (black, dashed). WOPC reported SAD (blue) and the bin 80th and 20th percentiles (blue, dashed). The swarm method (red) uses the $\mathrm{525}/\mathrm{1020}$ nm extinction ratio measured by SAGE II to find the SADR (seen in Fig. 8), and then it is scaled by the measured 1020 nm extinction. Dashed red lines show the 80th and 20th percentile values for all SAGE II values calculated for a bin. For the swarm method and v7.00 SAD, measurements were used which were found in a ±5° latitude–longitude square centered on Laramie, WY, USA (41° N, 105° W).

A possible application for these derived size distributions could be in providing SAD and VD estimates as a product for the Global Space-based Stratospheric Aerosol Climatology (GloSSAC) (NASA/LARC/SD/ASDC, 2022; Thomason et al., 2018; Kovilakam et al., 2023). GloSSAC is a global, gap-free aerosol climatology for the years 1979 through 2022. Data in this analysis are primarily from SAGE II, but in the period immediately after the eruption of Mt. Pinatubo, some data are reconstructed using CLAES, HALOE, and ground-based aerosol lidar data (Thomason et al., 2018). Aerosol extinction coefficients, at 525 and 1020 nm, are provided every half kilometer, every month, for latitudes centered at −77.5 to 77.5 in 5° increments.

Using the 525 and 1020 nm extinction and altitude of each data point, SAD and VD can be calculated from the corresponding SADR and VDR bin values by multiplying by the 1020 nm extinction value. Figure 9 shows those derived SAD and VD values calculated for 45° N (midlatitudes being the most applicable to WOPC size distributions) for the period around Mt. Pinatubo through the quiescent period. Among the interesting features of these figures is the obvious abrupt increase of both SAD and VD below 19 km. This is expected based on the way we have approached the altitude dependence in this analysis and the significant differences we observe between the two altitude ranges. If one were implementing this approach as a retrieval algorithm, it would be beneficial to use more altitude groupings than the two used here.

For comparison, Fig. 10 shows SAD above Laramie, Wyoming, derived from various methods. In black, SAGE II v7.00 SAD along with its uncertainty from Thomason et al. (2008) is shown. In blue, WOPC SAD values, median, and 80th$/$20th percentiles, using only the OPC values discussed in this paper, are shown. In red, the swarm-method-derived values, median, and 80th$/$20th percentiles are shown. For SAGE II, the $\mathrm{525}/\mathrm{1020}$ nm extinction ratio was calculated from measurements in a 5° latitude–longitude box around Laramie, Wyoming. For each SAGE II measurement in this region, the swarm-derived particle size distribution parameters associated with that extinction ratio were then used to find SADR (shown in Fig. 8) and then multiplied by the SAGE II 1020 nm extinction measurement to find SAD. Compared to the v7.00 SAD, there is an improvement using the swarm method in the low altitude upper extinction ratio SAD values but a divergence for extinction ratios $<\mathrm{3}.$ At the upper altitudes, both the swarm method and v7.00 SAD overestimate the WOPC SAD at extinction ratios >5; otherwise, the three estimates are quite similar. Investigating these differences is beyond the scope of this paper. The intention of this analysis has been to provide bimodal parameter values for SAGE II measurements, whereas the intension of Thomason et al. (2008) was to provide a reasonable range of SAD values for SAGE II measurements with no realistic underlying particle size distribution.

6 Conclusions

Herein, we have used SAGE II and WOPC data to infer some of the limitations to inferring aerosol size distribution and some bulk properties solely from SAGE II and similar measurements. Based solely on WOPC measurements, we have inferred a median relationship between extinction ratio ($\mathrm{525}/\mathrm{1020}$) and the ratios of surface area and volume densities to 1020 nm extinction (SADR and VDR) that are broadly well-behaved; they also exhibit substantial positive excursions that are effectively invisible to SAGE-like measurements and thus cannot be reproduced in any quantitative way. We have derived representative bimodal log-normal size distribution parameters as a function of the $\mathrm{525}/\mathrm{1020}$ extinction ratios using WOPC data. These data, in extinction ratio bins of width 0.2, were then used in a particle swarm optimization algorithm to generate bimodal size distribution parameters as a function of extinction ratio. The swarm-derived distribution parameters and the inferred SADR and VDR values are generally very close to the median values from the WOPC data. Overall, these bimodal size distributions may be useful in further applications, but care should be exercised since they are based solely on the behavior of data collected over Laramie, Wyoming, and may not be applicable at other latitudes. For instance, given the differences between the lower stratospheric values and the higher altitudes we observe herein, it is questionable in our minds how well these relationships would work in the lower tropical stratosphere where particle formation may be occurring. Additional issues may arise if sulfuric acid aerosol is not the dominant aerosol type, such as aerosol following the Australian fires of $\mathrm{2019}/\mathrm{2020}$ which was observed throughout the Southern Hemisphere stratosphere.