Nabro volcano aerosol in the stratosphere over Georgia, South Caucasus from ground-based spectrometry of twilight sky brightness

Ground-based spectral measurements of twilight sky brightness were carried out between September 2009 and August 2011 in Georgia, South Caucasus. The algorithm which allowed to retrieve the lower stratospheric and upper tropospheric aerosol extinction profiles was developed. The Monte-Carlo technique was used to correctly represent multiple scattering in a spherical atmosphere. The estimated stratospheric aerosol optical depths at a wavelength of 780 nm were: 6× 10−3 ± 2× 10−3 (31 August 2009–29 November 2009), 2.5 × 10−3 ± 7× 10−4 (20 March 2010– 15 January 2011) and 8 × 10−3 ± 3× 10−3 (18 July 2011– 3 August 2011). The optical depth values correspond to the moderately elevated stratospheric aerosol level after the Sarychev eruption in 2009, background stratospheric aerosol layer, and the volcanically disturbed stratospheric aerosol layer after the Nabro eruption in June 2011. Reconsideration of measurements acquired soon after the Pinatubo eruption in 1991 allowed to model the phenomenon of the “second purple light”, a twilight sky brightness enhancement at large solar zenith angles (97–102 ). MonteCarlo modelling reveals that the second purple light is caused by multiple scattering in the stratospheric aerosol layer.


Introduction
Stratospheric and upper tropospheric aerosols are important atmospheric constituents. They play an essential role in climate evolution because they influence the radiative balance and the global chemistry (e.g. Hendricks et al., 1999;Ramachandran et al., 2000). Volcanic eruptions inject large amounts of sulfur dioxide into the stratosphere where it is oxidized in the presence of water into sulfuric acid that condensates into sulfuric acid droplets forming the stratospheric aerosol layer (SAL).
The techniques used for SAL remote sensing can be classified in satellite-borne, ground-based and in situ measurements. Ground-based techniques are mainly represented by lidar measurements that provide aerosol backscatter profiles or aerosol extinction profiles in the case of Raman lidars, and sunphotometer and pyrheliometer measurements that can give the total tropospheric and stratospheric optical depths. In situ measurements provide information about particles composition. Lohmann et al. (2007) carried out long-term pyrheliometer measurements over Oregon. They studied direct irradiance and the influence of volcanic aerosols on the derived trends for years before and after the eruptions of El Chichón in 1982 and Pinatubo in 1991. Herber et al. (1996) used measurements by sunphotometers to detect the spectral dependence of aerosol optical depth. The measurements indicated an increase of aerosol optical depth after the eruptions of Cerro Hudson and Mount Pinatubo in 1991. Trickl et al. (2013) presented 35 yr of stratospheric aerosol measurements by lidar at Garmisch-Partenkirchen, focusing on the long-lasting background period after 1997. They investigated the residual lower-stratospheric aerosol layer 2564 N. Mateshvili et al.: Nabro volcano aerosol in the stratosphere over Georgia formation in the absence of major eruptions, and the impact of moderate volcanic eruptions. Jégou et al. (2013) combined satellite-based limb and lidar measurements, ground-based lidar measurements and in situ balloon measurements to investigate Arctic region stratospheric aerosols from the Sarychev volcano eruption in the summer of 2009. Andersson et al. (2013) investigated composition and evolution of volcanic aerosols from the eruptions of Kasatochi, Sarychev and Eyjafjallajökull in 2008-2010 based on CARIBIC (Civil Aircraft for Regular Investigation of the atmosphere Based on an Instrument Container) observations. Aerosol particles of 0.08-2 µm aerodynamic diameters were collected. They determined that the main components of the volcanic aerosol were sulphate, ash and carbonaceous material, where the source of the latter is supposed to be lowaltitude tropospheric air that is entrained into the volcanic jet and plume. In samples collected in the volcanic cloud from Eyjafjallajökull ash and sulphate contributed approximately equal amounts to the total aerosol mass (45 %). In samples collected following Sarychev and Kasatochi, ash was a minor part of the aerosol (1-7 %) while sulphate (50-77 %) and carbon (21-43 %) were dominating. Renard et al. (2008) considered balloon-borne and satellite measurements and tried to distinguish between the liquid and solid particles from the tropopause to the middle stratosphere in different geophysical conditions. They detected mineral particles of different origin in the 22-30 km altitude range.
Here we make use of another technique, the twilight sounding method (Mateshvili et al., 2005) where groundbased measurements of twilight sky brightness are performed to retrieve both lower stratospheric and upper tropospheric aerosol extinction profiles.
It is well known that the stratospheric aerosol presence manifests itself by the so-called "purple light", a reddening of the twilight sky when the Sun is a few degrees below the horizon. This effect is especially strong after major volcanic eruptions when the aerosol load in the stratosphere increases dramatically.
Photometric measurements of the twilight sky brightness at one or more wavelengths allow the retrieval of quantitative estimates of the stratospheric aerosol loading (e.g. Volz and Goody, 1962;Shakh, 1969;Volz, 1975;Mateshvili et al., 1998Mateshvili et al., , 2005Mateshvili and Rietmeijer, 2002). During a previous work (Mateshvili et al., 2005) stratospheric aerosol extinctions after the Pinatubo eruption in 1991 were retrieved from twilight measurements by using SDISORT (Dahlback and Stamnes, 1991), a pseudospherical radiative transfer code, that is part of the package Libradtran (www.libradtran. org).
The last major event was the Mount Pinatubo eruption in 1991 which produced a large quantity of sulfur dioxide that strongly enhanced the SAL. After this eruption, the SAL decayed to its background condition within a few years (e.g. Baumann et al., 2003).
The minor eruption of Nabro occurred in Eritrea on 13 June 2011. Bourassa et al. (2012) used OSIRIS (Optical Spectrograph and Infra-Red Imaging System, based on Odine satellite) data to estimate Nabro aerosol cloud optical depths. Uchino et al. (2012) have observed by lidar an increase of stratospheric aerosol caused by the Nabro eruption in July 2011 over Japan.
In this paper we present lower stratospheric and upper tropospheric aerosol extinction profiles retrieved from groundbased spectral measurements of twilight sky brightness. The measurement dataset covers the period from September 2009 to August 2011. Radiative transfer computations in a Monte Carlo approach are used to feed the retrieval algorithm and to better understand the role of single and multiple scattering in the twilight period. We also consider measurements acquired in 1991, after the Pinatubo eruption, to analyse the phenomenon of the "second purple light", the late reddening of twilight sky that is sometimes visible after the "first purple light", when the Sun sinks deeper behind the horizon (Mateshvili et al., 2005).

The measurements
The twilight sounding method had been described in detail in (Mateshvili et al., 2005). Here we will briefly summarise the main ideas ( Fig. 1).
During twilight when the solar zenith angle (SZA) is greater than 90 • , the lower part of the atmosphere above the observer is shadowed and the upper part is sunlit. The boundary between the shadowed and sunlit parts of the atmosphere shifts with the progress of twilight leaving the lower layers of the atmosphere in the Earth's shadow. This gives us a natural possibility of atmospheric sounding.
When the Sun is below the horizon, only the scattered light can be measured by a spectrometer. The intensity of the scattered light depends on the vertical extinction profiles of the different atmospheric species, as well as on aerosol and molecular phase functions.
In this paper we describe two kindred phenomena, both caused by the presence of the stratospheric aerosol layer: the first and the second purple light. The first purple light, more or less intensive, is observed at red visible -near infrared wavelengths 15-25 min after sunset or before sunrise even in volcanically quiet periods. After strong volcanic eruptions yellow and orange colours are also observed.
The second purple light is a quite rare phenomenon which is observed only after strong volcanic eruptions. A very deep red colour of the sky persists more than an hour after sunset. The second purple light was observed after the Krakatoa eruption (The eruption of Krakatoa and subsequent phenomena. Report of the Krakatoa Committee of the Royal Society, 1888). It was also observed after the Pinatubo eruption (Mateshvili et al., 2005). To discuss the origin of both the first and the second purple light we consider here two datasets, one of which was acquired recently and anotherafter the Pinatubo eruption.
The measurements were carried out with two different instruments. A CCD camera with a grating spectrograph was operating in 2009-2011. A filter photometer was used in the period 1990-1993. Table 1 summarises the both instrument characteristics. Figure 2 shows an example of an acquired by CCD spectral image of the twilight sky. Pixels from the central part of the frame were averaged between the two red lines, corresponding to a spatial average, to get a twilight spectrum (Fig. 3a).
A wavelength calibration was achieved by attributing the observed atmospheric water and oxygen absorption bands (Fig. 3a) to the modelled ones (Fig. 3b). The pixel numbers were in this way converted to wavelengths.
The dark current was measured before every spectrum acquisition. The dark frames were subtracted from the measurement frames. Every spectrum was averaged (red rectangle area on Fig. 2) over the wavelength interval [777.5-782.5 nm]. Hereafter the dependence of the monochromatic intensity on SZA is called "twilight curve". The choice of the wavelength will be discussed in Sect. 3. The CCD image presents a twilight spectral image (685-800 nm) acquired by the spectrometer. Abscissa corresponds to wavelength dimension, ordinate -to spatial dimension. The spatial dimension is shown in degrees, because the distance to the scattering volume varies with the SZAs. Red lines encompass the spatial area where the spectral image was averaged. The red rectangle indicates the area where the spectral image was averaged to receive the monochromatic twilight sky intensity for the corresponding time. The measurements were carried out in a SZA range of 88-100 • . Figure 4 shows an example of measured light intensity in arbitrary units presented on a logarithmic scale. The viewing zenith angle (VZA) and azimuth were constant during an observation. The choice of the VZA and azimuth determine the angle at which light scatters on aerosol and air molecules towards the observer. Atmospheric aerosol particles are significantly larger than air molecules which leads to strong forward scattering. Light scattered towards the observer during a twilight event contains two components scattered respectively by aerosols and air molecules. Smaller scattering angles contain a larger aerosol component of the scattered light (Mateshvili et al., 2005, Fig. 5). This means that large VZAs and the boresight azimuth close to the solar direction create favourable conditions for aerosol detection. The VZA was  Mateshvili et al. (2005).
chosen between 45-60 • . This choice depended on the season and was imposed by the local topography. As the scattering angle is important, it is necessary to know exactly the observational azimuth related to the solar one. To determine the viewing azimuth the spectrometer was pointed to the Sun just before a sunset. At that moment the viewing azimuth was equal to the solar azimuth. The relative azimuth angle between the viewing direction and the following solar directions can be calculated from each spectrum acquisition time. Both solar azimuth and SZAs were calculated using the SPICE NAIF program package (http://naif.jpl.nasa.gov/ naif/toolkit.html). Refraction effects are discussed in Sect. 3. Twilight sky brightness decreases ∼ 10 5 times in the SZA range from 88 to 100 • . This exceeds the CCD dynamic range. Therefore, each twilight observation consists of a few ensembles of measurements acquired with different exposure times. The ensembles are clearly visible in Fig. 4, they are separated by SZA gaps which correspond to time gaps necessary to change the exposure time. The measurements performed with different exposure times were divided by the corresponding exposure time values to be reconnected. The sensitivity of the spectrometer and CCD detector has not been calibrated in an absolute way (see Sect. 3 for discussion). Hence, all measured intensities are in arbitrary units.

The retrieval algorithm
The twilight sky brightness was analysed at 780 nm (observations in September 2009-August 2011). There is no gaseous absorption at this wavelength (except in the weak Wulf ozone band, see Sect. 4 for discussion) and therefore only the aerosol and molecular extinctions have to be considered. The molecular extinction can be calculated from the US standard atmosphere (see Sect. 4 for the atmospheric model choice discussion). The molecular extinction is weak at 780 nm, so the impact of the uncertainties due to deviations of the real atmospheric pressure and temperature profiles from the climatological values is negligible. Hence, the aerosol extinction profile was the only profile to be retrieved. The measurements were processed in the following way before applying the retrieval procedure. First, the twilight curves (monochromatic intensities versus SZA) were built. Since the light intensity changes by several orders of magnitude during twilight, it is easier to operate with the logarithm of the light intensity. To amplify the dynamics of the twilight curve and to remove all constant calibration factors, the measurement vector y was presented as the first differences of logarithm of intensity vs. SZA: where I is the measured light intensity. The SZA step SZA was 0.1 • . The aerosol phase function was approximated by the Heyney-Greenstein phase function, which is characterised by a single scattering albedo and an asymmetry factor. The single scattering albedo A s and the asymmetry factor g s for Atmos. Meas. Tech., 6, 2563-2576, 2013 www.atmos-meas-tech.net/6/2563/2013/  stratospheric sulfate aerosols can be estimated for volcanically quiet and volcanically disturbed conditions using the Mie theory and aerosol size distributions from the ECSTRA model (Fussen and Bingen, 1999). The aerosol complex refractive coefficient can be estimated using the code developed by Krieger et al. (2000) and the temperature value taken from the US Standard Atmosphere. The estimates for volcanically quiet and volcanically perturbed stratospheric aerosols are A s = 1 and g s = 0.65 . . . 0.75. The tropospheric aerosol single scattering albedo A t and asymmetry factor g t were adopted from the Ankara AERONET site, which is closest to Tbilisi (http://aeronet. gsfc.nasa.gov). The aerosol optical parameters typically vary in the ranges A t = 0.72. . . 0.87, g t = 0.62. . . 0.66. The average values A t = 0.8 and g t = 0.65 were used in the retrieval procedure.
The sphericity of the atmosphere needs to be taken into account for the correct modelling of the twilight event. The best tool to model a 3-D spherical atmosphere at any SZA is a Monte Carlo simulation of the photon path. But a Monte Carlo code is too slow to be directly applied as a forward model for any retrieval procedure. To bypass this problem a single scattering forward model with multiple scattering corrections obtained from the Monte Carlo code was developed.
For this purpose the reference aerosol extinction profiles ( Fig. 5b) typical for volcanically quiet (solid line) and volcanically disturbed (dash-dotted line) conditions were constructed from SAGE II (Stratospheric Aerosol and Gas Experiment) aerosol profiles. The Monte Carlo code Siro (Oikarinen et al., 1999) was used to model a twilight curve at 780 nm with the reference aerosol extinction profile as an input. The used number of photons was 10 6 , a trade off between acceptable computing time and the precision of the modelled curve. The calculations were repeated 10 times to obtain a standard deviation. This simulation allowed to estimate the fraction of the twilight intensity caused by single scattering and by multiple scattering. Figure 5a shows the modelled averaged ratio I ms (hereinafter called multiple scattering correction factor) between the total light intensity I and the single scattering I ss component as a function of SZA for two aerosol profiles (Fig. 5b). Dashed lines show the Monte Carlo modelling standard deviations. Figure 5a also shows that single scattering plays a dominant role in the SZA range 88-94 • . However, the relative multiple scattering contribution depends on the aerosol loading. This contribution may be important (Fig. 5a, case 2). This means that the twilight light intensity cannot be presented as a single scattering with the multiple scattering corrections because the multiple scattering is not a small constant factor for the SZA range of interest.
The Eq. (1) can be re-written as: (2) Figure 6 shows the modelled first differences of log(I ss ) and first differences of log(I ms ) (Eq. 2) for the two different aerosol profiles considered in Fig. 5b. The first differences of log(I ms ) (Fig. 6) is almost zero at SZA < 93 • and less sensitive to the aerosol profile at SZAs > 93 • in compare with the I ms presented in Fig. 5a. It can be considered as a small correction with respect to the single scattering curve (Fig. 6). We can conclude that a single scattering approximation with multiple scattering corrections is sufficient when the modelled quantity is the first differences of logarithm of intensity (Eq. 1).
The atmospheric refraction should be considered when modelling twilight sky brightness. The refraction mainly affects the long grazing direct sunrays before they undergo a single scattering towards the observer. Here the refraction Optical depth (twilight, volcanic multiple scattering correction)  effects were considered only in case of single scattering where they are more important. Refracted ray paths were modelled for different tangent altitudes and different SZAs as described in Auer and Standish (2000). Look-up tables of the refracted rays for every SZA were calculated and used in the forward model.
The Sun angular size was taken into account. The Sun was considered as a disk with an angular size a = 32 that was divided into 20 horizontal strips. The centre of the Sun corresponds to the actual SZA, whereas the upper and lower points of the Sun corresponds to the SZA ± a/2. Modelled twilight curves were integrated with SZA over the solar disk. The main consequence of the integration is twilight curve smoothing.
Normally, the simulated twilight sky light intensities should be convolved with the instrument vertical field-ofview. Modelling shows that this convolution can be neglected in case of the measurement vector is presented as (Eq. 1).

The retrieval procedure
An error-weighted least-squares fitting to retrieve the aerosol extinction profiles was performed by means of a Levenberg-Marquardt algorithm.
The cost function was presented as: where x is the state vector, the function f (x) represents the forward model, x a is the a priori state vector, S a is the a priori covariance matrix, S ε is the measurement error covariance matrix and y is the measurement vector (Eq. 1). The a priori state vector with its uncertainties σ a is presented in Fig. 7. The a priori covariance matrix was constructed as follows: where σ i,j a are uncertainties presented in Fig. 7, δz = 1 km is the difference between two altitude levels and h = 3 is a scale parameter that controls the altitude correlation length. To construct the measurement error covariance matrix both uncertainties of the measurements and uncertainties of the forward model should be considered.
The following measurement errors were taken into account. CCD noise was estimated using formulas from Merline and Howell (1995).
The signal averaged over a spectral image wavelength dimension and space dimension (Fig. 2, red square) is known within a standard deviation with respect to the derived mean value. This error was included in the retrieval algorithm.
Other measurement errors are systematic errors (systematic for a particular twilight curve) such as the SZA error, The modelling shows that the systematic errors mentioned above can be neglected if the measurement vector is presented as Eq. (1).
The most important uncertainty comes from the forward model. It has three components: the random error from the Monte Carlo simulation, the systematic error introduced by the forward model simplifications described in Sect. 3 and systematic errors induced by the atmospheric model and climatological values of single scattering albedo and asymmetry factor.
To estimate the Monte Carlo error, the associated computations have been repeated 10 times with the same aerosol extinction profile. The number of photons was 10 6 . The error was calculated as a standard mean square deviation. It makes no sense to increase the number of photons to achieve a smaller value of uncertainty, because the expected natural aerosol profile variability will cause variations of multiple scattering corrections larger than the Monte Carlo simulation uncertainties.
The main simplification in the forward model is the assumption that a single scattering radiative transfer code with multiple scattering corrections calculated using a Monte Carlo code combined with an a priori aerosol profile, can be used instead of a Monte Carlo code. Multiple scattering error increases with SZA (Fig. 6). Therefore, it is important to know the SZA upper limit below which such simplification can be applied.
Aerosol extinction profiles were retrieved from a few twilight measurements acquired before and after the Nabro volcano eruption. The retrieved aerosol profiles were used as an input in the Monte Carlo code to estimate the real multiple scattering corrections. The real and a priori multiple scattering corrections coincide within the limits of uncertainty below SZA ≈ 94 • for the background aerosol case and below SZA ≈ 93.5 • after the eruption. Therefore, these SZA upper limits will be used for the retrieval. The question arises if local hazes, fogs and local and remote clouds influence the measured twilight sky brightness. We should mention that all twilight observations were carried out in visibly clear sky conditions. Twilight measurements from SZA = 90 • to SZA = 94 • take about 24 min and therefore the bottom layer aerosol optical depth can be assumed relatively stable.
The modelling was performed using the fully spherical single scattering model with multiple scattering corrections described above. The first differences of logarithm of the modelled light intensity were calculated (see Eq. 1).
The following four scenarios were considered (see Mateshvili et al., 2012 for details): 1. an aerosol layer near the surface above the observer's position which was persistent during an entire observation period, 2. an aerosol layer near the surface above the observer's position which was persistent within SZA range 92-93 • , 3. a cloud placed at 9-11 km altitude and 20-30 km distance from the observer, 4. a cloud placed at 9-11 km altitude and 30-40 km distance from the observer.
The result of the first scenario is evident: a stable aerosol layer near the observer does not disturb the measurement vector (Eq. 1). It contributes only to absorption and thus decreases the modelled light intensities by a constant factor. This factor is removed by taking the first differences of intensity logarithm.
In the second scenario the presence of the haze increases the level of noise of the measurement vector. If the haze optical depth is low, the observation is still acceptable.
In the third scenario the cloud is situated close to the air volume responsible for the primary scattering. The cloud disturbs both twilight curve and its first differences of logarithm of intensity. Such observations can be easily eliminated because the cloud is quite close to the observer's position and should be visible above the observer's horizon.
In the fourth scenario the cloud is placed near the terminator. The grazing rays cannot penetrate through the lower atmospheric layers which are already opaque due to Rayleigh scattering and aerosol scattering and absorption. The presence of the cloud does not change the measurement vector.
We conclude that it is possible to neglect a remote cloud presence if it is not visible from the observational place.
Forward model parameters errors are not included in the retrieval algorithm. and 50th percentiles correspondingly for the ensembles of retrieved aerosol extinction profiles (black), OSIRIS aerosol extinction profiles (blue) and aerosol extinction profiles retrieved from twilight measurements, using a new multiple scattering correction factor (red). The stratospheric optical depth indicated above the figure was estimated using the corrected aerosol profiles.
In the forward model we used the US standard atmosphere. The error due to the deviations of the real atmospheric parameters from the standard ones was estimated by comparing the retrieved aerosol profiles for three cases: (1) the standard atmosphere, (2) mid-latitude summer atmosphere and (3) mid-latitude winter atmosphere. The differences between the three retrieved aerosol profiles are well within the standard deviation of the retrieval.
To estimate the errors caused by deviations of the asymmetry factor and the single scattering albedo from the climatological values (see Sect. 3) the aerosol extinction profiles were retrieved from the same observations using the extreme values of the optical parameters. The estimated uncertainties appeared to be strongly dependent on the particular aerosol extinction profiles. Therefore, the uncertainties were calculated for the averaged stratospheric optical depths (see Sect. 4, Eq. 9). Stratospheric aerosol asymmetry factors g s = 0.65; 0.7; 0.75, tropospheric aerosol single scattering albedos A t = 0.72; 0.8; 0.87 and tropospheric aerosol asymmetry factors g t = 0.62; 0.65 were used. The corresponding uncertainties were: 4-27 %, 4-10 %, 1-13%.
We cannot determine from our measurements the composition of aerosol; therefore we use some climatological data to estimate aerosol optical parameters. To estimate the uncertainties of our results due to the presence of mineral dust at high altitudes we assumed a total single scattering  albedo of 0.99 (Pueschel et al., 1992) for aerosol in the lower stratosphere and 0.93 (the single scattering albedo roughly estimated using results of Petzold et al., 1999) in the upper troposphere. The corresponding optical depth uncertainties were 2-13 %. We can conclude that the systematic uncertainties caused by aerosol optical parameters uncertainties vary significantly and should be estimated for a particular dataset.
All uncertainties discussed above are uncertainties of the intensity whereas in the retrieval algorithm the first differences of logarithm of the intensity are used as a forward model. This means that the uncertainties should be modified before to be included in the measurements covariance matrix S ε . The transformed signal is presented by Eq. (1). In this case, the mean squared standard deviation σ is: The retrieved aerosol extinction profile uncertainties were estimated from the retrieval Jacobian. The retrieval covariance matrix S is connected with the Jacobian K as follows (Rodgers, 2000): The diagonal elements of the covariance matrix represent the squares of the uncertainties for the corresponding extinction values.
The averaging kernel matrix is given by: The rows of A show the resolution of the retrieved profile, the columns give the response of the system to a delta function at the corresponding altitude. The number of degrees of freedom for signal d s was estimated as a sum of the diagonal elements of A.
The variance of the individual observation stratospheric optical depth was estimated as where S ij are the retrieval covariance matrix elements and n is the number of altitude points. The weighted average optical depth over a set of observations: where τ i are stratospheric optical depths with variances σ i calculated from individual observations, N number of observations in the set. The variance of the average is: 5 Results and discussion

Aerosol extinction profiles
Figures 8-10 present three typical examples of measured twilight curves fitted by the modelled ones (a), the retrieved aerosol extinction profiles (b) and corresponding averaging kernels (c, d). Figures 8d-10d show the rows of the averaging kernel matrix. The picks show maximums of sensitivity on the appropriate altitudes. The vertical resolutions of the retrieved aerosol extinction profiles can be estimated as halfwidths of the picks. The experimental and the modelled twilight curves (a) have been shifted to coincide at SZA = 92 • . This is necessary because in the retrieval algorithm it is not the twilight curves, but their first differences of logarithm of intensities that are used. Figure 8 corresponds to the background aerosol loading, Figs. 9 and 10 to the volcanically disturbed aerosol layer observed after the Nabro eruption. The blue line on Panel (a), Fig. 9 shows also a background twilight curve (24 June 2011) for a comparison with the volcanically disturbed one. For background stratospheric aerosol loading conditions the uncertainties of the retrieved aerosol extinction are high (Fig. 8b). The averaging kernel matrix (Fig. 8c) shows that the vertical resolution was about 2 km in the upper troposphere and about 6 km in the stratosphere (Fig. 8d).
Uncertainties decrease at the layer altitude for the enhanced stratospheric aerosol layer after the Nabro eruption (Figs. 9-10b), which simply means that aerosol influences the signal more strongly. The averaging kernel matrix ( Fig. 9-10c, d) show that the vertical resolution is very nonuniform and increases up to 1.5-2 km just below the layer. Aerosol extinctions below 5 km are strongly affected by the upper layers and are not reliable.
To investigate the retrieved aerosol profiles variability, for each altitude the 15.87th, 50th and 84.13th percentiles were calculated from an ensemble of retrievals. The 50th percentile is the median value; the 15.87th and 84.13th percentiles represent the distribution width. The same percentiles from an ensemble of OSIRIS aerosol extinction profiles (http://osirus.usask.ca) are presented by blue curves. The dataset was divided in four parts: fall 2009-winter 2010 (Fig. 11a), spring-summer 2010 (Fig. 11b), fall 2010-winter 2011 (Fig. 11c) and summer 2011 (Fig. 11d). Time periods of the OSIRIS data acquisition correspond to the periods of twilight measurements indicated above the panels (Fig. 11a, b, d). In the case of Fig. 11c time period of the OSIRIS data acquisition was September-November 2010. The region of the OSIRIS data acquisition was 40-50 • E, 40-44 • N. Figure 11b, c correspond to the background stratospheric aerosol conditions, and Fig. 11a, d correspond to the volcanically disturbed SAL after the Sarychev and Nabro eruptions.
The post-Nabro enhanced SAL persisted at about 17 km altitude during July-beginning of August 2011. The observed altitude of the layer is in a good agreement with the lidar measurements ( Table 2). The enhanced stratospheric aerosol layer was observed just above the tropopause. The tropopause altitude for Tbilisi in July 2011 was 16 km. It was derived from the real temperature and pressure profile extracted from ECMWF data (The European Centre for Medium-Range Weather Forecast, http://www.ecmwf.int/ products/data/archive/).
Above we have shown that in the volcanically quiet conditions the uncertainties of the retrieved individual aerosol profiles are large. Moreover, averaging kernels show that in these conditions the vertical resolution in the lower stratosphere is low. Therefore, it is better to consider stratospheric aerosol optical depths instead of aerosol extinction profiles. In our paper (Table 3) we compare aerosol optical depths obtained from twilight measurements both for volcanically quiet and disturbed periods with OSIRIS stratospheric aerosol optical depths (Bourassa et al., 2012).
We believe that there is no sense to validate the aerosol profiles retrieved in volcanically quiet period, only corresponding stratospheric aerosol optical depths. Nevertheless in Fig. 11b, c we present such a comparison. Figure 11a and b show a good agreement with our results. Aerosol profiles presented on Fig. 11a are attributed to the period of the SAL decay after the Sarychev Peak (48 • N, 153 • E) eruption in the middle of June 2009 (Mattis et al., 2010, http://www.nasa.gov/multimedia/imagegallery/ image_feature_1397.html) when there was a moderately enhanced aerosol loading in the stratosphere and therefore the retrieved stratospheric aerosol extinction uncertainty was lower. Profiles presented on Fig. 11b, c correspond to the volcanically quiet period. Discrepancies with OSIRIS profiles are larger on Fig. 11c than on Fig. 11a, but averaged stratospheric aerosol optical depths are in good agreement (Table 3). Figure 11d shows aerosol profiles with enhanced SAL after the Nabro eruption. There are no OSIRIS measurements in the region 40-50 • E, 40-44 • N during the considered period. Therefore we had to take larger region 0-100 • E, 35-44 • N, i.e. consider OSIRIS measurements which were acquired quite far from our observational place. OSIRIS aerosol profiles show smaller aerosol extinctions in the SAL than twilight profiles. Closer inspection shows that many OSIRIS aerosol extinction profiles do not extend below 18-19 km whereas the maximal extinctions often were lower (see Table 2). This is one reason why the retrieved aerosol extinction profiles show larger extinctions of the SAL. The other reason can be an underestimation of multiple scattering in the aerosol layer by our retrieval code. Indeed, the multiple scattering corrections used in the retrieval algorithm were estimated for the background aerosol layer. Volcanic multiple scattering correction was estimated by Monte Carlo code using the aerosol profile that represents 50th percentile for the ensemble of retrieved aerosol extinction profiles (Fig. 11d, black solid line). The post-Nabro aerosol extinction profiles were retrieved using the new correction. This led to smaller aerosol extinctions in the SAL (see Fig. 12, red line). Still there is about 50 % systematic shift between OSIRIS and our 50th percentile (Fig. 12, red line) aerosol extinction profile.
We tried also to compare our extinction profiles with GOMOS ones, but could not find GOMOS measurements close enough to ours. Nevertheless, GOMOS aerosol extinction profiles (see Fig. 13, the profiles were retrieved in the frame of the AERGOM project) show that the Nabro aerosol cloud was very dense, with extinctions up to 2 × 10 −3 . . . 3 × 10 −3 km −1 at 17 km altitude which is consistent with our results.
To estimate stratospheric aerosol loading all retrieved aerosol profiles were integrated above 16 km. The optical depths are shown in Table 3. The uncertainties were calculated as it was described in Sect. 4. Although the tropopause altitude is lower than 16 km in winter, the same low limit of integration over altitude was taken for all observations to make the result comparable. The post-Nabro optical depth was calculated twice: from aerosol profiles retrieved using background multiple scattering correction and ones retrieved using volcanic multiple scattering correction. The stratospheric aerosol loading was 2-3 times higher after the eruption than in background conditions. Bourassa et al. (2012) presented the weekly zonal mean stratospheric aerosol optical depths at 750 nm estimated from OSIRIS measurements. The stratospheric aerosol optical depth varied in the range 8 × 10 −3 . . . 10 × 10 −3 at Tbilisi latitude after the Nabro eruption which is very close to our result. Their background aerosol level (2 × 10 −3 . . . 4 × 10 −3 ) also agrees with our estimates. The relatively elevated stratospheric aerosol optical depth level 6.1 × 10 −3 is in agreement with the OSIRIS results (6 × 10 −3 . . . 8 × 10 −3 ) and is caused by the decaying aerosols resulting from Sarychev eruption.

The second purple light
As we have seen above, stratospheric aerosol disturbs significantly twilight sky brightness in the SZA range 91-95 • . It is interesting to model the disturbance of the twilight sky brightness in presence of volcanic aerosol at higher SZAs.
The phenomenon of the second purple light observed after the Pinatubo eruption was discussed in (Mateshvili et al., 2005) but we failed to model it. More specifically this phenomenon consists of a twilight sky brightness enhancement at large SZAs of 97-102 • (Fig. 14, red dashed line). This phenomenon is reported to be observed also after the Krakatoa eruption in 1889 (William Ashcroft, contemporary English painter made a lot of sketches of spectacular twilights after the eruption, see http://www. scienceandsociety.co.uk/results.asp?image=10316149\ &wwwflag=\&imagepos=8\&screenwidth=1280).
The data acquired in city conditions cannot be used for this purpose because the measurements at SZAs larger than 95 • were strongly affected by light pollution. To understand how stratospheric aerosols can affect twilight sky brightness at larger SZAs, we consider the measurements acquired in rural conditions. The Monte Carlo code Siro (Oikarinen et al., 1999) was used to model twilight sky intensities at large SZAs. Figure 14 shows the measured pre-Pinatubo (black dashed line) and post-Pinatubo (8 November 1991, red dashed line) twilight curves. The measured post-Pinatubo twilight curve is compared with three curves modelled using different aerosol profiles presented on Fig. 15. Figure 15 shows the retrieved from the measured twilight curve (Fig. 14, red dashed line) aerosol profile (black solid line), modifications of the retrieved aerosol profiles (solid red and dash-dotted blue) and SAGE II aerosol profiles (green lines). The coincidence criteria between twilight and SAGE II observations were one day and 1000 km.
The black solid line (Fig. 14) is modelled using the retrieved aerosol profile (Fig. 15, black solid line). There is no second purple light signature on the modelled curve; it is very close to the twilight curve measured in the background conditions (black dashed). It is clearly necessary to derive conditions in which the second purple light will be visible. Therefore, a new aerosol profile was constructed. The blue dot-dashed curve (Fig. 14) was modelled using modified aerosol profile (Fig. 15, blue dot-dashed line shows the modification). The modelled twilight curve is in good agreement with the measurements at SZAs smaller than 97 • and the corresponding aerosol profile does not contradict the SAGE II profiles. At larger SZAs the twilight sky brightnesses are higher than ones calculated using the unmodified aerosol profile, showing weak second purple light signature.
The red solid line (Fig. 14) is modelled using strongly modified retrieved aerosol profile (Fig. 15, red solid line shows the modification). This example demonstrates that the second purple light is caused by multiple scattering in the stratospheric aerosol layer. The second purple light is prominent, but the modelled curve does not fit well the measured one at small SZAs. Figure 16 shows the same modelled twilight curve as Fig. 14 (solid red lines on both Figs) and the contributions of the different orders of scattering to the first and the second purple lights.
A limitation of our model seems to be the reason why the same aerosol profile does not fit both the first and the second purple lights. In our model we suppose that the aerosol extinction profile is fully spherically symmetrical. This limitation seems not to be important in case of the aerosol extinction profile retrieval from the first purple light. The reason is that the first purple light is mainly caused by primary scattering in the SAL near the observational place. The second purple light is caused mainly by double and triple scattering (Fig. 16). When the aerosol extinction below the stratospheric aerosol layer is small the light primarily scattered in the SAL far from the observer's place can pass below the layer and scatter again in the SAL near the observer's place. In this case the spherical dissymmetry of the aerosol layer becomes important and the observed twilight sky brightnesses cannot be modelled correctly. We use the artificially constructed aerosol profile (Fig. 15, red solid line) to show that indeed the second purple light (Fig. 14, red solid line) is caused by the aerosol distribution in the lower stratosphere and not by presence of aerosol in the mesosphere.

Summary and conclusions
In this paper we have considered the stratospheric aerosol extinction profile retrieval procedure from twilight sky light intensities. Ground-based spectral measurements of twilight sky brightness at different SZAs were carried out in Tbilisi, Georgia, using the CCD spectrometer based on the SBIG CCD camera ST9E. The Monte Carlo code Siro (Oikarinen et al., 1999) was used to model twilight sky brightness variations as function of SZA. It was shown that for the selected wavelength of 780 nm and the considered SZA range of 90-94 • the single scattering plays a crucial role and the multiple scattering can be considered as a relatively small correction. Since the absolute calibration for a particular wavelength is a constant factor, it can be easily excluded by taking first differences of logarithm of intensity from both the data and the model. The Levenberg-Marquardt algorithm was applied to retrieve aerosol extinction profiles from the twilight sky brightness as function of SZA.
The retrieval errors depend on stratospheric aerosol loading. The volcanically disturbed SAL extinctions can be estimated with 50-60 % error, whereas the errors of the individual background stratospheric aerosol extinction profiles vary in 200-500 % range.
The retrieved aerosol profiles were grouped into volcanically quiet and disturbed time intervals. The 15.87th, 50th and 84.13th percentiles were calculated from an ensemble of retrievals for each altitude. The same percentiles (Fig. 11) were calculated from an ensemble of OSIRIS instrument onboard the Odin satellite aerosol profiles (http://osirus.usask. ca). The agreement between retrieved aerosol profiles and the OSIRIS ones were better in the volcanically disturbed case.
Due to high uncertainties of the background aerosol profiles, we prefer to consider stratospheric aerosol optical depths averaged over a set of observations instead of an individual extinction profiles. Our estimates of the stratospheric aerosol optical depths are shown in the Table 3. They correspond to the moderately elevated SAL level after the Sarychev eruption in 2009, background SAL, and the volcanically disturbed SAL after the Nabro eruption in June 2011. Bourassa et al. (2012) estimated the stratospheric aerosol optical depth in the same periods. These values are close to our estimates.
The averaging kernels show that the vertical resolution in the stratosphere depends also on aerosol loading. In background conditions stratospheric aerosol vertical resolution value is about 6 km, and decreases to 2 km in presence of volcanic aerosol.
The detected Nabro SAL altitude was 17 km. The layer altitude is in good agreement with lidar measurements (Table 2).
A few months after the Pinatubo eruption in November 1991 a "second purple light" phenomenon was observed (Mateshvili et al., 2005). The retrieved aerosol profile was used to model the "second purple light", an enhancement of twilight sky brightness at large SZAs of 97-102 • (Fig. 14). The modelling results show that the "second purple light" is caused by multiple scattering in the stratospheric aerosol layer. The necessary condition to observe the "second purple light" is a high transparency zone below the SAL.
We conclude that the used technique based on twilight sky brightness measurements allow to separate stratospheric and tropospheric aerosol optical depths and estimate stratospheric aerosol loading with reasonable uncertainties. The retrieved aerosol extinction profiles for volcanically perturbed periods have acceptable uncertainties. The twilight sounding method is therefore a good complement to other groundbased methods.