A new regularized inversion method for the retrieval of stratospheric aerosol size distributions applied to 16 years of SAGE II data (1984–2000): method, results and validation

Abstract. We apply a regularization method for the optical inversion of SAGE II aerosol extinction profiles and derive the particle number density N , the mode radius r and width s of an effective lognormal aerosol size distribution. The constraint applied to the inversion scheme allows us to appreciably enhance the stability of the solution. Therefore, because of the disposal of a more stable inversion scheme and of the wide extend of SAGE II data in time and space, we were able to improve the estimation of the aerosol parameter profiles with respect to previous published retrievals and, hence, our knowledge of the aerosol distribution characteristics in space and time. After presenting the inversion method and retrieved profiles concerning the particle number density profile over the time period 1984–2000, we validate our results by means of data derived from both in situ and remote spectral measurements. We also discuss the limits of the comparison between the various types of measurements due to their respective particularities. The validation gives a satisfying agreement with other data sources for N and r as long as the mode radius is not too small compared to the shortest SAGE II wavelength, whereas s appears to be less easily retrieved with a good accuracy. Key words. Atmospheric composition and structure (aerosols and particles; middle atmosphere – composition and chemistry; volcanic effects)


Introduction
Due to their role in the physico-chemistry of the atmosphere and hence in the global circulation and climate (e.g.Russell et al., 1996), stratospheric aerosols are known to be an important component of the atmosphere.Therefore, a lot of effort is devoted to the comprehension of their microphysics (Thomason, 1992;Weisenstein et al., 1997) and dynamics (Trepte and Hitchman, 1992;Trepte et al., 1993), and to the Correspondence to: C. Bingen (Christine.Bingen@oma.be)quantification of parameters able to characterize them: concentration, size, composition, etc.
The various types of measurements bring specific information and are able to refine our knowledge of aerosols on complementary ways.Space experiments like SAGE II (Chu et al., 1989) provide large-scale information (Trepte et al., 1994) and a global insight on the extinction coefficient of stratospheric aerosols (Hitchman et al., 1994;Thomason et al., 1997;Wang et al., 1995).However, they cannot discriminate particles where the size is very small compared to the wavelength, due to fundamental limitations of the scattering theory in the Rayleigh limit.Therefore, they are unable to provide reliable information about very thin particles because the absorption of UV radiation in the stratosphere prevents the use of spectral channels with very short wavelengths.On the other hand, in situ measurements (e.g.Deshler et al., 1992Deshler et al., , 1993;;Deshler, 1994;Sugita et al., 1999) and observations from ground based lidar stations (e.g.Del Guasta et al., 1994;Guzzi et al., 1999;Stein et al., 1994) are very useful complementary approaches, giving direct access to accurate information about the aerosol characteristics.The disposal of adequate impactors and particle counters allows quite a reliable description and quantification of very thin aerosol particles and condensation nuclei.But their reach is essentially local.
Despite the wide range of experimental techniques able to give a good insight into aerosol characteristics, it is still difficult to obtain a refined knowledge of the stratospheric aerosol microphysics because their composition and size vary with altitude, time and volcanic loading, due to nucleation, coagulation, condensation or evaporation, transport and sedimentation.Moreover, the aerosol size distribution retrieval by optical inversion of spectral extinction data is an ill-conditioned problem, so that retrieved aerosol parameters like the particle number density, mode radius and dispersion of the particle size distribution are found to be strongly correlated (Echle et al., 1998).Therefore, most of the published works report only a few estimates of those parameters, or else make use of integrated parameters (effective radius, surface area density, volume density) less sensitive to the ill-posedness of the inversion problem.For instance, Thomason et al. (1997) derived surface area densities from SAGE II extinction profiles, using a simplified relationship between both quantities, and proposed a global climatology of stratospheric aerosols based on this parameter.Brogniez and Lenoble (1988) use SAGE II events at various altitudes to infer the size distribution characteristics, using effective radius and variance and assuming the aerosol size distribution to be monomodal and lognormal.The comparison of their results with values derived from ground-based lidar and balloonborne polarimeter measurements (Brogniez et al., 1992) reveals the difficulty of comparing aerosol parameters derived from differents experimental techniques and the possible discordance between retrieved aerosol parameters, even when corresponding radiative quantities (extinction, optical thickness, reflectance or polarization ratio) are found to be in good agreement.Later, Anderson and Saxena (1996) used a modified randomized minimization search technique for inferring the surface area density, effective radius of aerosols, total columnar number density and mass loading and presented profiles spanning latitudes between 30 • N and 60 • N during the period March 1991 -March 1994.In order to improve the stability of the solution, the uncertainty of the extinction is reduced by assuming that the wavelength dependence of the extinction follows the Angström's law.Finally, Hervig and Deshler (2002) chose the surface area density and the extinction (reconstructed from the size distribution in the case of in situ measurements) for a comparison of aerosol measurements made by SAGE II, HALOE and balloon-borne optical particle counters.
The purpose of this work is to contribute to the improvement of the optical inversion techniques and to a better quantification of aerosol parameters, more particularly for non integrated quantities such as the mode radius and width, and the particle number density.Therefore, we applied an original scheme to the SAGE II data, based on a regularized inversion method.After a description of the algorithm, we present and validate our results using various data sets, coming from another space experiment and from in situ measurements as well.

Inversion method
Light scattering and absorption by a population of atmospheric spherical particles leads to an extinction coefficient β given by where N (z) is the particle total number density at altitude z, f (r) is the size distribution of particles of radius r and Q(r; λ) is the extinction cross section at wavelength λ, calculated from the Mie scattering theory (van de Hulst, 1957).
It is generally accepted that stratospheric aerosols consist of a mixture with a weight fraction of about 75% H 2 SO 4 and 25% H 2 O, characterized by a refractive index of 1.43 in the visible wavelength range.The particle size distribution is usually expressed as a lognormal function with mode radius ρ and width σ .Multimodal distributions have not been considered here due to the limited spectral range, the reduced number of SAGE II channels (4 channels at λ = 1.020, 0.525, 0.453 and 0.385 µm) and the measurement error level.
Using experimental values β e (z; λ) of the extinction coefficient, the purpose of the optical inversion is to retrieve the aerosol parameter profiles N (z), ρ(z) and σ (z) from Eqs. ( 1) and ( 2).The problem can be simplified by elimination of N(z) through the use of the normalized extinction profiles: where we used the 1.020 µm SAGE II channel as reference wavelength λ r .From the Mie theory and the Eqs.( 1) and ( 2), we can compile a comprehensive table of theoretical values of the normalized extinction from a set of ρ and σ values.
The inversion strategy consists of searching, for each altitude z, values of both parameters which minimize the squared difference between the model β t n (ρ, σ ; λ) and the normalized experimental data β e n (z; λ) obtained by applying Eq. ( 3) on β e (z; λ).Minimizing the merit function leads to the unconstrained solution ρ 0 (z), σ 0 (z) .The λ i refers to the four SAGE II wavelengths and β e n (z; λ i ) is the estimated experimental error on β e n .As already mentioned by several authors (Echle et al., 1998;Fussen et al., 2001a), the minimization of Eq. ( 5) is a numerically delicate problem due to the partial interdependence of the parameters to be retrieved and to a rather weak sensitivity of the model extinction to the wavelength.In particular, the study of the merit function shows that a very flat topology for M may be encountered in the neighbourhood of the solution, leading to large fluctuations in the result, especially for σ .
In order to attenuate the ill-posedness of the problem, Fussen et al. (2001a) recently proposed a vertical regularization method constraining the {ρ(z), σ (z)} to vary continuously with z.Here, we adapt this inversion method in order to enhance the numerical stability of the optimization.We propose to express the altitude dependence of {ρ(z), σ (z)} as the positive quantities where the maximum values ρ, σ are bounded within the working ranges [ρ min , ρ max ] and [σ min , σ max ], respectively.The number n of fit coefficients a i , b i was chosen to equal 4, as a good trade-off between reduction of the vertical resolution, stability of the solution and computing time.
The constrained problem comes down to finding the peak values ρ, σ and the a i , b i coefficients for which the merit function summed over the whole altitude range and computed using Eqs.( 6) and ( 7), is minimal.The minimization of Eq. ( 8) was achieved using a Levenberg-Marquardt algorithm.In order to check the stability of the solution of the constrained problem, the nonlinear optimization was also processed using the simplex search method and the BFGS Quasi-Newton method (Press et al., 1992).We compared the performances of the three methods and the Levenberg-Marquardt algorithm turned out to be the most satisfactory.The inversion scheme was then applied to all SAGE II events.Usually, the iterative process converged rapidly and the obtained solutions were robust with respect to the first-guess choice.
From the knowledge of ρ(z) and σ (z), we retrieved the number density N(z) by using Eqs.( 1), ( 3) and (4).Binned profiles were then built from all ρ(z), σ (z) and N (z) profiles, using time and latitude intervals of, respectively, one month and 10 degrees.

Basic data set
The SAGE II experiment, launched in October 1984 and still operational, has been described extensively in the literature (e.g.Chu et al., 1989).We have used here the version 6.0 extinction profiles available from October 1984 up to March 2000.In this version, extinction profiles were reprocessed using a better estimation of the altitude and a refined quantification of the extinction error 1 .Indeed, we observed that {ρ(z), σ (z), N(z)} profiles that were found to take unrealistic values at high altitudes (z > 30 km) when using the version 5.931, showed a significant improvement with version 6.0.

Data sources
In order to validate the optical inversion method, we compared the binned aerosol parameter profiles obtained from 1 No publication is available at this time about the reprocessing of SAGE II data.Information can be found in the documentation of the SAGE II data, version 6.0, or on the website of the LaRC: http://www-sage2.larc.nasa.gov/SAGE II with several other data sources found in the literature.Deshler et al. (1992Deshler et al. ( , 1993Deshler et al. ( , 1994) ) carried out several balloon-borne measurements during 1991 and 1992, using optical particle counters.They found mostly bimodal size distributions whose modes were attributed to new particle formation and to condensation on pre-existing aerosols following the Pinatubo eruption.
During the same period, Pueschel et al. (1994) observed aerosol properties from aircraft measurements using impactors.Several campaigns were planned from February 1991 to March 1992, at various latitudes between 26 • N to 89 • N. The aerosol size distribution was studied at altitudes between 16.5 and 20.2 km and between 9.5 and 12.6 km altitude, using a bimodal distribution.
More recently, Sugita et al. (1999) described measurements carried out over Kiruna (68 • N) and Aire sur l'Adour (44 • N) between September 1993 and February 1995, using a balloon-borne optical particle counter.Aerosol concentration profiles are reported between about 10 and 30 km altitude.
All these in situ measurements were found to be very useful for the validation of our results because they are based on totally different experimental techniques.However, discrepancies can be expected between our results and these local profiles, precisely due to the totally different character of the measurement.The reasons for these discrepancies have to be emphasized for a correct evaluation of the validation results.This point is discussed in the next section.
Finally, Fussen et al. (2001a) described very recently the vertical behaviour of particle size distributions computed from extinction profiles measured during the period August 1992 to May 1993 by the Occultation RAdiometer experiment (ORA) within the latitudinal range 40 • S, 40 • N. A spectral inversion of the measurements from the 340, 385, 435, 442, 600, 1013 nm channels of ORA (Fussen et al., 2001b) resulted in extinction profiles ranging from 10 to 50 km.So, the ORA experiment supplies an extended set of data in the period following the Pinatubo eruption for the validation of our results.

Discussion
Figure 1 presents the temporal evolution of the particle number density profile between 1984 and 2000 at the equator; it clearly shows a slow decay during the post-volcanic periods (especially after eruptions of El Chichon in November 1982, Ruiz in November 85 and Pinatubo in June 1991), due to both coagulation and sedimentation processes.
An estimation of the corresponding error is given in Fig. 2; this shows that the uncertainty remains within the 0 − 100% range between 18 and 30 km, and increases more significatively outside these limits.Due to the retrieval scheme of the particle number density, this error is greatly influenced by the error on the mode radius and width.The uncertainty on the corresponding mode radii is usually found to be smaller than 25% between 18 and 30 km altitude, and smaller than 35% between 12 and 18 km.Above 30 km, the error reaches 50%.1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 Year On the other side, the uncertainty on the mode width is less than 100% between 18 and 25 km.

Comparison with in situ measurements
In situ measurements are characterized by a great variability of the experimental values, due to the strong influence of local air mass motions, and also due to the difficult control of the experimental conditions.For these reasons, it is usual to find variations of the measured particle concentration that can reach one order of magnitude within an altitude interval of 1 km.As an example, Deshler (1994) observed strong variations between four successive measurements carried out over Kiruna within a delay of one month, due to the displacement of the polar vortex position and the local incursion of midlatitude air masses.Another illustration is given by the comparison between lidar measurements carried out by Stein et al. (1994) at Sodankylä (67 • N, 26 • E) on 13 February 1992, and Deshler's measurements (Deshler, 1994) on the same date above Kiruna (68 • N, 21 • E).Despite the proximity of these observations, both data sets show discrepancies up to about 20% on the mode radius, and about a factor of 2 on the particle number density.As a conclusion, different local measurements can show significative discrepancies, even when they are almost coincident.This should be kept in mind for validations using in situ data.In particular, when the profile to validate shows a systematic deviation over a restricted altitude range with respect to the in situ reference profile, it can reflect a local incursion of air masses from different origins.
A comparison of our results with size distribution parameters obtained by Deshler et al. (1992Deshler et al. ( , 1993) ) and Pueschel et al. (1994) is presented in Fig. 3, and corresponding numerical values are reported in Table 1.
In order to make the comparison with the SAGE II binned profiles valid, we used a mean value of the available data when they corresponded to the same bin and the same value of the SAGE II altitude grid.It is important to notice that all these reference data are situated within a rather restricted time and latitude range, so we do not expect a great variation of the mode radius and width through our corresponding data set.In most cases, the cited authors use a bimodal description.The fine mode involves very thin particles with respect to the SAGE II spectral channels, so that their radius cannot be discriminated using the SAGE II radiative measurements.For this reason, we only took into account the coarse mode radius, consisting of particles that can be well resolved in the SAGE II wavelength range.
Even if the particle size cannot be determined in the Rayleigh limit of the scattering theory, very small particles still affect the extinction value and, hence, the retrieval of the aerosol parameters.This is a possible explanation for the relative overestimation, below 0.2 µm, of the mode radius with respect to values derived from in situ measurements.For larger particles, the different estimations of the mode radius are in agreement.Taking into account the expected variability on the particle number density in the case of in situ data, the correlation between number density values is fairly satisfactory, whereas the comparison between mode width data sets shows a moderate agreement, due to the already mentioned reasons.
Another interesting way of using in situ data is to consider the vertical profiles of partial number densities, corresponding to cumulated particle classes measured by the particle counter.Using these raw profiles allows us to get rid of the bias introduced by inverting the in situ experimental data.Moreover, they allow us to validate the vertical dependence of our aerosol number density profiles.
It is easy to verify that for a lognormal particle size distribution described by Eq. ( 2), the aerosol partial number density corresponding to the number density of particles with a radius larger than a value r * , is given by Figure 4 shows a comparison of partial density profiles retrieved from our aerosol parameters and from in situ measurements published by Sugita et al. (1999), and carried out over Aire sur l'Adour in September 1993.
For all particle classes, the vertical dependence of both data sets are in good agreement.The difference between partial density values is generally less than a factor 5; it is 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 Year interesting to note that such a high factor is observed for particle classes involving very thin particles (r ≥ 0.17µm and r ≥ 0.20µm) where, we suspect, as discussed previously, an overestimation of the values obtained from remote measurements with respect to in situ data.For the other particle classes, the difference between both data sets is less important.Therefore, when taking into account the specifics of in situ measurements and the inherent inability of remote spectral measurements to distinguish very thin particles, the agreement between the equivalent particle size distributions that we retrieved from our approach and the considered reference profiles are found to be satisfactory.

Comparison with remote measurements
Finally, our results are compared in Fig. 5 with equivalent profiles obtained from ORA measurements by a similar method of retrieval after a binning of all the valid profiles over 10 • latitude intervals and 1 month time periods by means of an error weighted mean.
Within the altitude range of 18 to 40 km, the particle number density and mode radius profiles show a good agreement between ORA and SAGE II.Below this range, the increasing discrepancy between both experiments is explained by a lack of available SAGE II data.For both experiments, an increase of the mode radius above 30 km should be noted.As expected, the comparison of the mode width profiles is less satisfactory.this confirms that it is difficult to get an accurate estimation of this parameter, due to the reduced sensitivity to σ of the merit function M at fixed ρ.

Conclusions
We derived particle number density, mode radius and width of the stratospheric aerosol size distribution, by optical inversion of SAGE II extinction profiles.The profiles we obtained cover the entire latitude and altitude ranges (respectively 80 • S to 80 • N, and about 12 to 40 km) and 186 months (from October 1984 to March 2000) of the SAGE II flight duration.By using a vertical regularization method, we could partly get rid of the instabilities caused by the optical inversion problem.It allowed us to derive large scale number density profiles whereas most authors compute the surface area density that can be determined less ambiguously.Moreover, the present results were derived from version 6.0 of the SAGE II data which showed a significant improvement with respect to the previous version.
Our results were first validated by means of data coming from in situ measurements.The specifics of this kind of measurement are discussed and, taking into account their intrinsic properties, we found a good agreement between these reference data and our results.Small values of the mode radius  The particle classes take into account all the aerosol particles with radii larger than respectively 0.17, 0.20, 0.25, 0.33, 0.38, 0.44 and 0.65 µm showed an overestimate, with respect to corresponding in situ data, possibly related to the insensitivity of UV-visible measurements to the spectral signature of very small particles in the limit of Rayleigh scattering.Furthermore, validation of vertical profiles of partial number density gives good results.
On the other hand, the derived particle density and mode radius profiles show a good agreement with ORA profiles obtained with a similar inversion method.
Finally, due to the ill-conditioning associated with the optical inverse problem, it remains difficult to derive reliable values of the mode width from extinction measurements.
With this extended set of results, we are now able to study the microphysical properties of aerosols in the wide range of volcanic situations observed by SAGE II; also the aerosol transport.In the near future, these new data will be used to improve and extend the Extinction Coefficient for STRatospheric Aerosol (ECSTRA) climatology published previously (Bingen and Fussen, 2000;Fussen and Bingen, 1999).Furthermore, we will investigate the dynamical evolution of the particle size distribution in post-volcanic conditions in order to estimate the relative contributions of sedimentation, transport and coagulation, the main drivers of the stratospheric aerosol relaxation after important eruptions.

Fig. 1 .
Fig. 1.Temporal and vertical evolution of log 10 (N [cm −3 ]) at the equator, derived from SAGE II extinction profiles for the period October 1984 up to March 2000.

Fig. 2 .Fig. 3 .
Fig. 2. Temporal and vertical evolution of the relative uncertainty on the particle number density illustrated in Fig. 1.

Fig. 4 .
Fig. 4. Comparison of partial number densities for aerosol profiles obtained from this work (solid lines) and in situ measurements published by Sugita et al. (1999) (dots) and carried out on 29 September 1993 over Aire sur l'Adour.The particle classes take into account all the aerosol particles with radii larger than respectively 0.17, 0.20, 0.25, 0.33, 0.38, 0.44 and 0.65 µm

Fig. 5 .
Fig. 5. Comparison of mean vertical profiles for aerosol parameters obtained from this work (•; error range in dark grey) and from ORA (•; error range in light grey); the average is weighted by the error found on the parameters, and concerns the time period August 1992 -May 1993 and the latitude interval [10 • N-20 • N].

Table 1 .
Comparison between aerosol parameter values from various reference sources and binned values from SAGE II (this work); when the cited observations refer to a bimodal distribution, only the coarse mode, well detected by the SAGE II wavelength range, is mentioned