Long-term (1999-2019) variability of stratospheric aerosol over Mauna Loa, Hawaii, as seen by two co-located lidars and satellite measurements

As part of the Network for the Detection of Atmospheric Composition Change (NDACC), ground-based measurements obtained from the Jet Propulsion Laboratory (JPL) stratospheric ozone lidar and the NOAA stratospheric aerosol lidar at Mauna Loa, Hawaii over the past two decades were used to investigate the impact of volcanic eruptions and pyro-cumulonimbus smoke plumes on the stratospheric aerosol load above Hawaii since 1999. Measurements at 355 nm and 532 nm conducted by these 5 two lidars revealed Ångström exponents of -1.6 for background plumes and -0.6 for a PyroCb plume recorded on September 2017. Measurements of the Nabro plume by the JPL lidar in 2011/2012 showed a lidar ratio of (64.0± 12.7) sr at 355 nm around the center of the plume. The new GloSSAC, CALIOP Level 3 and SAGE III-ISS stratospheric aerosol datasets were compared to the ground-based lidar datasets. The intercomparison revealed a generally good agreement, with vertical profiles of extinction coefficient within 50% of discrepancy between 17 km and 23 km above sea level (ASL), and 25% above 23 km 10 ASL. The stratospheric aerosol depth derived from all these datasets shows good agreement, with the largest discrepancy (20%) being observed between the new CALIOP Level 3 and the other datasets. All datasets consistently reveal a relatively quiescent period between 1999 and 2005, followed by an active period of multiple eruptions (e.g., Nabro) until early 2012. Another quiescent period, with slightly higher aerosol background, lasted until mid-2017, when a combination of extensive wildfires and multiple volcanic eruptions caused a significant increase in stratospheric aerosol loading. This loading maximized at the 15 very end of the time period considered (fall 2019) as a result of the Raikoke eruption, the plume of which ascended to 26 km altitude in less than three months.

atmospheric temperature and ozone profiles, but has also gained relevance during recent years because of their potential use as a geoengineering tool to reduce the impact of global warming (Rasch et al., 2008). Dominated by sulfate aerosols, stratospheric aerosols are typically found in the form of a layer that extends from the tropopause up to 35 km above sea level (ASL) in the tropics, and about 30 km ASL at mid-latitudes (Hitchman et al., 1994;Kremser et al., 2016). This stratospheric aerosol layer (SAL), also known as Junge layer, was discovered in 1960 by means of balloon-borne measurements (Junge and Manson, 5 1961).
Although the stratospheric sulfur burden has been dominated by periodic volcanic injections of large amounts of SO 2 and volcanic ash (Kremser et al., 2016) during the past decades, some recent studies (Solomon et al., 2011) reported an increase in the background aerosol load during less active volcanic periods that could be attributed to other sources, including anthropogenic SO 2 emissions (Brock et al., 1995;Rollins et al., 2018) and large wildfire-driven thunderstorms (Peterson et al.,10 2018).
Since its discovery, several techniques have been used to monitor the evolution of the SAL, including ground-based lidars (e.g. Fiocco and Grams, 1964;Barnes and Hofmann, 2001;Trickl et al., 2013;Khaykin et al., 2017;, balloon-borne in-situ measurements (Deshler et al., 2003), and satellite-borne lidars (Vernier et al., 2009) and spectrometers (McCormick and Veiga, 1992). While satellite-borne instruments are able to provide global coverage over several years, their 15 limited lifetime requires additional measurements that can close eventual gaps existing between different missions, and can provide a common reference to investigate possible instrumental biases that can arise from the difference in the measurement techniques and wavelengths used on each mission.
As part of the Network for the Detection of Atmospheric Composition Change (NDACC), the Jet Propulsion Laboratory (JPL) lidar group has been performing stratospheric ozone and temperature measurements at Mauna Loa, Hawaii since 1994 20 (e.g. Leblanc et al., 2006;Kirgis et al., 2013). This long-term ground-based lidar dataset provides not only stratospheric ozone and temperature records, but also a unique opportunity to evaluate current and past global stratospheric aerosol datasets. In contrast to most stratospheric aerosol systems, the JPL stratospheric ozone lidar operates in the ultraviolet spectral region (308 nm-387 nm) and has two different Raman receivers that allow the retrieval of aerosol extinction coefficient (α a ) profiles without the need to assume a lidar ratio.
The paper is organized as follows. Section 2 provides a brief description of the instruments and datasets used in this study, including a description of the JPL and NOAA lidars and the three satellite-based datasets. Section 3 describes the methods applied to the ground-based lidar datasets in order to retrieve aerosol backscatter and extinction profiles as well as aerosol optical depth, and the corrections needed in order to obtain comparable datasets. Section 4 provides an overview of the measurements conducted at Mauna Loa in the period between 1999 and 2019, and different optical properties retrieved from the synergy 5 between JPL and NOAA lidars. Section 5 presents lidar-derived optical properties. Section 6 presents an intercomparison of the two ground-based lidars with the three satellite-based datasets under evaluation. Finally, a summary of the key finding on this paper is presented in Section 7.  McDermid et al. (1995). Since then, the system underwent a few major modifications and some minor technical issues. The major modifications include the migration of the system from a mobile trailer to a building and the replacement of the Raman shifting cell used to generated the 353 nm by a Nd:YAG laser in March 2001 and the upgrade on 15 the data acquisition system in April 2019.
In its current configuration, the system transmitter consists of a Spectra Physics PIV-400 Nd:YAG laser operating at a repetition rate of 50 Hz followed by a third harmonic generator (THG), emitting pulses of about 140 mJ at 355 nm. For the generation of the "on-wavelength" needed for the ozone differential absorption lidar (ozone DIAL) retrieval, a set of two XeCl Excimer lasers is used. The first Excimer laser (Coherent LPXpro 220) is used as an oscillator while the second is used as 20 a single-pass amplifier (Lambda Physik LPX 220i). This configuration is operated at a repetition rate of 200 Hz with a pulse energy of about 300 mJ at 308 nm. Both laser beams are expanded by a factor of 5 in order to reduce beam divergence, and redirected to the atmosphere by motor-controlled mirrors used for alignment purposes.
The system receiver, mainly unchanged since the system description presented in McDermid et al. (1995), consists of a 1-m-aperture Dall-Kirkham telescope with a focal ratio of f/8. A PARC Model 192 chopper is placed between the telescope 25 and the receiver in order to block backscatter signal from altitudes below 10 km above ground level (AGL). The chopper blocking, together with the electronic gating of the photomultiplier tubes (PMTs), helps to reduce the signal induced noise (SIN) generated by the strong backscatter from low altitudes. The system receiver consists of 6 photon-counting channels: two high-intensity Rayleigh backscatter channels (355H and 308H), two low-intensity Rayleigh backscatter channels (355L and 308L) and two nitrogen Raman backscatter channels (387M and 332M). Up to April 2019, the signals were digitized 30 with Tennelec/Nucleus multichannel scaler (MCS) boards with a resolution of 300 m. After the last system upgrade, a Licel transient recorder with photon-counting and analog detection capabilities has been used. The new acquisition system has a vertical resolution of 15 m.

NOAA Mauna Loa Stratospheric Aerosol Lidar
The NOAA Stratospheric Aerosol Lidar is located in the same building as MLSOL and has been operating with its current configuration since May 1994 (Barnes and Hofmann, 2001). The lidar is based on a Spectra Physics GCR-6 Nd:YAG laser (30 Hz, 40 W at 1064 nm) with frequency doubling and tripling (532 nm and 355 nm), although the 355 nm has not been used routinely. One 61 cm telescope is dedicated to 532 nm measurements, one 61 cm to 1064 nm and a 74 cm to Raman nitrogen 5 (607 nm) and water vapor (660 nm). PMTs are used in the photon counting mode for all channels. The system uses a PC 80486 and the data acquisition electronics are MCS II boards made by Tennelec. Measurements are made during the night, usually once a week.

The Global Space-based Stratospheric Aerosol Climatology (GloSSAC)
The GloSSAC (data version 1.1) is a 38-year climatology of stratospheric aerosol properties based on measurements from the 10 SAGE instruments series through August 2005, and a combination of the Optical Spectrograph and InfraRed Imager System (OSIRIS) and CALIOP datasets thereafter (Thomason et al., 2018). The main product reported by this dataset is a series of extinction coefficient profiles at 525 nm and 1020 nm. The data corresponds to zonally averaged extinction profiles with a latitude grid of 5 • and a vertical grid spacing of 0.5 km from 5 km to 39.5 km.
In the time frame covered by this study, the instruments contributing to the GloSSAC dataset are SAGE II, OSIRIS and

CALIOP Level 3 Stratospheric Aerosol Profile
The new CALIOP level 3 (L3) stratospheric aerosol profile product (Kar et al., 2019), released in August 2018, reports monthly mean profiles of aerosol extinction, particulate backscatter, attenuated scattering ratio (SR) and stratospheric aerosol optical depth on a spatial grid of 5 • in latitude, 20 • in longitude and 900 m in altitude.
As part of this dataset, two different aerosol products are reported. One is labeled as "background" and the other is named 5 "all aerosols". While the first corresponds to profiles retrieved after removing clouds, aerosols and polar stratospheric clouds (PSCs), the second only screens out clouds and PSCs. For this study, we use the "all aerosols" data product from dataset version 1.0. Auxiliary atmospheric parameters required for the aerosol extinction coefficient retrieval, including molecular and ozone absorption corrections, are obtained from MERRA-2 (Modern-Era Retrospective analysis for Research and Applications, Version 2) (Gelaro et al., 2017). The lidar ratio assumed for the retrieval is 50 sr.

SAGE III -ISS
The SAGE III instrument (Cisewski et al., 2014), launched in February 2017 and carried by the ISS, is a moderate resolution spectrometer covering wavelengths from 290 nm to 1550 nm. Data collection is performed in three different modes, namely solar occultation, lunar occultation, and limb scatter measurement. The expected science products include vertical profiles of ozone, nitrogen dioxide and water vapor, along with multi-wavelength aerosol extinction.

15
In this study, zonally averaged extinction profiles at 521 nm for latitudes between 15°N and 25°N are used (data version 5.1). Native vertical grid spacing is 0.5 km.
3 Lidar products and data analysis 3.1 Backscatter coefficient retrieval Different approaches exist for the calculation of backscatter coefficient profiles out of ground-based lidar measurements, in-20 cluding the Klett inversion technique (Klett, 1985) and the scattering ratio approach. Generally speaking, these methods require the knowledge of a reference atmospheric density profile which is usually derived from co-located radiosonde launches or atmospheric models. In the case of MLSOL, and because this system has Rayleigh and nitrogen Raman channels with high enough signal-to-noise ratio (SNR) to cover the stratospheric aerosol layer, the scattering ratio approach is used, which can be determined based on the ratio of the Rayleigh and Raman signals according to the procedure described elsewhere (Gross et al., 25 1995;Langenbach et al., 2019).
The main advantage of using the Raman channel as the atmospheric density reference profile for MLSOL is the smaller sensitivity to timing uncertainty. Since before April 2019 the data acquisition of MLSOL had a range resolution of 300 m, the determination of the range zero bin with an accuracy better than 150 m could not be done with the traditional approach of looking for the first bin with non-zero readings. An error of 150 m (half a bin in the old MLSOL configuration) in the assumption aerosol optical depth (AOD) in the stratosphere if the Klett method is used. This is mainly due to the error introduced at the time of the so-called range correction. Although alternative methods can be applied to determine sub-bin timing in coarse resolution systems, this is only possible if access to the configuration under study is possible. Since several changes have been introduced in the systems over the last 20 years, an accurate timing characterization of old setups is not possible. In contrast, the use of the ratio of the Rayleigh and Raman channels cancel the effect of the quadratic range dependence that characterizes lidar signals 5 without the need to know the range zero bin. The only requirement is a well-known relative timing difference between the two channels. Since the Rayleigh and Raman backscatter share the same laser and acquisition timing, this difference is typically very close to zero. As a downside, and due to the limited SNR of the Raman channel, the retrieved profiles are generally noisier than the ones obtained by the Klett algorithm.
The SR retrieval from MLSOL measurements starts with the calculation of the average of the recorded signals for each 10 channel during the length of the measurement period. Typical measurement periods for this lidar correspond to two-hour measurements after astronomical twilight. After an average profile is obtained for each channel, a correction of the lidar signals for count pile-up (saturation) is applied. This is performed by assuming a nonparalyzable model and dead-times derived based on the comparison between high and low intensity channels available on each system. Then, the signal background corresponding to moon light, airglow and electronic noise is calculated as the mean of the high altitude tail of the recorded 15 signals, where no contribution of the laser scattering is expected. This background is then subtracted from each signal. Finally, a correction for Rayleigh extinction is applied to the Rayleigh and Raman channel signals to be used in the calculation of the scattering ratio (Eq. 1 and Eq. 2). (1) (2) 20 where P λ (z) is the Rayleigh-extinction corrected lidar signal, N λ,corr (z) is the saturation and background corrected lidar signal, and T m,λ (z) is the one-way Rayleigh atmospheric transmission for λ = 355 nm and 387 nm. The extinction profiles required for this correction are derived from the closest pressure and temperature profiles available from MERRA-2. While the effect of the aerosol extinction could also be included in the calculation or corrected using an iterative approach (Friberg et al., 2018), its contribution is small for the cases presented in this study and compared to other uncertainty sources.

25
The normalization is performed by dividing the signals by the average of the Rayleigh (P 355 (z ref )) and Raman signals (P 387 (z ref )) at a reference altitude range assumed to be free of aerosols (z ref = 35 − 37 km ASL in this study).
The scattering ratio is finally calculated as the ratio of the normalized and corrected Rayleigh and Raman signals (Eq. 3).
Once the scattering ratio is calculated, the aerosol backscatter coefficient (β a (z)) can be retrieved if a molecular backscatter profile (β m (z)) is assumed (Eq. 4). This reference profile can be derived, as in the case of the molecular extinction, from the MERRA-2 dataset. Uncertainty is calculated following the procedures detailed in Russell et al. (1979).

5
In the case of the NOAA lidar the Klett inversion approach is followed, and the SR is calculated based on a radiosonde and model-derived atmospheric density reference profile instead of using a Raman channel. The normalization altitude used is nearly always 38 − 40 km ASL which has improved the consistency when compared to the earlier archived profiles in the NDACC database.

10
While the Raman channel of MLSOL allows the independent retrieval of the aerosol extinction coefficient profiles and lidar ratio at 355 nm based on technique presented in Ansmann et al. (1990), the application of this method for stratospheric retrievals has been limited to relatively thick stratospheric plume due to SNR constraints. As in the case of the backscatter coefficient retrieval, the molecular reference profile required for the inversion is obtained from MERRA-2 temperature and pressure profiles.
Based on balloon-borne measurement (Jäger and Deshler, 2002), the extinction Ångström exponent relating the Rayleigh and

15
Raman wavelengths is assumed to be -1.6. Once a stratospheric extinction profile is derived, the lidar ratio of the stratospheric plume under study is calculated as the ratio of the extinction and backscatter profiles. In this work, the Raman technique is applied to measurements conducted on the Nabro plume, and the results compared to similar previous measurements (Sec. 5).
Unfortunately, because the Raman-based extinction coefficient retrieval approach is not able to provide acceptable results during most of the period under study because of the generally low aerosol load conditions, the extinction coefficient profiles 20 presented in Sec. 4 and Sec. 6 are based on the assumption of a constant lidar ratio. Since the satellite-borne datasets provide extinction profiles at either 525 nm or 532 nm, the MLSOL backscatter measurements are first translated to 532 nm using the backscatter coefficient color ratio derived from collocated NOAA lidar measurements (Sec. 5.1). Then, NOAA lidar and MLSOL backscatter profiles are converted to extinction profiles by assuming a lidar ratio of 50 sr, as typically referenced in previous studies (Trickl et al., 2013;Khaykin et al., 2017).

25
Finally, the stratospheric AOD is calculated by integrating the extinction coefficient profiles between 17 km and 33 km.
The calculation of the extinction coefficient uncertainty is performed based on a Monte Carlo analysis, adopting a 3 % standard deviation as the uncertainty on the temperature and pressure profiles from MERRA-2.   Table 1).
Three different periods defined based on the impact of the volcanic activity are also shown (see text for details). Monthly-averaged AOD retrieved from MLSOL extinction profiles (lower panel). 4.1 The first "quiescent" period (1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006) Several studies focusing on mid-latitude refer to the time period between 1997 and 2003 as a "quiescent" or background period (Trickl et al., 2013;Sakai et al., 2016;Khaykin et al., 2017). In the case of the mid-latitude station located in Tomsk, Russia, 10 a "quiescent" period between 1999 and 2006 was reported . While all definitions imply some degree of arbitrariness, the separation of the time-series in periods based on the intensity of the volcanic activity impact helps to organize the discussion. In the case of the MLSOL data series, corresponding to a tropical station, a similar background period between January 1999 and January 2006 can be defined. Although a few eruptions took place during that time frame, and with some occurring at equatorial latitudes, no strong impact in the total stratospheric AOD was observed by MLSOL.
15 Table 1. List of eruptions with VEI≥4 in the tropics and northern hemisphere between January 1999 and November 2019. Information regarding the location of the erupting volcano and the maximum plume altitude (MPA) reported by the Global Volcanism Program is included for reference (https://volcano.si.edu/). Raikoke eruption is included in the list, although the VEI has not been determined yet. The first two eruptions that occurred during this period, corresponding to the Ulawun and Shiveluch volcanos, did not show a clear impact in the extinction profiles derived from MLSOL. After October 2002, some plumes, presumably corresponding to the Ruang and Reventador eruptions, were observed by a relatively short period of time, but with small impact in the stratospheric AOD. During the background period defined for the MLSOL station, the lidar-derived AOD is 2.9 × 10 −3 ± 0.1 × 10 −3 . The variability reported corresponds to 1σ. as the result of the Nabro eruption (Sawamura et al., 2012). Between eruptions, the AOD quasi-background level reported by MLSOL is about 5 × 10 −3 ±0.4 × 10 −3 , almost 70 % higher than during the quiescent period described in the previous section.
Despite the fact that stratospheric aerosols interact with incoming solar radiation and stratospheric chemistry, and are subject to evaporation and sedimentation (Kremser et al., 2016), they proved to be a useful quasi-passive trace to investigate stratospheric circulation (Trepte et al., 1993). In particular, the diabatic plume ascent induced by residual Brewer-Dobson circulation 5 (BDC) has been studied based on several satellite observations (Trepte and Hitchman, 1992;Vernier et al., 2009;Fairlie et al., 2014). While analyzing the vertical evolution of the volcanic plumes shown in Fig. 1, a typical "tape-recorder-like" signature as in the case of water vapor observations (Mote et al., 1996) can be seen. Over the lifetime of the eruption plumes, bottom and top of the plume elevate at a similar rate. In some cases, as previously reported by Vernier et al. (2011), the space left by the plumes is filled with overshooting tropospheric clean air.

10
By analyzing the aerosol plume height evolution over time, a rough estimation of the apparent BDC ascent rate can be calculated. Figure 2 presents a close-up of the MLSOL measurements of the Sarychev plume as well as the altitude evolution of the center of the plume as a function of time. For the calculation of the center of the plume, the background extinction calculated based on the mean extinction of the six months previous to the eruption is first subtracted from the extinction time series. Since the plume has a vertical profile shape similar to a gaussian function, a function fit is used to estimate the plume 15 center. The result of this calculation, shown as black crosses in Fig. 2, allow us to estimate the ascent rate between 18 km ASL and 19 km ASL to be about 0.6 km per month (0.025 cm s −1 ). As the plume further rises, the ascent speed seems to diminish, but a quantitative estimation based on MLSOL data is difficult due to the reduced SNR. This result is in the same order of magnitude as those retrieved based on water vapor and carbon monoxide tape-recorder signatures (Minschwaner et al., 2016) and double the ascent rate derived from CALIOP measurements of the Soufriere Hills plume . For the 20 period comprehended between the end on 2009 and beginning of 2010, the ascent rate derived from water vapor records was estimated to be about 0.02 cm s −1 at 51 hPa (about 19 km ASL).

The second "quiescent" period (2013-2019) and Raikoke eruption
After the decay of the Nabro eruption plume and before the Raikoke eruption, a second relatively quiescent period can be defined based on the aerosol extinction coefficient and AOD records. Although during this period, defined between early 2013 and July 2019, some eruptions with VEI≥4 occurred, only a slight impact was observed on both, the extinction and AOD records. The observed mean AOD for the period is 4.4 × 10 −3 ± 0.7 × 10 −4 which is 50 % higher than the ones measured 5 during the first quiescent period and similar to the AOD observed in the period measured between volcanic eruptions in the volcanic active period.
In addition to volcanic eruptions, large wildfires can also contribute to an enhancement of the stratospheric aerosol load (e.g. Khaykin et al., 2018;Zuev et al., 2019). On 12 August 2017, five near-simultaneous extreme pyrocumulonimbus (pyroCb) events took place in British Columbia, Canada. According to recent studies (Peterson et al., 2018), these events injected a mass were analyzed. For all these overpasses, a thin stratospheric plume at about 15 km ASL was observed 150 km north of MLO 15 (Fig. 3). In all cases, the plumes were characterized by a low average particle depolarization ratio (< 0.1), which is compatible with smoke particles (Kim et al., 2018). The rapid equatorward transport of the plume seems in agreement with the wind field reported by the MERRA-2 reanalysis. Additionally, compatible results can be seen in Fig. 3c    Based on CALIOP L1 attenuated backscatter profiles, the back-trajectory of the Raikoke plume observed at the end of 5 September 2019 by MLSOL was estimated (Fig. 6). The tracking of this plume started with its first detection at Manua Loa on 24 September 2019 and ended on 17 July 2019 over the Kamchatka peninsula, Russia. At that point, the track is lost as there are many different plumes in the region. Over a period of slightly over two months, the plume ascended 7 km, from 19 km ASL to 26 km ASL. The tracking of the plume was conducted mainly through a manual inspection process of CALIOP L1 profiles considering the main stratospheric circulation patterns. This was possible due to the well-defined shape of the plume 10 and its strong backscattering properties when compared to other previously observed volcanic plumes.
While this well defined and elevated plume is the most prominent feature observed by MLSOL during the period between July and November 2019, an enhancement of the aerosol load between the tropopause and 21 km ASL was also observed during this period. In order to further investigate whether this plume correspond to the Raikoke eruption or the Ulawun eruption, CALIOP L3 longitudinally averaged SR cross-sections (latitude vs. altitude) from this period are presented in Fig. 7.   Since no significative stratospheric injection events were registered during the first half of 2019, the SR cross-section from June 2019 can be considered as background condition for this analysis (Fig. 7, BKG). Starting in July 2019, an enhancement of the aerosol load is clearly visible between the tropopause and 19 km ASL, with a small gap around 15°N. This gap corresponds to the division between the plume of the Ulawun (south) and Raikoke (north) eruptions. By August 2019, this gap is closed as both plumes mixed together, making them indistinguishable. Between 30°N and 40°N, a fraction of the Raikoke plume is 5 visible rising above the rest of the plume, confirming the back-trajectories presented in Fig. 6. By September 2019, the mixed Ulawun-Raikoke plume reaches 21 km ASL and increases its SR at low latitudes. The secondary Raikoke plume displaces southward reaching 20°N at an altitude of over 25 km ASL. Finally, October 2019 cross-section starts to show a decay in the SR of the plume, which is compatible with the AOD measurements by MLSOL recorded for that month (Fig. 1).

Lidar-derived optical properties 10
Most lidar-derived stratospheric aerosol long-term records consist of backscatter coefficient profiles reported at 532 nm. While some satellite-based datasets, like CALIOP, operate at the same wavelength and also provide aerosol backscatter coefficient as its natural product, other instruments like those of the SAGE series provide extinction profiles at other wavelengths (e.g.

nm).
In the latter case, the comparisons with ground-based lidar datasets require knowledge of the lidar ratio, which is typically derived using Raman lidar measurements. On the other hand, the recently launched Aeouls wind lidar mission (Stof-15 felen et al., 2005;Flamant et al., 2008) will provide aerosol backscatter and extinction profiles in the lowermost stratosphere at 355 nm. Its validation and intercomparison with other datasets (like CALIOP) will then require a good knowledge of the backscatter and extinction wavelength dependence.

The color ratio of the backscatter coefficient
Based on the collocated backscatter coefficient profiles of MLSOL, NOAA lidar and CALIOP, the backscatter color ratio (and 20 backscatter Ångström exponent) between 355 nm and 532 nm was calculated as function of altitude and time. Figure 8 presents the temporal average of the color ratio for the volcanic and second quiescent period. In the case of color ratio derived from CALIOP and MLSOL measurements, the monthly means from both datasets are used. In the case of the color ratio derived from the two ground-based lidars, only same-day measurements are included in the calculation. Because measurements from CALIOP are available only during the volcanic active period and second quiescent period, the measurements from NOAA lidar 25 during the first quiescent period is excluded to simplify the comparison.
Generally speaking, a large color ratio (close to unity) indicates the presence of large particles and is directly associated with a small Ångström exponent (Jäger and Deshler, 2002). This was well documented during the Pinatubo eruption, when Ångström exponent values derived from balloon-borne measurements were very close to zero during the first months after the eruption and slowly increased as the plume dissipated (Jäger and Deshler, 2002). In this case, because the eruptions occurred In the case of high aerosol events, the color ratio of specific plumes can be derived based on the simultaneous measurements of NOAA lidar and MLSOL. Figure 9 presents two examples, one corresponding to the plume of the Nabro eruption and a the British Columbia (pyroCb) smoke plume (right). In the color ratio plots, the corresponding average color ratio derived for the volcanic active and second quiescent period (see Fig. 8, right plot) is included for reference (black, dashed). The altitude interval affected by these two features is highlighted in grey.
second example corresponding to the British Columbia smoke plume. In both cases, the plumes show higher color ratios than the surrounding. In the case of the Nabro example, the average color ratio above the plume is approximately 0.4, while in the plume the average is about 0.5. The values observed above the plume are also below the average for other cases (not shown).
For the PyroCb plume, the color ratio has a peak value of about 0.8 (corresponding to an Ångström exponent of -0.6), and values above the plume are in good agreement with the mean corresponding to the second quiescent period. 5 5.2 Lidar ratio of stratospheric volcanic aerosols measured after Nabro eruption (2011) Due to SNR limitations, direct measurements of stratospheric aerosol extinction profiles by Raman lidars are not common and typically restricted to strong volcanic events like the 1991 Mt. Pinatubo (Ferrare et al., 1992;Ansmann et al., 1993), Sarychev Peak (Mattis et al., 2010) and Nabro eruptions (Sawamura et al., 2012). For this reason, most of the long-term stratospheric aerosol lidar studies rely on a lidar ratio derived out of balloon-borne in-situ measurements (Jäger et al., 1995). In this study, 10 and thanks to the the large receiver, laser power and elevation of MLSOL, direct aerosol extinction measurements during the Nabro eruption are presented. These measurements allow, in combination with the derived backscatter profiles, the retrieval of the lidar ratio associated with the volcanic plume.
The results presented in Fig. 10, correspond to measurements conducted on 19 July 2011 by MLSOL. A well defined plume with a peak backscatter coefficient of 0.75 × 10 −3 sr −1 km −1 was found at 18.7 km. The corresponding extinction was 15 measured to be 4.8 × 10 −2 km −1 , leading to a lidar ratio of (64.0 ± 12.7) sr at 355 nm around the center of the plume. The The retrieval of the extinction coefficient out of the CALIOP measurements requires the assumption of a lidar ratio. In the case of the fraction of the GloSSAC dataset where CALIOP is included, this is assumed to be 53 sr, while in the CALIOP L3 data product and the profiles derived from MLSOL and NOAA lidar, a lidar ratio of 50 sr is used.
Another point to take into account while comparing the datasets is the differences in their spatial resolution. The GloSSAC dataset provides a zonal average of the extinction profile for a latitude bin of 5 • . In contrast, the CALIOP L3 dataset has 10 a longitudinal resolution of 20 • . In the case of MLSOL and NOAA lidar, and given their ground-based nature and uneven temporal sampling, extinction profiles correspond to a very small atmospheric volume making them more susceptible to small scale variability. With regard to the vertical resolution, and in order to calculate the corresponding AOD, an interpolation of the GloSSAC and CALIOP L3 datasets needs to be introduced in order to match the MLSOL and NOAA lidar grids.
Finally, it is important to notice that the backscatter retrievals from MLSOL are performed at 355 nm. In order to compare 15 them with the GloSSAC, CALIOP and NOAA lidar datasets, a conversion following the results from Sec. 5.1 is applied. For the conversion of the GloSSAC extinction dataset from 525 nm to 532 nm, an constant Ångström exponent of -1.6 is used (Jäger and Deshler, 2002).
The extinction coefficient time series for the five datasets under consideration are presented in Fig. 11. As can be seen, there is a general qualitative agreement between all the datasets. In the case of GloSSAC and CALIOP, and due to the higher amount 20 and more evenly distributed measurements included in the time series calculation, a smoother time-series can be observed. In At the bottom of the stratospheric aerosol layer, between 17 km ASL and 23 km ASL, aerosol plumes corresponding to different stratospheric injection events are clearly visible in all datasets. The largest contributions can be easily correlated to volcanic eruptions that occurred in the northern hemisphere (Table 1)  top altitude that can be associated with the quasi-biannual oscillation (Hommel et al., 2015). In this case, the maximum of the SAL was observed to be at about 32 km ASL, while its minimum can reach as low as 26 km ASL.
In order to provide a better overview of the differences between the datasets under study, the mean relative differences between MLSOL and the other four datasets as a function of altitude and time are presented in Fig. 12. The difference between the two ground-based lidars (Fig. 12d) shows a slight temporal-dependent extinction coefficient difference with higher values period and most of the second quiescent period. After the second half of 2018, slightly higher values are reported by NOAA lidar. For the difference between MLSOL and GloSSAC (Fig. 12a), a similar temporal-dependent variation is observed, with GloSSAC showing higher extinction values during the first quiescent period and lower values than MLSOL during the volcanic active and second quiescent periods. Below 20 km ASL, the effect of the zonal average on the GloSSAC dataset can be appreciated during volcanic injection events. The zonal average introduces a small shift in the plume temporal shape which 5 in turn translates into higher values reported by GloSSAC followed by higher values reported by MLSOL. While comparing MLSOL with CALIOP L3 dataset (Fig. 12b), higher extinction values are generally shown by CALIOP L3. This is opposite to what GloSSAC shows for the same time period (for which CALIOP and OSIRIS measurements are used). In order to discard the influence of the zonal averaging in GloSSAC, the CALIOP L3 datasets were zonally averaged. The results (not shown), did not reveal a significant influence of the zonal average that can explain this difference. Finally, the difference with SAGE-III ISS Figure 13 provides a more quantitative analysis of the differences between datasets by averaging the results provided in Figs. 11 and 12 over the three periods under study (first quiescent, volcanic active and second quiescent periods). Although the general agreement between datasets is good, with a maximum relative difference usually around 25 % above 23 km, 15 the differences tend to grow towards the bottom of the vertical profiles. Although in the case of MLSOL and NOAA lidar this difference likely originates from different retrieval approaches, including noise subtraction and normalization, systematic errors caused by misalignment cannot be ruled out. The relative impact of these error sources is currently under investigation.
Other differences, as the one shown with the new CALIOP L3 data product, seem to be more systematic, with CALIOP exhibiting consistently higher extinction values when compared with the other two satellite-based datasets and MLSOL. When 20 comparing with NOAA lidar, the CALIOP dataset also tends to show higher extinction values below 23 km, with the exception of the second half of the second quiescent period, where the agreement is good. The agreement of MLSOL with GloSSAC is good, with a relative difference consistently below 20 %.
As a final metric to evaluate the agreement between datasets, the corresponding stratospheric AOD at 532 nm is calculated for the period comprehended between January 1999 and November 2019 and presented in Fig. 14 Figure 13. Temporally-averaged extinction coefficient profiles at 532 nm for the three pre-defined periods (upper panels) presented together with the mean relative difference between MLSOL and each other available dataset for each period (lower panels). The mean relative difference is calculated as in Fig. 12. In the case of the second quiescent period, the average is split in two sections to take into account the partial temporal coverage of GloSSAC (only first half) and SAGE III (only second half).
of the agreement between datasets is presented in Table 2, where the mean relative differences between AOD time-series (X row,column ) are presented.

Summary and conclusions
Leveraging on the experience of previous studies on lidar stratospheric aerosol retrievals (e.g. Vernier et al., 2009)    backscatter and extinction retrieval process over the whole processing period. As an additional outcome for this study, a new version of the JPL MLSOL aerosol product is expected to be uploaded to the NDACC database in the near future.
A first analysis of the extinction time-series provided the opportunity to investigate the typical "tape-recorder-like" signature on the different volcanic plumes as well as the modulation of the SAL top by the QBO. By analyzing the evolution of the Sarychev plume centroid as function of time, an estimate of the ascent rate of the BDC was obtained (0.025 cm s −1 ) as it was 5 done in the past with satellite-based water vapor and aerosol observations. The result showed to be within 20% when compared with Minschwaner et al. (2016) and double the one reported in Vernier et al. (2011).
Based on this reprocessed long-term dataset and collocated NOAA lidar measurements, an independent retrieval of the stratospheric color ratio (and Ångström exponent) was obtained as a function of altitude and time. The results show a good agreement with previous observations based on balloon-borne measurements (Deshler et al., 2003) which increases the confi-10 dence on the values used for several previous studies. In a similar way, and thanks to the high SNR Raman channels available at MLSOL, the lidar ratio at the center of the aged Nabro plume was measured to be (64.0 ± 12.7) sr at 355 nm. Although critical for the derivation of extinction profiles and AOD time-series out of typical elastic backscatter lidars, lidar ratio measurements are sparse and typically show high variability and uncertainty. The observed lidar ratio for the presented case was shown to be in relative good agreement (overlapping 1σ uncertainties and a 16 % higher mean value for the MLSOL-derived observation) with previously reported measurements from the same event (Sawamura et al., 2012).
Recently released spaceborne long-term datasets of stratospheric aerosol extinction profiles, namely, CALIOP L3, GloS-SAC and SAGE III-ISS, were evaluated by comparing them with MLSOL and NOAA lidar profiles. While a generally good agreement between all datasets (GloSSAC, CALIOP L3 and SAGE III-ISS) and ground-based lidars was observed, the new 5 CALIOP L3 stratospheric data product showed consistently higher values in the lower stratosphere compared to MLSOL and SAGE III-ISS. Further intercomparison between this newly-released CALIOP data product and other ground-based datasets is required in order to assess the cause and the extent of this difference.
Finally, the stratospheric AOD at 532 nm between 17 km and 33 km ASL was derived from all datasets under intercomparison. The AOD time series show an increase of the AOD over time while comparing the two defined quiescent periods.

10
This observed increase in the background aerosol levels is in agreement with previously reported observations (Khaykin et al., 2017