Recovery of the first ever multi-year lidar dataset of the stratospheric aerosol layer, from Lexington, MA, and Fairbanks, AK, January 1964 to July 1965

We report the recovery and processing methodology of the first ever multi-year lidar dataset of the stratospheric aerosol layer. A Q-switched ruby lidar measured 66 vertical profiles of 694 nm attenuated backscatter at Lexington, Massachusetts, between January 1964 and August 1965, with an additional nine profile measurements conducted from College, Alaska, during July and August 1964. We describe the processing of the recovered lidar backscattering ratio profiles to produce mid-visible (532 nm) stratospheric aerosol extinction profiles (sAEP₅₃₂) and stratospheric aerosol optical depth (sAOD₅₃₂) measurements, utilizing a number of contemporary measurements of several different atmospheric variables. Stratospheric soundings of temperature and pressure generate an accurate local molecular backscattering profile, with nearby ozone soundings determining the ozone absorption, which are used to correct for two-way ozone transmittance. Two-way aerosol transmittance corrections are also applied based on nearby observations of total aerosol optical depth (across the troposphere and stratosphere) from sun photometer measurements. We show that accounting for these two-way transmittance effects substantially increases the magnitude of the 1964/1965 stratospheric aerosol layer's optical thickness in the Northern Hemisphere mid-latitudes, then ∼ 50 % larger than represented in the Coupled Model Intercomparison Project 6 (CMIP6) volcanic forcing dataset. Compared to the uncorrected dataset, the combined transmittance correction increases the sAOD₅₃₂ by up to 66 % for Lexington and up to 27 % for Fairbanks, as well as individual sAEP₅₃₂ adjustments of similar magnitude. Comparisons with the few contemporary measurements available show better agreement with the corrected two-way transmittance values.

Within the January 1964 to August 1965 measurement time span, the corrected Lexington sAOD₅₃₂ time series is substantially above 0.05 in three distinct periods, October 1964, March 1965, and May–June 1965, whereas the 6 nights the lidar measured in December 1964 and January 1965 had sAOD values of at most ∼ 0.03. The comparison with interactive stratospheric aerosol model simulations of the Agung aerosol cloud shows that, although substantial variation in mid-latitude sAOD₅₃₂ are expected from the seasonal cycle in the stratospheric circulation, the Agung cloud's dispersion from the tropics would have been at its strongest in winter and weakest in summer. The increasing trend in sAOD from January to July 1965, also considering the large variability, suggests that the observed variations are from a different source than Agung, possibly from one or both of the two eruptions that occurred in 1964/1965 with a Volcanic Explosivity Index (VEI) of 3: Trident, Alaska, and Vestmannaeyjar, Heimaey, south of Iceland. A detailed error analysis of the uncertainties in each of the variables involved in the processing chain was conducted. Relative errors for the uncorrected sAEP₅₃₂ were 54 % for Fairbanks and 44 % Lexington. For the corrected sAEP₅₃₂ the errors were 61 % and 64 %, respectively. The analysis of the uncertainties identified variables that with additional data recovery and reprocessing could reduce these relative error levels. Data described in this work are available at https://doi.org/10.1594/PANGAEA.922105 (Antuña-Marrero et al., 2020a).

1 Introduction

The abrupt enhancements to the stratospheric aerosol layer from historical large magnitude volcanic eruptions (e.g. Deshler, 2008) cause substantial radiative forcing of the Earth's climate system. Reducing their uncertainty remains an important priority since volcanic forcings can be the strongest driver of natural climate variability (e.g. Hansen, 1978; Robock, 2000). One of the coordinated multi-model experiments within the current international ISA-MIP activity (Interactive Stratospheric Aerosol Model Intercomparison Project; Timmreck et al., 2018) involves simulations of the volcanic aerosol clouds from the largest volcanic eruptions in the last century: Mt. Agung in 1963, El Chichón in 1982, and Mt. Pinatubo in 1991. One of the main motivations within this HErSEA multi-model experiment (Historical Eruption SO₂ Emissions Assessment) is to gather stratospheric aerosol observations in the periods after major tropical eruptions to provide new constraints to evaluate the model simulations. Another is to seek to understand whether the current diversity in the sulfur emission amount and altitude distribution that stratospheric aerosol models use when simulating the Pinatubo aerosol cloud is also seen for other major tropical eruptions such as Agung (see Sect. 3.3.2 of Timmreck et al., 2018). The first of the ISA-MIP modelling groups to present results from all three of the HErSEA eruption cloud experiments was recently published (Dhomse et al., 2020). Another recent study focused on assessing the variability in and global distribution of the Agung aerosol cloud (Niemeier et al., 2019).

Whereas the models participating in ISA-MIP simulate volcanic aerosol clouds interactively, the historical climate model simulations (Hegerl and Schwierz, 2011; Gillett et al., 2016) use prescribed volcanic forcing datasets (e.g. Sato et al., 1993; Ammann et al., 2003; Luo, 2016; Thomason et al., 2018). The observational data constraining the Agung aerosol cloud in both the interactive models and for the volcanic forcing datasets have tended to be based on column optical properties measured at the surface. These are primarily the extensive synthesis of surface radiation observations summarized by Dyer and Hicks (1968), with additional turbidity anomaly data from astronomical measurements of the atmospheric attenuation of starlight (Stothers, 2001). Although the literature includes several papers reporting profile measurements of the Agung aerosol cloud from balloon measurements (Rosen, 1964, 1968), lidars (Fiocco and Grams, 1964; Clemesha et al., 1966), and searchlights (Elterman and Campbell, 1964), no profile dataset of Agung backscatter ratio or aerosol extinction has yet been available to the scientific community. Whereas the Jamaica lidar (Clemesha et al., 1966) also measured the Agung cloud, the first multi-year lidar measurements of the Agung eruption were conducted from Lexington, Massachusetts, from January 1964 to August 1965 (Grams and Fiocco, 1967, hereinafter GF-67). No digital record of these lidar measurements existed until now, the data apparently only being presented in figures within published scientific papers. Only a few quantitative estimates of the cloud's optical properties from the lidar dataset have been found: aerosol extinction of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ km⁻¹ at 16 km and the aerosol optical depth of 0.015 (Deirmejian, 1971).

Table 1 Technical features of the lidars operated at Lexington and College, Fairbanks.

A problem with these early ruby lasers was the fluorescent emission which followed the laser pulse. These lidars incorporated a rotating shutter in the transmitting unit, synchronized with the Q-switching. The sensing unit for the backscattered signal consisted of an astronomical telescope with an interference filter and photomultiplier tube synchronized to another rotating shutter to avoid exposure to the intense returns from short distances. The photomultiplier was cooled with methanol and dry ice to reduce the dark current (Grams, 1966).

However, after initial failed searches of digital archives at several institutions, we discovered that the original lidar backscatter ratio profile measurements from the Lexington and Alaska 1964/1965 soundings are fully tabulated in Gerald W. Grams' PhD thesis conducted under the supervision of Giorgio Fiocco (Grams, 1966, hereinafter identified as G-66). Fortunately, at those times it was quite common for some observational datasets to be tabulated within PhD theses, grant reports, etc., a practice that after several decades is becoming required again, with many journals now mandating authors to make available the data they use via a recognized open-access data archive.

Dhomse et al. (2020) used preliminarily processed lidar data from Lexington, MA, one of the two sites reported in Grams and Fiocco (1967) to compare model aerosol extinction at 16 km with lidar observations, finding good agreement. Dhomse et al. (2020) and Niemeier et al. (2019) also noted the large change in the volcanic forcing for the Agung periods with the volcanic aerosol datasets used in the Coupled Model Intercomparison Projects 5 and 6 (CMIP5; Taylor et al., 2012; CMIP6; Eyring, et al., 2016; Zanchettin et al., 2016). The importance of reducing uncertainty by reconciling the datasets with additional stratospheric aerosol observations was also identified within these studies. Only an initial (preliminary) version of the Lexington 550 nm aerosol extinction dataset was used in Dhomse et al. (2020), with the analysis here producing a vertical profile dataset between 12 and 24 km. An important aspect of the dataset here is the two-way transmittance corrections applied to the aerosol backscatter ratio when deriving the aerosol extinction and optical depth datasets, with a detailed and transparent assessment of the relative error in each quantity also included.

This work is a contribution to the data rescue activity of the Stratospheric Sulfur and its Role in Climate (SSiRC), a SPARC (Stratosphere–Troposphere Processes and their role in Climate) initiative (SSiRC, 2020). The 1964/1965 lidar data recovered here follows on from another important volcanic aerosol dataset recovery of two ship-borne lidar datasets that measured the progression of the highly uncertain “tropical core” of the Pinatubo aerosol cloud from 17 July to 13 September 1991, 4–12 weeks after the 15 June 1991 Pinatubo eruption (Antuña-Marrero et al., 2020b). Those datasets were an identified priority within the SSiRC data rescue activity since they provide new constraints within the period when the Stratospheric Aerosols and Gas Experiment II (SAGE II) satellite could only observe the upper part of the Pinatubo aerosol cloud due to the saturation of the aerosol extinction retrieval (e.g. Thomason, 1992; McCormick and Veiga, 1992).

2 Materials and methods

2.1 Lidar instrumentation

The first successful laser radar ranging experiment was conducted at the Research Laboratory of Electronics, Massachusetts Institute of Technology, at Lexington, Massachusetts, and consisted of analysing the return signal from a short pulse (micro-second) laser covering the 60–140 km altitude range (Smullin and Fiocco, 1962). The research team, led by Giorgio Fiocco, continued developing applications of the lidar for atmospheric research. Scattering layers were detected in the upper atmosphere between 110 and 140 km (Fiocco and Smullin, 1963) and were interpreted to originate from meteoric fragments entering the outer atmosphere (Fiocco and Colombo, 1964). After some changes and improvements, stratospheric aerosols were detected between 10 and 30 km altitude, and the first lidar measurements of the stratospheric aerosol layer began (Fiocco and Grams, 1964).

The schematic diagram and a photo of the instrument are in Figs. 3 and 4 of G-66, respectively. Also listed are the main features of the lidars used for the measurements at Lexington and College, Alaska, reproduced in Table 1. Both lidars used a Q-switched ruby laser at the 694 nm wavelength.

2.2 Lidar measurements

Lidar observations were conducted at Lexington, Massachusetts (42^∘25^′ N, 71^∘15^′ W), and at College (64^∘53^′ N, 148^∘3^′ W), located in the city of Fairbanks, Alaska, hereinafter identified as Fairbanks. The measurements were supported by the NASA Grant NGR-22-009-131. One of the semi-annual reports mentions more than 100 measurements conducted (Fiocco, 1966a). However the number of profiles appearing in Grams' PhD dissertation was 75. A total of 9 d of measurements from 26 July to 21 August 1964 were conducted at Fairbanks. At Lexington, Massachusetts, 23 d of measurements from 14 January to 20 May 1964 and 43 d from 11 October 1964 to 21 July 1965 were made. At both sites, measurements were restricted to dark nighttime conditions. A single laser shot was registered by photographing the return signal on an oscilloscope covering up to 40 km and then digitized by hand. The digitized return signals from a set of laser shots were then averaged in 1 km bins (G-66; GF-67).

2.3 Backscattering ratios in the original lidar dataset

It is well known that solving the lidar equation for the single-wavelength elastic lidar is an ill-posed problem. The return signal is the result of scatter from both molecules and aerosol particles; hence additional information is necessary to separate their contributions (e.g. Kovalev, 2015). Considering this fact, we may understand the magnitude of the challenge confronted by Giorgio Fiocco and then Gerald W. Grams, Giorgio Fiocco's PhD student, when they conducted the processing of the first ever set of lidar return signals from stratospheric aerosols.

First, we describe the procedure applied in G-66 to derive the backscattering ratio at 694 nm and range z, SR(694, z). The average photoelectron flux (electrons per second) registered by the photomultiplier, which is proportional to the backscattered signal, was represented in Eq. (3.8) of G66 as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}n\left(z\right)}{\mathrm{d}t}}=K{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{T}}\left(z\right)}{z^{\mathrm{2}}}},\end{array}$$

where z is the altitude, σ_T(z) the total radar cross-section per unit volume of atmospheric constituents at z, and K is a constant resulting from all the terms not depending on the altitude in the optical radar equation, including T_2w , the two-way atmospheric transmittance (see G-66 for more details). The assumption of a constant value for T_2w in the stratosphere was based on the atmospheric attenuation model proposed by Elterman (1964). The model provided magnitudes of the molecular and aerosol scattering, as well as the ozone absorption, showing that almost all attenuation of the laser beam occurs in the troposphere. The model gave an estimate of the variability in the term T_2w, at 700 nm in the stratosphere, between 10 and 30 km, which was below 3 %. The correction of the return signal, associated with the two-way transmittance of the laser beam throughout the atmosphere, was then neglected, and it was assumed that the atmospheric extinction term was constant. This is a good assumption for times of low stratospheric aerosol loading. For enhanced stratospheric aerosol, e.g. after volcanic eruptions, however, aerosol extinction becomes important, reduces the stratospheric transmission, and makes it range dependent.

The return signal from a set of laser shots was averaged in time and in altitude to a resolution of 1 km between 12 and 30 km. Next, the ratios between the averaged signal at each level and the values at the same level of the right side of Eq. (1) were calculated for each profile between 12 and 30 km. A final step for each profile consisted in normalizing the ratios calculated in each profile between 12 and 24 km, with the average ratios between 25 and 30 km producing the derived SR(694, z) under the assumptions already cited. The normalization procedure assumed the contribution from aerosols was negligible above 24 km, leading to an underestimate of stratospheric aerosol since there would have been aerosol at these altitudes (Russell et al., 1979).

The SR(694, z) values derived from the lidar measurements conducted at Lexington and Fairbanks were reported in tabular format in Gerald W. Grams' PhD thesis (G-66) and cited in the acknowledgements section of GF-67. It was the unique reference of its existence that was the clue that guided us in our search for the lidar measurements.

2.4 Algorithms used in the processing

The SSiRC data rescue activity is committed, whenever it is possible, to re-calibrate each dataset and determine its levels of uncertainties (SSiRC, 2020). Because some stratospheric aerosol lidar datasets have already been identified and located, we endeavour to reprocess them using a standardized algorithm to guarantee the best possible consistency among the different lidar datasets. Below, we describe the processing algorithm as a first step in that direction.

The lidar backscattering ratio, SR(λ, z), is commonly defined as the ratio between the total backscatter β_T(λ,z) and the molecular backscatter β_m(λ,z) at the altitude z and wavelength λ. β_T(λ,z) is the sum of β_m(λ,z) and the aerosol backscatter, β_a(λ,z):

$$\begin{array}{}\text{(2)}& \mathrm{SR}\left(\mathit{\lambda },z\right)={\displaystyle \frac{{\mathit{\beta }}_{\mathrm{m}}\left(\mathit{\lambda },z\right)+{\mathit{\beta }}_{\mathrm{a}}(\mathit{\lambda },z)}{{\mathit{\beta }}_{\mathrm{m}}\left(\mathit{\lambda },z\right)}}.\end{array}$$

β_m(λ,z) is derived using the following equation:

$$\begin{array}{}\text{(3)}& {\mathit{\beta }}_{\mathrm{m}}(\mathit{\lambda },z)={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{m}}(\mathit{\lambda },z)}{S_{\mathrm{m}}}}={\displaystyle \frac{\mathrm{3}{\mathit{\sigma }}_{\mathrm{m}}(\mathit{\lambda },z)}{\mathrm{8}\mathit{\pi }}},\end{array}$$

where $S_{\text{m}}=(\mathrm{8}\mathit{\pi }/\mathrm{3})k_{\text{bw}}$ is the molecular extinction to backscatter ratio for the molecular scattering, commonly approximated by $\mathrm{8}\mathit{\pi }/\mathrm{3}$ (Collins and Russell, 1976) after neglecting the dispersion of the refractive index and the King factor of the air represented by k_bw. The volume molecular scattering coefficient, σ_m(λ,z), is determined by the following equation:

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{\mathrm{m}}(\mathit{\lambda },z)={\displaystyle \frac{N_{\mathrm{A}}\mathrm{Pr}\left(z\right)}{R_{\mathrm{a}}\mathrm{Temp}\left(z\right)}}{Q}_{\mathrm{s}}\left(\mathit{\lambda }\right),\end{array}$$

where $N_{\text{A}}=\mathrm{6.02214}\times {\mathrm{10}}^{\mathrm{23}}$ (mol⁻¹) is Avogadro's number, R_a= 8.314472 (J K⁻¹ mol⁻¹) is the gas constant, and Q_s(λ) the total molecular scattering cross-section per molecule for the standard air. The derived equation for Q_s(λ) for standard air is (Hostetler et al., 2006)

$$\begin{array}{}\text{(5)}& Q_{\mathrm{s}}\left(\mathit{\lambda }\right)=\mathrm{4.5102}\times {\mathrm{10}}^{-\mathrm{27}}{\left[{\displaystyle \frac{\mathit{\lambda }\left(\mathrm{nm}\right)}{\mathrm{550}}}\right]}^{-\mathrm{4.025}-\mathrm{0.05627}\times {\left[\frac{\mathit{\lambda }\left(\mathrm{nm}\right)}{\mathrm{550}}\right]}^{-\mathrm{1.647}}}.\end{array}$$

The single-wavelength elastic lidar systems provide profiles of attenuated total backscatter. The “true” total backscatter is calculated from ${\mathit{\beta }}_{\mathrm{T}}(\mathit{\lambda },z)={\mathit{\beta }}_{\mathrm{T}}^{\mathrm{A}}\left(\mathit{\lambda },z\right)T_j(\mathit{\lambda },z)$, where ${\mathit{\beta }}_{\mathrm{T}}^{\mathrm{A}}(\mathit{\lambda },z)=\left[{\mathit{\beta }}_{\mathrm{m}}^{\mathrm{A}}\right(\mathit{\lambda },z)+{\mathit{\beta }}_{\mathrm{a}}^{\mathrm{A}}(\mathit{\lambda },z\left)\right]$. Substituting into Eq. (2), taking into account that ${\mathit{\beta }}_{\mathrm{m}}(\mathit{\lambda },z)={\mathit{\beta }}_{\mathrm{m}}^{\mathrm{A}}(\mathit{\lambda },z)T_j(\mathit{\lambda },z)$, and re-arranging specifically for the case of the 694 nm backscattering ratio give

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathit{\beta }}_{\mathrm{a}}^{\mathrm{A}}(\mathit{\lambda },z)& ={\displaystyle \frac{\left[\mathrm{SR}(\mathrm{694},z)-\mathrm{1}\right]{\mathit{\beta }}_{\mathrm{m}}(\mathit{\lambda },z)}{T_j(\mathit{\lambda },z)}}\\ & ={\displaystyle \frac{\left[\mathrm{SR}(\mathrm{694},z)-\mathrm{1}\right]{\mathit{\beta }}_{\mathrm{m}}(\mathit{\lambda },z)}{T_{\mathrm{m}}(\mathit{\lambda },z)T_{{\mathrm{O}}_{\mathrm{3}}}(\mathit{\lambda },z)}}.\end{array}\end{array}$$

Whereas the term T_j(λ,z) in Eq. (6) usually specifies the full two-way transmittance correction, at this stage of the processing methodology, we correct only for the attenuation due to molecular backscatter and ozone absorption:

$$\begin{array}{}\text{(7)}& T_j(\mathit{\lambda },z)=e^{-\mathrm{2}{\int }_{z_o}^z{\mathit{\alpha }}_j(\mathit{\lambda },z)\text{d}z}.\end{array}$$

The α_j(λ,z) term is thus the vertical profile of 694 nm extinction due only to molecules, α_m(λ,z), and ozone absorption, ${\mathit{\alpha }}_{{\mathrm{O}}_{\mathrm{3}}}(\mathit{\lambda },z)$, with z_o being the altitude of the lidar. This postponement of the aerosol attenuation correction is due to the 500 nm wavelength of the contemporaneous aerosol optical depth (AOD) measurements being much closer to the target wavelength of 532 nm, the method preserving the two-way molecular and ozone transmittance corrections to be applied at the measurement wavelength of 694 nm.

The conversion to aerosol extinction is carried out after first translating the aerosol backscatter from 694 to 532 nm and then applying the corresponding wavelength exponent kb(z,t) calculated from in situ size distribution measurements of the Northern Hemisphere mid-latitude Pinatubo aerosol cloud (Jäger and Deshler, 2002, 2003):

$$\begin{array}{}\text{(8)}& {\mathit{\beta }}_{\mathrm{a}}(\mathrm{532},z)={\left[{\displaystyle \frac{\mathrm{532}}{\mathrm{694}}}\right]}^{\mathrm{kb}(z,t)}{\mathit{\beta }}_{\mathrm{a}}(\mathrm{694},z).\end{array}$$

Note again that the derived 532 nm aerosol backscatter, β_a(532,z), has at this point only been corrected for two-way molecular scattering and ozone absorption transmittance effects.

The aerosol extinction, α_a(532,z), at this point still uncorrected for two-way aerosol transmittance, is then calculated by the following expression:

$$\begin{array}{}\text{(9)}& {\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z)={\mathrm{EB}}_{\mathrm{c}}(z,t){\mathit{\beta }}_{\mathrm{a}}(\mathrm{532},z),\end{array}$$

where EB_c(z,t) are altitude- and time-dependent coefficients to convert aerosol backscatter to aerosol extinction (at λ=532 nm), derived from the same Pinatubo aerosol size distribution measurements (Jäger and Deshler, 2003).

Both the EB_c and kb factors are derived from log-normal size distribution fits to balloon-borne optical particle counter measurements of the Pinatubo aerosol cloud from Laramie, Wyoming, USA. Each of the conversion factors in Jäger and Deshler (2002, 2003) represent averages over four height ranges – tropopause–15, 15–20, 20–25, and 25–30 km – and are provided for the 4-month periods November–February, March–June, and July–October of each year after the eruption. We used EB_c(z,t) and kb(z,t) for the same 4-month periods after the March 1963 Agung eruption, based on matching the same time offset after the Mt. Pinatubo 1991 eruption.

The algorithms for the solution of the single wavelength lidar equations apply the two-way transmittance correction to the raw lidar return signal, together with squared distance correction, before the backscattering ratio is calculated. In our case the only measurement information we have is the lidar backscatter ratios, which have been derived without conducting the two-way transmittance correction (G-66) for any species. However, only the molecular and ozone two-way transmittance corrections, T_m(694,z) and $T_{{\mathrm{O}}_{\mathrm{3}}}(\mathrm{694},z)$, were included in Eq (6) by the reasons explained below.

As explained above, the aerosol two-way transmittance correction, T_a(532,z) , was deliberately postponed until the final step to derive ${\mathit{\alpha }}_{\text{a}}^{T_{\text{a}}}(\mathrm{532},z)$ due to the available contemporaneous measurement information for AOD being at λ=500 nm (see Supplement). We converted the measured AOD at 500 to 532 nm using Ångstrom exponents covering the cited wavelength range from 1995 to 2019 from the nearest Aerosol Robotic Network (AERONET, 2020) stations. Although the tropospheric aerosol layer in the eastern USA will have had different physical and chemical properties in the 1960s (e.g. Went, 1960; Husar et al., 1991), this is only a small change in wavelength, with the method then introducing much less error than had we converted from 500 to 694 nm to apply the aerosol attenuation correction. We note that the calculated monthly mean total AOD (TAOD) at 532 nm from Blue Hill Observatory, MA, from 1961 to 1966 is in the range from 0.1 to 0.4, consistent with the elevated background TAOD reported for the eastern USA during the 1960s (Husar et al., 1981; Sect. S1 in the Supplement).

We produced a first guess $T_{\text{a}}(\mathrm{532},z{)}_{*}$ for each measurement day at each site in the range from 11 to 24 km. The $T_{\text{a}}(\mathrm{532},\mathrm{11}\text{km}{)}_{*}$, a unique value from the lidar altitude to 11 km, was calculated using Eq. (7) and the TAOD value for the month in which the measurement was conducted. $T_{\text{a}}(\mathrm{532},z{)}_{*}$ values from 12 to 24 km were calculated using Eq. (7) and the uncorrected α_a(532,z).

The first guess aerosol scattering corrected by the total two-way transmittance, ${\mathit{\alpha }}_{\text{a}}^{T_{\text{a}}}{\left(\mathrm{532},z\right)}_{*}$, was derived by applying the correction for the two-way aerosol transmittance T_a(532,z):

$$\begin{array}{}\text{(10)}& {\mathit{\alpha }}_{\text{a}}^{T_{\text{a}}}(\mathrm{532},z)={\displaystyle \frac{{\mathit{\alpha }}_{\text{a}}(\mathrm{532},z)}{T_{\text{a}}(\mathrm{532},z)}}.\end{array}$$

Since we are using the measured TAOD, which included the stratospheric AOD, to calculate the first guess $T_{\text{a}}(\mathrm{532},z{)}_{*}$, we applied a second step after producing a first guess ${\mathit{\alpha }}_{\text{a}}^{T_{\text{a}}}(\mathrm{532},z{)}_{*}$ profile, for which we then calculate the stratospheric AOD (sAOD_*), integrating ${\mathit{\alpha }}_{\text{a}}^{T_{\text{a}}}(\mathrm{532},z{)}_{*}$ between 12 and 24 km. Then sAOD_* is used in Eq. (11) to estimate the tropospheric AOD (tAOD) for each measurement:

$$\begin{array}{}\text{(11)}& \text{tAOD}=\text{TAOD}-{\text{sAOD}}_{*}.\end{array}$$

The former $T_{\text{a}}(\mathrm{532},\mathrm{11}\text{km}{)}_{*}$ values for each measurement will be replaced by new ones derived using the calculated tAOD corresponding to each measurement. Then the final profile of T_a(532,z) for each measurement will consist of the new T_a(532,11 km) calculated using tAOD and the already derived T_a(532,z) values from 12 to 24 km that were calculated using the uncorrected α_a(532,z) in Eq. (7). Those profiles of T_a(532,z) are applied in Eq. (10) producing the definitive values of ${\mathit{\alpha }}_{\text{a}}^{T_{\text{a}}}(\mathrm{532},z)$.

The method we applied uses monthly means of the TAOD from the sun photometer, converted from 500 to 532 nm (described above), assuming this to be tropospheric (tAOD) in the first step, to produce for each site a $T_{\mathrm{a}}(\mathrm{532},z{)}_{*}$ profile so as to account for the tropospheric aerosol transmittance from the lidar altitude across the troposphere up to 11 km and the stratospheric aerosol transmittance in the lower stratosphere from 12 to 24 km. However, when the stratospheric AOD (sAOD_*) is calculated in the next step (from the first guess ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z{)}_{*}$), the resulting first guess total AOD (tAOD + sAOD_*) will be higher than the observed TAOD. The second step is aimed at estimating the magnitude of a consistent value of tAOD for each measurement to constrain tAOD + sAOD ≤ TAOD.

Although more iterations of those final steps would be possible, with the high magnitude of the estimated ${\mathit{\alpha }}_{\text{a}}^{T_{\text{a}}}(\mathrm{532},z)$ mean error, around 60 %, compared to an estimated 15 %–20 % maximum improvement achieved by the iteration procedure, we do not believe those additional calculations would be worthwhile in this case.

2.5 Complementary datasets used

The correction for the attenuation of the lidar signal by the two-way transmission by atmospheric molecules, ozone, and aerosols is often considered negligible and ignored based on signal to noise ratio considerations or for simplicity (e.g. G-66; GF-67). We were motivated to make that correction by the fact that the accuracies of the different instruments available for measurements of the stratospheric aerosols from the 1963 Mt. Agung eruption are still under scrutiny and discussion (e.g. Deirmendjian, 1965, 1971; Dyer, 1971a, b; Clemesha, 1971; Stothers, 2001; Timreck et al., 2018). Our goal was to produce the consistently processed aerosol extinction and optical depth from the rescued measurements based on the contemporary state of the art measurements in the 1960s. For the different data sources and processing algorithms, we calculate the two-way transmittance corrections by atmospheric molecules, ozone, and aerosols as described in Sect. S1.

2.6 Numerical and statistical methods

For each of the two datasets we calculate percentage differences (Δα_*US) between α_a(532,z)_US calculated using the same β_m(694,z) profile from the 1962 US Standard Atmosphere for all the days and the ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$ calculated using the β_m(694,z) profiles derived from the daily soundings,

$$\begin{array}{}\text{(12)}& \mathrm{\Delta }{\mathit{\alpha }}_{{\mathrm{a}}_{*}}={\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{\mathrm{US}}-{\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*},\end{array}$$

and the percent differences $\mathrm{\Delta }\mathit{\alpha }{\mathit{\%}}_{{\mathrm{a}}_{*}}$ by the following expression:

$$\begin{array}{}\text{(13)}& \mathrm{\Delta }{\mathit{\alpha }}_{{\mathrm{a}}_{*}}\mathit{\%}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{\mathrm{US}}-{\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}}{{\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{\mathrm{US}}}}\times \mathrm{100}.\end{array}$$

Similarly we defined the differences, Δα_at2w, and the percent differences, Δα_at2w %, between the ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$ calculated using the β_m(694,z) profiles derived from the daily soundings and its corrected values, α_a(532,z)_t2w, resulting from accounting for the two-way atmospheric transmittance.

Also, for cumulative aerosol optical depth in the 12 to 24 km layer, we define ${\mathit{\tau }}_{\mathrm{a}}(\mathrm{532},{)}_{*}$ and τ_a(532,z)_US calculated from the ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$ and α_a(532,z)_US, respectively, as

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }{\mathit{\tau }}_{{\mathrm{a}}_{*}}={\mathit{\tau }}_{\mathrm{a}}(\mathrm{532},z{)}_{\mathrm{US}}-{\mathit{\tau }}_{\mathrm{a}}(\mathrm{532},z{)}_{*},\end{array}$$

and the percent differences, $\mathrm{\Delta }\mathit{\tau }{\mathit{\%}}_{{\mathrm{a}}_{*}}$, by the following expression:

$$\begin{array}{}\text{(15)}& \mathrm{\Delta }{\mathit{\tau }}_{{\mathrm{a}}_{*}}\mathit{\%}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{a}}(\mathrm{532},z{)}_{\mathrm{US}}-{\mathit{\tau }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}}{{\mathit{\tau }}_{\mathrm{a}}(\mathrm{532},z{)}_{\mathrm{US}}}}\times \mathrm{100}.\end{array}$$

2.7 Relative error estimates

The present evaluation of the relative errors in the different processing steps of the single wavelength elastic lidar followed the algorithms developed by Russell (1979). Whenever it was possible, we calculated the different terms of the equation based on the available dataset error. In several cases we combined information from the rescued metadata associated with the measurements and from available additional information in the literature.

2.7.1 Backscattering ratio relative error

We use Eq. (19) from Russell (1979) quantifying the contributions from the different sources to the relative error in backscattering ratio $\frac{\mathit{\delta }\mathrm{SR}}{\mathrm{SR}}$:

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\left({\displaystyle \frac{\mathit{\delta }\mathrm{SR}}{\mathrm{SR}}}\right)}^{\mathrm{2}}& ={\left({\displaystyle \frac{\mathit{\delta }{N}_{\mathrm{s}}}{N_{\mathrm{s}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{T}}}{T_{\mathrm{T}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{m}}}{{\mathit{\beta }}_{\mathrm{m}}}}\right)}^{\mathrm{2}}+\left({\displaystyle \frac{\mathit{\delta }{\mathit{\beta }}_{{\mathrm{m}}^{*}}}{{\mathit{\beta }}_{{\mathrm{m}}^{*}}}}\right)\\ & -\left({\displaystyle \frac{{\mathrm{C}}_{{\mathrm{FF}}^{*}}^{\mathrm{2}}}{{\mathit{\beta }}_{\mathrm{m}}{\mathit{\beta }}_{{\mathrm{m}}^{*}}}}\right)+{\left({\displaystyle \frac{\mathit{\delta }\mathrm{S}{R}_{\min }}{{\mathrm{SR}}_{\min }}}\right)}^{\mathrm{2}},\end{array}\end{array}$$

where SR(λ,z) is the total backscattering ratio; N_s is the signal measured; T_T is the two-way transmittance from aerosols, molecules, and ozone; β_m is the molecular backscattering; ${\mathit{\beta }}_{{\mathrm{m}}^{*}}$ is molecular backscatter at the normalization level; SRmin(λ,z) is total backscattering ratio at the normalization level; and $C_{{\mathrm{FF}}^{*}}^{\mathrm{2}}$ is the covariance between measured β_m and ${\mathit{\beta }}_{{\mathrm{m}}^{*}}$.

For estimating the magnitude of the signal measurement error we rely on the information provided by G-66. He estimated statistical fluctuation of the signal, the shot noise of the photodetector, and other sources of the order of 0.2 % to 3 %. For both Lexington and Fairbanks we assume $\left(\frac{\mathit{\delta }{N}_{\mathrm{s}}}{N_{\mathrm{s}}}\right)=\mathrm{3}$ %.

As cited above, according to G-66 if no T_T correction was conducted, then the term ${\left(\frac{\mathit{\delta }{T}_{\mathrm{T}}}{T_{\mathrm{T}}}\right)}^{\mathrm{2}}=\mathrm{0}$.

Because in the calculation of SR(λ,z), values of N_d(z) from the 1962 US Standard Atmosphere were used (G-66), it was assumed $\frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{m}}(\mathit{\lambda },Z)}{{\mathit{\beta }}_{\mathrm{m}}(\mathit{\lambda },z)}=\mathrm{3}$ % for both sites (e.g. Russell et al., 1979). In addition we assumed $\left(\frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{m}}}{{\mathit{\beta }}_{\mathrm{m}}}\right)=\left(\frac{\mathit{\delta }{\mathit{\beta }}_{{\mathrm{m}}^{*}}}{{\mathit{\beta }}_{{\mathrm{m}}^{*}}}\right)$ and $C_{{\mathrm{FF}}^{*}}^{\mathrm{2}}=\mathrm{0}$ after assuming measurement errors are uncorrelated. The use at the lidar levels of interpolated β_m values from the lower-resolution ones calculated using the US 1962 Standard Atmosphere support the former assumption.

The term δSR_min was calculated according to Table (1b) in Russell (1979) for SR_min=1.01 and the respective latitudes of both sites. Then following Russell et al. (1979), we assume

$$\begin{array}{}\text{(17)}& \mathit{\delta }{\mathrm{SR}}_{\min }=\mathrm{0.07}\left({\mathrm{SR}}_{\max }-\mathrm{1}\right).\end{array}$$

2.7.2 Aerosol backscattering relative errors

Equation (18) in Russell (1979), estimating the relative error in β_a(694,z), can be approximated in our case by

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& {\left({\displaystyle \frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{a}}(\mathrm{694},z)}{{\mathit{\beta }}_{\mathrm{a}}(\mathrm{694},z)}}\right)}^{\mathrm{2}}=\\ & {\left({\displaystyle \frac{{\mathit{\beta }}_{\mathrm{m}}}{{\mathit{\beta }}_{\mathrm{a}}}}\right)}^{\mathrm{2}}\left\{{\left(\mathrm{SR}\right)}^{\mathrm{2}}\left[{\left({\displaystyle \frac{\mathit{\delta }\mathrm{SR}}{\mathrm{SR}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{T}}}{T_{\mathrm{T}}}}\right)}^{\mathrm{2}}\right]+{\left({\displaystyle \frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{m}}}{{\mathit{\beta }}_{\mathrm{m}}}}\right)}^{\mathrm{2}}\right\}.\end{array}\end{array}$$

The estimated error for the two-way transmission corrections in Russell et al. (1979) provides the following expression,

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& {\left({\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{T}}}{T_{\mathrm{T}}}}\right)}^{\mathrm{2}}=\\ & \mathrm{4}\mathit{\left\{}\right[\mathit{\delta }{\mathit{\tau }}_{\mathrm{a}}(\mathit{\lambda },z){]}^{\mathrm{2}}+\left[\mathit{\delta }{\mathit{\tau }}_{\mathrm{m}}\right(\mathit{\lambda },z){]}^{\mathrm{2}}+[\mathit{\delta }{\mathit{\tau }}_{{\mathrm{O}}_{\mathrm{3}}}(\mathit{\lambda },z){]}^{\mathrm{2}}\mathit{\}},\end{array}\end{array}$$

and considering the standard error of determination of τ_a, ${\mathit{\tau }}_{{\mathrm{O}}_{\mathrm{3}}}$, and τ_m being respectively 50 %, 20 %, and 10 %, the following estimates are produced: δτ_a=0.5τ_a, $\mathit{\delta }{\mathit{\tau }}_{{\mathrm{O}}_{\mathrm{3}}}=\mathrm{0.2}{\mathit{\tau }}_{{\mathrm{O}}_{\mathrm{3}}}$, and δτ_m=0.1τ_m . However, in our calculation of β_a, only the ozone and molecular two-way transmittances were used.

For this section of the procedure, we consider $\left(\frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{m}}}{{\mathit{\beta }}_{\mathrm{m}}}\right)=\mathrm{10}$ %. We neglected the error in computing Q_s using Eq. (5) because its maximum relative error is 0.2 % for a spectral region of 350–1600 nm (Hostetler et al., 2006), well below the errors in $\left(\frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{m}}}{{\mathit{\beta }}_{\mathrm{m}}}\right)$.

Next, we determined the relative error in β_a(532,z) associated with the conversion from β_a(694,z) in Eq. (7) using the wavelength exponents, kb(z,t), for aerosol backscatter in the range of wavelengths between 694 and 532 nm (Jäger and Deshler, 2003). The errors were estimated from their Fig. 1 with ${\left(\frac{\mathit{\delta }\text{kb}}{\text{kb}}\right)}^{\mathrm{2}}=\mathrm{10}$ %:

$$\begin{array}{}\text{(20)}& {\left({\displaystyle \frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{a}(\mathrm{532},z)}}{{\mathit{\beta }}_{\mathrm{a}}(\mathrm{532},z)}}\right)}^{\mathrm{2}}={\left({\displaystyle \frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{a}}(\mathrm{694},z)}{{\mathit{\beta }}_{\mathrm{a}}(\mathrm{694},z)}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathit{\delta }\text{kb}}{\text{kb}}}\right)}^{\mathrm{2}}.\end{array}$$

Figure 1 (a) Contour plot of the vertical profiles of 532 nm aerosol extinction α_a(532,z) calculated using the same β_m(694,z) profile from the 1962 US Standard Atmosphere for all the days. (b) The contour plot of the vertical profiles of 532 nm aerosol extinction α_a(532,z) was calculated using the daily β_m(694,z) profiles from the sounding at Nantucket, MA. The red stars indicate the dates the measurements were conducted. The measurement gaps longer than 1 month – March, and July to September both in 1964 – have been left blank.

Table 2 Relative differences between the ${\mathit{\alpha }}_{{\mathrm{a}}_{\mathrm{US}}}$ and ${\mathit{\alpha }}_{{\mathrm{a}}_{*}}$, as well as ${\mathit{\tau }}_{{\mathrm{a}}_{\mathrm{US}}}$ and ${\mathit{\tau }}_{{\mathrm{a}}_{*}}$.

2.7.3 Aerosol extinction relative errors

In the case of the α_a, its relative errors are

$$\begin{array}{}\text{(21)}& {\left({\displaystyle \frac{\mathit{\delta }{\mathit{\alpha }}_{\mathrm{a}}}{{\mathit{\alpha }}_{\mathrm{a}}}}\right)}^{\mathrm{2}}={\left({\displaystyle \frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{a}}}{{\mathit{\beta }}_{\mathrm{a}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathit{\delta }{\mathrm{EB}}_{\mathrm{c}}}{{\mathrm{EB}}_{\mathrm{c}}}}\right)}^{\mathrm{2}}.\end{array}$$

The last term in the right side represents the error in the EB_c for λ=532 nm. In the case of the ones we used (Jäger and Deshler, 2002, 2003) the error has been estimated at ±40 % according to Deshler et al. (2003). For ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}$, the aerosol extinction corrected by the aerosols' two-way transmittance using the estimates of its relative error described above is

$$\begin{array}{}\text{(22)}& {\left({\displaystyle \frac{\mathit{\delta }{\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}}{{\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}}}\right)}^{\mathrm{2}}={\left({\displaystyle \frac{\mathit{\delta }{\mathit{\alpha }}_{\mathrm{a}}}{{\mathit{\alpha }}_{\mathrm{a}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\mathit{\delta }{T}_{\mathrm{2}\mathrm{wa}}}{T_{\mathrm{2}\mathrm{wa}}}}\right)}^{\mathrm{2}}.\end{array}$$

Using the cited set of equations and the assumptions described above we evaluated the error for each altitude in each measurement.

3 Results

3.1 The 532 nm aerosol extinction profiles and optical depth

In Fig. 1 we show contour plot of the vertical profiles of 532 nm aerosol extinction, α_a(532,z), for Lexington calculated using the same β_m(694,z) profile from the 1962 US Standard Atmosphere for all the days and the daily β_m(694,z) profiles derived from the sounding at Nantucket, MA. If measurement gaps are longer than 1 month (March 1964, and July to September 1964), they have been left blank. The temporal and vertical contour plots of the aerosol extinctions were generated using a linear time interpolation.

In general, the contour plots show a high level of variability in the aerosol extinction for Lexington both in time and altitude associated with the complex thermodynamic processes in the upper troposphere–lower stratosphere. Three main maximums are identified across the entire period. The first is between 16 and 18 km at the beginning of the record in mid-January 1964, the second between 14 and 16 km by November 1964, and the third at the same altitude but in the transition between March and April 1965. Evident is the decaying altitude of the maximums in time typical of the volcanic aerosol clouds in the lower stratosphere. However, the occurrence of the absolute maximum at this time cannot be attributed to the volcanic aerosols from Mt. Agung, as will be discussed below. No long-term analysis of this type could be conducted in Fig. 2 for Fairbanks because of the very short period of time it covers. However, the cross-section of ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$ for Fairbanks reveals maximum values between 14 and 16 km, with the absolute maximum around mid-August centred at 15 km. The magnitudes of ${\mathit{\alpha }}_{\mathrm{a}(\mathrm{532},z{)}_{\mathrm{US}}}$ are slightly higher than the ones from ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$ for both sites, and it is also true for ${\mathit{\tau }}_{{\mathrm{a}}_{\mathrm{US}}}$ and ${\mathit{\tau }}_{{\mathrm{a}}_{*}}$. This is quantified in Table 2. The magnitude of the mean percent difference increase in both variables is around 1 %.

This difference disagrees with G-66 in that he found retrievals using the 1962 US Standard Atmosphere slightly lower than the more realistic ones using soundings, but the differences are within calculated errors. He arrived at that conclusion from “a cursory examination” of the local variations of molecular number density, N_d(z), estimated with the Temp(z) profiles from ozone soundings at Bedford, MA (Hering and Borden, 1967). He reported N_d(z) variability rarely exceeded 5 % of the N_dUS(z) values at altitudes between 10 and 30 km.

Figure 2 The same as Fig. 1 but for Fairbanks, AK.

Figure 3 Differences between the number molecular density, N_d(z), from soundings and from the 1962 US Standard Atmosphere in the region from 12 to 24 km. Panel (a) represents N_d(z) from Nantucket soundings used for Lexington and (b) N_d(z) from Fairbanks.

Figure 4 Stratospheric aerosol optical depth (sAOD₅₃₂) for Lexington (blue stars), Fairbanks (red diamonds), and for the Northern Hemisphere (black asterisks) for the period in which the measurements were conducted. The sAOD₅₃₂ values were calculated from the ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$ derived using local soundings. Blue stars and red diamonds on the top axis of the figure are the dates on which the measurements were conducted.

To estimate the effects of the differences between the magnitudes of N_dUS(z) and N_d(z) in the backscattering ratios, we calculate the differences between the ratios defined by

$$\begin{array}{}\text{(23)}& \mathrm{\Delta }{N}_{\mathrm{d}}\left(z\right)={\displaystyle \frac{N_{\mathrm{dUS}}\left(z\right)}{M_{\mathrm{dUS}}}}-{\displaystyle \frac{N_{\mathrm{d}}\left(z\right)}{M_{\mathrm{d}}}}.\end{array}$$

M_dUS and M_d are the mean values of N_dUS(z) and N_d(z) between 25 and 30 km, replicating the procedure used by G-66. In Fig. 3 the differences ΔN_d(z) for 66 soundings at Nantucket and the 9 for Fairbanks are plotted. For Lexington, ΔN_d(z) values are both negative and positive, but higher values of N_dUS(z) dominate. For Fairbanks N_dUS(z) always is greater.

The errors in lidar retrievals of ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$, attributed to the use of Temp(z) and Pr(z) from a model atmosphere to retrieve N_d(z), are of the order of 3 % and decrease to 1 % when the source of Temp(z) and Pr(z) are soundings (Russell et al., 1979). Again in Table 2, the magnitudes of the absolute differences between the US 1962 Standard Atmosphere and the soundings at Lexington and Fairbanks for α_a(532,z) are of the order of 3 %, agreeing with the error attributed if case models are used instead of soundings to derive β_m(λ,z).

Figure 5 Contour plots of the vertical profiles of 532 nm aerosol extinction ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$ at Lexington, showing values uncorrected and corrected for two-way transmittance, ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z)$. The red stars indicate the dates on which the measurements were conducted.

Figure 4 shows ${\mathit{\tau }}_{{\mathrm{a}}_{*}}$ both for Lexington (blue stars) and for Fairbanks (red diamonds). The means for the entire period of measurements available at each site are 0.0215 and 0.0099, respectively. Also shown is a monthly mean τ_a for the Northern Hemisphere (Sato et al., 1993). The mean ${\mathit{\tau }}_{{\mathrm{a}}_{*}}$ at Fairbanks is half that at Lexington, providing evidence of the decreasing aerosol amount with increasing latitude. Because of the variability in ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$, ${\mathit{\tau }}_{{\mathrm{a}}_{*}}$ values from Lexington vary widely from the Fairbanks mean to the Sato et al. (1993) magnitude, the current reference for this period. However, as we will see in the next section, better agreement is found when the measurements are corrected with two-way transmittance attenuation.

Taking into account the small difference between the results using the US 1962 Standard Atmosphere and the soundings to derive β_m(λ,z), the first simpler option can reliably be used. However, we decided to use the soundings to minimize the errors and to capture the more realistic features of the aerosol cloud.

Table 3 The same as Table 2 but for the comparison of ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$ vs. ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z)$ and ${\mathit{\tau }}_{{\mathrm{a}}_{*}}$ vs. ${\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}$. See text for details.

3.2 Aerosols extinction contour plots and optical depth corrected by aerosol two-way transmittance attenuation

Figure 5 shows the contour plots of ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{\ast }$ for uncorrected and corrected two-way transmittance, ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z)$, for Lexington. The initial values of TAOD were used to obtain a first estimate of ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{{*}_{\text{t2w}}}$. This ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z{)}_{*}$ is only used to calculate sAOD for each day and is subtracted from TAOD to produce the tropospheric-corrected value (tAOD); the calculation is repeated to determine new profiles of the two-way aerosol transmittance and correct ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{\ast }$ generating ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z)$. Panel (a) shows the contour plot of uncorrected values of ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{\ast }$ and panel (b) the contour plot of ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z)$. The magnitudes of ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z)$ are higher than ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{\ast }$. The two-way transmittance correction is dominated by the aerosols, in particular the tropospheric aerosols. The maximum extinction is at the third maximum, 1.071 $\times {\mathrm{10}}^{-\mathrm{2}}$ km⁻¹ located at 15 km, on 27 March 1965. Similarly in Fig. 6, the Fairbanks contour plots for ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{\ast }$ and ${\mathit{\alpha }}_{\text{a}}^{T_{\text{a}}}(\mathrm{532},z)$ show a notable difference in magnitude. Here the absolute maximum extinction occurred on 16 August 1964 at 15 km with a magnitude of 3.8 $\times {\mathrm{10}}^{-\mathrm{3}}$ km⁻¹.

Figure 6 Contour plots of the vertical profiles of 532 nm aerosol extinction, ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$, at Fairbanks, showing values uncorrected and corrected for two-way transmittance, ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z)$. The red stars indicate the dates on which the measurements were conducted.

Table 3 contains the relative and absolute means and maximums for $\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}$, $\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}\mathit{\%}$, $\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}$, and $\mathrm{\Delta }{\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}\mathit{\%}$ calculated using Eqs. (14) to (17) respectively but for ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{\ast }$ vs. ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z)$ and ${\mathit{\tau }}_{{\mathrm{a}}_{*}}$ vs. ${\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}$. The magnitude of $\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}$ produced by the two-way transmittance correction is of the order of 10⁻³ km⁻¹ for Lexington and 10⁻⁴ km⁻¹ for Fairbanks or an increase of 67 % and 26 %, respectively. These increases are due mainly to the two-way aerosol transmittances dominated by the tropospheric AOD with magnitudes more than twice as high at Lexington than at Fairbanks. The increase in magnitude reveals more details of the vertical distribution of the ${\mathit{\alpha }}_{\text{a}}^{T_{\text{a}}}(\mathrm{532},z)$ and in the case of Lexington the presence of a fourth maximum during May 1964 whose vertical location matches the decreasing trend at the core of the stratospheric aerosol cloud.

Figure 7 Time series of stratospheric AOD (sAOD₅₃₂) for Fairbanks for τ_a(532,z) and ${\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z)$.

In Figs. 7 and 8 the increases in ${\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}$ with respect to ${\mathit{\tau }}_{{\mathrm{a}}_{*}}$ for Fairbanks and Lexington are shown. At Lexington the ${\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}$ magnitudes are approximately the values of τ_a from Sato et al. (1993) for the Northern Hemisphere (black line). This agreement is an important confirmation of the Sato et al. (1993) magnitudes for τ_a from Agung in the Northern Hemisphere. Again in Table 3, the magnitudes of the increase in ${\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}$ are of the order of 10⁻² for Lexington and 10⁻³ for Fairbanks, representing 66 % and 26 % increases, respectively.

At Lexington the absolute maximum value of ${\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}$, 0.071, occurs on 30 March 1965, 3 d after the absolute maximum extinction was registered at 15 km. At Fairbanks the absolute maximum value of ${\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}$, 0.018, was registered on 16 August 1964, the same day the absolute maximum extinction was registered at 15 km at this site.

Since the pioneering lidar work by Fiocco and Grams (1964), multiple researchers have contributed to the development of the processing algorithms to retrieve aerosol optical properties and errors (Russell et al., 1979; Klett, 1981, 1985; Kovalev, 2015). These works explain the limitations of retrieving the full set of optical variables characterizing the stratospheric aerosols from the dataset of Fiocco and Grams (1964). However, assuming a Junge size-distribution model and Mie scattering with a refractive index of 1.5, Fiocco and Grams did produce estimates of the aerosol content of the stratosphere at 16 km: number concentration, surface area, and the aerosol density per unit volume of air. They also used the mean profile they derived to estimate the total particles per cubic centimetre, total surface area, and a total mass, integrating the concentrations obtained between 12 and 24 km (GF-67). The only available optical property estimates, based on some of the cited particle concentration estimates at 16 km and in the column, are the aerosol extinction at 16 km ($\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ km⁻¹) and the aerosol optical depth (0.015), both at 694 nm (Deirmejian, 1971).

For comparing with the values reported above, we estimate α_a(694,z) from α_a(532,z), as well as τ_a(694,z) from τ_a(532,z), using the wavelength exponents for aerosols from Mt. Pinatubo in the range of wavelengths from 532 to 694 nm (Jäger and Deshler, 2002). We made the estimates for both Lexington and Fairbanks as no clear assignation of the values to either site is made in G-66 or GF-67. At 16 km, the mean value of α_a(694,z) was 10⁻³ km⁻¹ for Fairbanks and $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ km⁻¹ for Lexington, matching the order of magnitude estimated by Deirmejian (1971).

From 1963 to mid-1965, in addition to the 1963 Mt. Agung eruption, two other volcanoes were reported to have erupted in the Northern Hemisphere with Volcanic Explosivity Index (VEI) 3. They were Trident Volcano in Alaska at 58^∘ N and 155^∘ W and the Vestmannaeyjar volcano (also known as Surtsey) south of Iceland at around 63^∘ N and 20^∘ W. The first was reported to have erupted in April 1963, and its plume reached 15 km (Decker, 1967). The second remained in eruption between November 1963 and February 1964, with its plume in November 1963 reaching an altitude around 4.5 km above the tropopause more than once, located approximately at 10.5 km (Thorairinsson, 1965). They were attributed as contributing to the replenishment of aerosols in the mid-latitude lower stratosphere, following the increase in the atmospheric turbidity determined using twilight measurements (Cronin, 1971).

Twilight measurements revealed three peaks in atmospheric turbidity between the March 1963 Agung eruption and the end of 1965 shown in Fig. 1 from Volz (1970). The first turbidity peak in that figure with the highest magnitude was registered by the end of 1963, when no lidar measurements were available, but its decaying is seen in the sAOD during the first half of 1964 in our Fig. 4. The second turbidity peak, having approximately the same magnitude as the third, is located in the last months of 1964, coincident with the second sAOD peak in Fig. 4. The third turbidity peak also coincides with the third sAOD peak. Updated information reveals the extension in time of the Vestmannaeyjar, from late 1963 to the middle of 1964 (GVP, 2013a), and the occurrence of two additional eruptions of Trident Volcano, the first between 17 October and 17 November 1963 and the second on 31 May 1964 (GVP, 2013b), all of them at VEI 3. That sustained input of the aerosols in the Northern Hemisphere stratosphere explains the second and third peaks with similar magnitudes in the turbidity (Fig. 1 in Volz, 1970, and in the sAOD in our Fig. 4).

Figure 8 Time series of stratospheric AOD (sAOD₅₃₂) for Lexington for τ_a(532,z) and ${\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z)$ compared to that from the NASA Goddard Institute for Space Studies volcanic aerosol dataset (Sato et al., 1993).

3.3 Relative errors

Table 4 reports the results for the estimated relative errors in the aerosol extinction with and without the aerosol two-way transmittance correction for both sites. In addition, the relative errors of the backscattering ratio and aerosol backscatter at 694 nm and the aerosol backscatter at 532 nm are reported. The relative errors for ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}\le \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ km⁻¹ were excluded in the statistics.

Note the increases in the mean relative errors from $\left(\frac{\mathit{\delta }\mathrm{SR}}{\mathrm{SR}}\right)$ to $\left(\frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{a}}}{{\mathit{\beta }}_{\mathrm{a}}}\right)$, 12 % to 48 % at Fairbanks and 13 % to 36 % at Lexington, the higher increases occurring during the full processing. It is explained by the factor ${\left(\frac{{\mathit{\beta }}_{\mathrm{m}}}{{\mathit{\beta }}_{\mathrm{a}}}\right)}^{\mathrm{2}}$ in Eq. (18). Because the processing algorithm relies on Eq. (6) to derive β_a from β_m, the squared ratio will be lower than 1 if SR< 2, increasing as SR decreases and reaching the value ${\left(\frac{{\mathit{\beta }}_{\mathrm{m}}}{{\mathit{\beta }}_{\mathrm{a}}}\right)}^{\mathrm{2}}={\mathrm{10}}^{\mathrm{4}}$ for $\mathrm{SR}=\mathrm{1.01}.$ Only with SR≥2 is the ratio lower than 1, which in the case of Fairbanks happens at one level on 1 d. In the case of Lexington, 10 % of the SR values are higher than 2.

Table 4 Relative error estimates of the backscattering ratio, aerosol backscatter at 694 nm and 532 nm, and aerosol extinction with and without correction for aerosol two-way transmittance at 532 nm for Lexington and Fairbanks. Errors for ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}\le \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ km⁻¹ were not included in the statistics. All errors are percentages.

In Table 4, the second highest increase in the mean relative error occurs in the calculation of $\left(\frac{\mathit{\delta }{\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}}{{\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}}\right)$ from $\left(\frac{\mathit{\delta }{\mathit{\alpha }}_{\mathrm{a}}}{{\mathit{\alpha }}_{\mathrm{a}}}\right)$. At Fairbanks the increase is 7 %, from 54 % to 61 %, and at Lexington the increase is 20 %, from 44 % to 64 %. The error is associated with the magnitudes of the relative errors from $\left(\frac{\mathit{\delta }{T}_{\mathrm{a}}}{T_{\mathrm{a}}}\right)$, conducted at this step for the reasons explained above. At Fairbanks the mean value of $\left(\frac{\mathit{\delta }{T}_{\mathrm{a}}}{T_{\mathrm{a}}}\right)$ is 8 %, while it is 44 % at Lexington, associated with the expression δτ_a=0.5τ_a. It should be taken into account that the total AOD at both sites is dominated by the magnitude of the tropospheric AOD, which is higher at Lexington.

Figure 9 (a) Contour plot of relative error estimates for Lexington. (b) Contour plot of the aerosol extinction at 532 nm corrected by two-way transmittance for Lexington . Note the two data gap periods greater than 1 month for the consecutive measurements: March, and July to September both in 1964. They are identified with vertical dotted red lines at the 7th and 23rd measurements. In the top panel the areas in white in the relative error contour plot represent relative errors for ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}\le \mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ km⁻¹. They were not included in the statistics in Table 4.

Figure 10 The same as Fig. 9 but for Fairbanks.

The time vs. altitude contour plots of the $\left(\frac{\mathit{\delta }{\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}}{{\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}}\right)$ relative errors and of the ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(z,n)$ are shown in Figs. 9 and 10 for Lexington and Fairbanks, respectively. The regions with maximum magnitudes of ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}$ at both sites are associated with the lower relative errors as expected. At Lexington, for ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}>\mathrm{8}\times {\mathrm{10}}^{-\mathrm{3}}$ km⁻¹ the relative errors are ≤ 30 %. It is also evident that relative errors equal to or lower than 50 % dominate both in time and altitude. In the case of Fairbanks, for ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}>\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ km⁻¹ the relative errors are ≤ 40 %. The relative errors of ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}$, in Table 4, produce ${\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}$ relative errors above 100 %. Those estimated values of the relative errors for ${\mathit{\tau }}_{\mathrm{a}}^{T_{\text{a}}}$, together with the ones in Table 4, are substantially larger than other sets of volcanically perturbed stratospheric aerosols lidar measurements.

The high error magnitudes in the $\left(\frac{\mathit{\delta }{\mathit{\beta }}_{\mathrm{a}}}{{\mathit{\beta }}_{\mathrm{a}}}\right)$ at 694 nm estimation could be reduced in the case that SR values increase. In several of the 75 SR profiles a renormalisation processing could increase SR magnitude. This is reasonable since the normalization altitude range (no aerosol present) was 25 to 30 km, where there certainly would be some aerosol present. An inspection of the plots of SR vs. altitude in Figs. 14, 15, and 16 in G-66 shows the presence of aerosols between 25 and 30 km. Furthermore, in some of the profiles SR is above 1 at all levels (1.0 indicates no aerosol). In addition, the introduction of the two-way transmittance correction in the processing of SR will increase SR from the raw returned lidar signal.

Options are available to find the raw lidar data to conduct the reprocessing described above. These include searching for the filmed images of the oscilloscopes used as registers and/or the original punch cards (probably transferred to tapes) both reported in G-66. A last resort would be the digitalization of the SR from the figures cited above. The original signal profiles could then be reconstructed by inverting the normalization procedure applied to produce the SR profiles.

3.4 Attribution of the 1963 Agung aerosol cloud within the Lexington lidar dataset

In this section, we seek to understand whether some of the sAOD variations observed by the Lexington lidar may originate from sources other than the March 1963 Agung eruption (such as the two stratosphere-injecting 1963 VEI 3 eruptions discussed in Sect. 3.2: Trident, Alaska, and Vestmannaeyjar, Iceland). Specifically, we compare the Lexington extinction dataset to four different model-based volcanic forcing datasets for the Agung aerosol cloud. Three of the four Agung forcing datasets are from two different interactive stratospheric aerosol models: two different SO₂ emissions scenarios from the UM-UKCA model (Dhomse et al., 2020, 2021) and a third simulation from the 2D-AER model (Arfeuille et al., 2014), as applied within the CMIP6 volcanic aerosol dataset (Luo et al., 2016). The fourth simulation is from an idealized model representation of the Agung cloud, based on a simple parameterization for the progression of the tropical reservoir of volcanic aerosol, and its dispersion to mid-latitudes (Ammann et al., 2003), used to represent historical volcanic forcings in some CMIP5 climate model historical integrations (see Driscoll et al., 2012).

Figure 11 Interactive stratospheric aerosol model (Dhomse et al., 2020) representations of the Agung aerosol cloud sAOD₅₅₀ compared with that of the Lexington dataset, the Ammann et al. (2003) volcanic forcing dataset, and the CMIP6-AER2D volcanic aerosol dataset (Luo, 2016).

The progression of volcanic aerosol clouds from major tropical eruptions reaching the stratosphere was established by Dyer et al. (1970, 1974). They synthesized the extensive set of observations on the Agung aerosol cloud (Dyer and Hicks, 1968) and used the analyses of the global dispersion of radionuclides from Pacific thermonuclear tests in the 1950s (e.g. Machta and List, 1959). The continual slow upwelling circulation in the tropics and the sub-tropical barrier at the edge of the tropical pipe combine to cause the long-lived tropical stratospheric reservoir (Dyer, 1974; Grant et al., 1996), which is the reason why tropical eruptions have such prolonged radiative cooling compared to mid-latitude eruptions. The Brewer–Dobson circulation (Brewer, 1949; Dobson, 1956) has a strong seasonal cycle, transporting air preferentially towards the winter pole, causing an increasing mid-latitude sAOD trend during autumn and a decreasing mid-latitude sAOD trend during spring (in both hemispheres). Each of the model lines in Fig. 11 shows this circulation-driven seasonal variation in sAOD, with the transport of the Agung aerosol remaining in the tropical reservoir predicted to increase during October and November, reaching a peak in January to March in both 1964 and 1965. The model-predicted variations are consistent with the initial observed sAOD values of 0.04–0.05 in January and February 1964 being higher than most of the 0.01–0.04 sAOD values observed in October and November 1964, and the expected variations from the models suggest the suddenly higher sAOD values ∼ 0.05 may be from a different source than Agung. However, whereas sAOD values would be expected to increase going into winter, the December, January, and February sAOD at Lexington are mostly lower than during the autumn, which indicates an additional source of stratospheric aerosol may have continued to add to the Agung cloud sAOD throughout the autumn of 1964. Furthermore, the 1965 Lexington observations show a continuing increase in sAOD into the springtime, whereas the models predict the sAOD from Agung would have reduced by a factor of 2 during the first 6 months of 1965. The analysis suggests that another source of sAOD influential during this period (either the two VEI 3 volcanic eruptions in 1963/1964 or some other source of material into the stratosphere) must have caused the observed increase in stratospheric AOD during 1965 with a potentially substantial influence also during autumn 1964.

Figure 12 compares the vertical structure of the 9 Tg representation of the Agung aerosol cloud from Dhomse et al. (2020) at 42^∘ N, compared to the Lexington observations, confirming that these model simulations capture the altitude of the cloud during the early period of the measurements (January to May 1964). However, although the magnitude of the simulated aerosol extinction compared well with the original Lexington dataset (Dhomse et al., 2020), with the two-way transmittance corrections applied here, the 9 Tg simulation is low-biased compared to the lidar measurements, even in this earlier period, suggesting the 12 Tg UM-UKCA aerosol simulation (not shown) would likely compare better (Dhomse et al., 2020, 2021). None of the four model-generated Agung forcing datasets can explain the observed increase in extinction during January to July 1965. The sudden peaks in April and June 1965 have quite a different vertical structure compared to the early 1964 measurements, the sAOD in 1965 having a substantial component from the altitude range 18–20 km. This vertical profile analysis again suggests the episodic sAOD enhancements in spring 1965 were from a different source than the 1964 measurements.

Figure 12 Contour plots of ${\mathit{\alpha }}_{\mathrm{a}}(\mathrm{532},z{)}_{*}$ from Dhomse et al. (2020) at 42^∘ N and measured aerosol extinction corrected for two-way transmittance ${\mathit{\alpha }}_{\mathrm{a}}^{T_{\text{a}}}(\mathrm{532},z)$ from lidar for Lexington.

4 Data availability

The datasets of the original rescued backscattering ratios and the calculated aerosol backscatter (both at 694 and 532 nm) and the aerosol extinction at 532 nm (both corrected and uncorrected for two-way aerosol transmittances) at Lexington and Alaska are available at https://doi.org/10.1594/PANGAEA.922105 (Antuña-Marrero et al., 2020a).

5 Summary

We have carried out a data recovery of the first ever multi-year lidar dataset of the stratospheric aerosol layer, the Lexington and Fairbanks measurements profiling the portion of the Agung volcanic aerosol cloud that dispersed to Northern Hemisphere mid-latitudes and high latitudes, respectively. The results show the high level of variability in the stratospheric aerosol extinction for Lexington between January 1964 and July 1965 that is mainly attributed to the 1963 Mt. Agung eruption. At Lexington the highest aerosol extinction values and aerosol optical depths (1.1 $\times {\mathrm{10}}^{-\mathrm{2}}$ km⁻¹ and 0.076, respectively) were actually observed at the end of March 1965 in the final phase of the 1.5-year-long record. Based on contemporary and updated reports of additional volcanic eruptions in the Northern Hemisphere between 1963 and 1965, we tentatively suggest a potential explanation for the apparently contradictory temporal sAOD trends to be the VEI 3 eruptions of Vestmannaeyjar, Iceland, and Trident, Alaska. Further research combining observational data and modelling should be conducted to elucidate the individual contributions from each of those eruptions to the stratospheric aerosol layer at this location of the Northern Hemisphere.

The level of the relative errors are unusually high considering that under high loads of volcanic aerosols in the stratosphere, the signal to noise ratio is high in the returned lidar signal. The analysis of the contributions of the variables along the different steps of the processing algorithm identified the two main sources of error. The main one, accounting for a little more than 30 % of the relative error, is associated with the division of the molecular backscatter by the aerosol backscatter, which is directly linked to low magnitudes of the backscattering ratio. Those low magnitudes are produced by two factors: the first is the lack of two-way transmittance corrections in the backscattering ratio calculation from the raw squared distance-corrected signal, and the second is that the normalization altitudes, considered to be empty of aerosols, were too low and actually did contain aerosol. We suggested alternatives to search for the original signal profile records or to reconstruct the original signal profiles from the plotted backscattering ratio records, including the normalization region from 25 to 30 km. A future search for original records should also take into account the 34 missing files from the 100 referred to by Fiocco for Lexington (Waymouth et al., 1966).

In general the results reported should be considered as the first estimates. We report the comparison of the aerosol extinction values and aerosol optical depths that we calculated with information available up to the present, showing reasonable results. Improvements in the two factors cited above lead to an increase in magnitude of the aerosol extinction and optical depth in several of the profiles.

We have also compared the Lexington sAOD time series to four different model representations of the 1963 Agung aerosol cloud and illustrated how the model predictions suggest that the sAOD above Lexington from Agung must have decreased from January to July 1965, whereas the 1965 lidar observations show a clear increase in sAOD through the spring into summer. The UM-UKCA Agung aerosol simulations suggest the Agung cloud likely descended to a lower altitude in 1965 than in 1964. In contrast, the lidar measurements show more sudden aerosol extinction enhancements reaching up to 20 km in altitude during 1965. Contemporary records of two VEI 3 high-latitude eruptions (in Alaska and Iceland) suggest their volcanic clouds reached the stratosphere in both cases, and model comparisons strengthen the potential attribution of the January to July 1965 sAOD₅₅₀ increase to a source other than the 1963 Agung eruption.