Interactive comment on “ 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

Somewhere, it should be mentioned that the primary quantity measured by a stratospheric lidar (and also produced in the dataset) is the backscatter ratio or the aerosol backscattering coefficient, not the extinction coefficient. Extinction is a derived / secondary quantity. It relies much more on assumptions (about the extinction to backscatter rato, also called lidar ratio) than backscatter. Extinction is usually small, but is necessary to derive the best possible backscatter profile.

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.pangaea. de/10.1594/PANGAEA.922105 (Dataset in Review) . 45

Introduction:
The abrupt enhancements to the stratospheric aerosol layer from historical large magnitude volcanic eruptions (e.g. Deshler, 2008) cause substantial radiative forcings to the Earth's climate system, and reducing their uncertainty remains an important priority for international scientific research, volcanic forcings being the strongest driver of natural climate 50 variability (e.g. Hansen, 1978;Robock, 2000). One of the co-ordinated 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 1982and Mt. Pinatubo in 1991. The main motivation for this HErSEA multi-model experiment (Historical Eruption SO2 Emission Assessment) is to gather stratospheric aerosol observations to evaluate the model 55 simulations, and understand the current diversity in the sulphur emission amount and altitude distribution interactive stratospheric aerosol models use when simulating the Pinatubo aerosol cloud (see section 3.3.2 of Timmreck et al., 2008).
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), with another recent study focusing on assessing the variability and global distribution of the Agung aerosol cloud (Niemeier et al., 2019). 60 Whereas the models participating in ISA-MIP simulate volcanic aerosol clouds interactively, the historical climate model simulations that provide the main basis for attributing past climate variability (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), reference aerosol optical properties used to enact volcanic surface cooling. The observational data constraining the Agung aerosol cloud in both the interactive models and for the volcanic forcing datasets has hitherto tended to be based on 65 column optical properties measured at the surface, 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 data set of Agung backscatter ratio 70 or aerosol extinction has yet been a vailable to the scientif ic community. Whereas the Jamaica lidar (Clemesha et al., 1966) also measured the Agung cloud, the first multi-year data set of lidar measurements of the volcanic aerosol from that eruption was conducted from Lexington, Massachusetts from January 1964 to August 1965 (Grams and Fiocco, 1967, hereinafter GF-67). No d igital record of these lidar measurements existed until now, the data apparently only presented in Figures of the lidar backscattering ratio profiles within published scientific papers, providin g only quantita tive information 75 about the altitude of the Agung aerosol cloud. Although the descent in the peak of the backscatter ratio profile from Lexington is analysed within GF-67, only limited estimates of the cloud's aerosol extinction exist (2 x 10 -3 km -1 at 16 km and the aerosol optical depth of 0.015 (Deirmejian, 1971) were been produced.
We discovered however, after initial search of digital archives at several institutions, that the original lidar backscattering ratio profile measurements from the Lexington and Alaska 1964/5 soundings are fully tabulated in Gerald W. Grams PhD 80 thesis conducted under the supervision of Prof. Giogio 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 or 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 recognised open-access data archive. https: //doi.org/10.5194/essd-2020-246 Open Access Earth System Science Data Discussions Preprint. Discussion started: 3 November 2020 c Author(s) 2020. CC BY 4.0 License. Dhomse et al. (2020) used preliminarily processed lidar data from Lexington, MA, one of the two sites reported in  to compare model aerosol extinction at 16 km with lidar observations, finding good agreement. They also noted the large differences between the CMIP5 and CMIP6 volcanic aerosol datasets for the Agung and El Chichón periods, pointing out the importance of reducing this uncertainty by reconcilin g the datasets with additional stratospheric aerosol observations.
There only an init ial (preliminary) single-level version the Lexington aerosol extinction dataset was used, our analysis here completing the processin g of the full vertical profile (between 12 and 24km), with t wo -way transmittance corrections 90 applied to the aerosol backscatter ratio, aerosol extinction and optical depth datasets, also with a detailed and transparent assessment of the relative error in each metric.
This work is a contribution to the Data Rescue activity of the Stratospheric Sulfur and its Role in C limate (SSiRC) a SPARC init iative (SSiRC, 2020), following a recent dataset recovery of two ship-borne lidar datasets that measured the progression of the highly uncertain "tropical core" of the Pinatubo aerosol cloud in July 17 th to September 13 th 1991, 4-12 95 weeks after the 15 th June 1991 Pinatubo eruption (Antuña -Marrero et al., 2020b). Those datasets were identified prio rity for the data recovery, being in is the period when the Stratospheric Aeroso ls and Gas Experiment II (SAGE II) satellite could only observe the stratospheric aerosol above around 22 km due to the extreme opacity of the aerosols (McCormick and Veiga, 1992). 100 2. Materials and Methods.

Lidar instrumentation:
The first successful laser radar ranging experiment was conducted at the Research Laboratory of Electronics, Massachusetts 105 Institute of Technology, at Lexington, Massachusetts, and consisted of analysing the return signal from a very high frequency (nano-second) laser and for the 60-140km altitude range (Smullin and Fiocco, 1962). The research team , led by Prof. 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 110 improvements, stratospheric aerosols were detected 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 figures 3 and 4 of G-66 respectively. There are listed also the main features of the lidar instruments used for the measurements at Lexington and Co llege, Alaska, reported in its  An additional set of relevant features of the instrument follows. The fluorescent emission after the laser pulse has been emitted, was prevented incorporating a small rotating shutter into the transmitting unit, synchronized with the Q-switchin g device. The sensing unit for the backscattered signal consisted in an astronomical telescope, with an interference filter and a photomultiplier synchronized to another rotating shutter to avoid its exposition to the intense returns from short distance s. 120 The photomultiplier was cooled with methanol and dry ice, to reduce the levels of its dark current (G-66).

Lidar measurements:
Lidar observations were conducted at Lexington,Massachusetts (42° 25'N,71° 15'W) and also at College, (64° 53'N, 125 148°3'W) located in the city of Fairbanks, Alaska, hereinafter identified as Fairbanks. The measurements were supported by the NASA . One of the semi-annual reports mention more than 100 measurements conducted (Fiocco, 1966a). However the amount of total profiles appearing in Grams PhD dissertation was 75. Nine days of measurements from July 26 to Au gust 21, 1964, were conducted in Fairbanks. At Lex ington, Massachusetts, 23 days of measurements from January 14 to May 20, 1964, and 43 days from October 11, 1964to July 21, 1965 were made. At both 130 sites measurements were conducted at night to avoid the contribution of daylight to the signals registered by the photomultipliers.
The lidar signal returns at both sites were registered photographically from oscilloscopes covering up to 40 km and then digitized. Then the digitized lidar return signals from a set of daily laser shots were averaged in 1 km bins (G-66; GF-67).

Backscattering ratios in the original lidar dataset:
It is well known that solving the equation for the single wavelength elastic lidar is an ill-posed problem. The single returned signal, is the resu lt of two main species, molecules and particles, making it necessary to the use of additional information to estimate the solution (eg, Kovalev, 2015). That is the reason why still today processing single wavelength lidar profiles 140 remains a challenge. Considerin g this fact we may understand the magnitude of the challenge confronted by Prof. Girogio Fiocco and then BSc Gerald W. Grams, Prof. Fiocco's Ph D student, when they conducted the processing of the first ever set of lidar returned signals from stratospheric aerosols.
We now describe the procedure applied in G-66 to derive SRo(z). The average photoelectron flux registered by the photomultiplier was represented by the expression: 145 Where z is the altitude, n z is the number of photons at the altitude z, NA is the molecular number density at altitude z, obtained from the US 1962 Standard Atmosphere. K is a constant resulting from all the terms not depending on the altitude in the optical radar equation, includin g T 2w 2 , the two-way atmospheric transmittance (see G-66 for more details). The assumption of a constant value for T 2w 2 was based on the atmospheric attenuation model proposed by Elterman (1964). The 150 model provided magnitudes of the molecular and aerosol scattering and the ozone absorption, showing that almost all attenuation of the laser bean occurs in the troposphere. The model allowed an estimate at 700 nm of the variability of the term T 2w 2 between 10 and 30 km which was belo w 3%. The correct ion of the returned signal, associated with the two way transmittance of the laser beam throughout the atmosphere, was then neglected and it was assumed that the atmospheric 155 attenuation term was constant.
The returned 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 the calculated in each profile between 12 and 24 km. To that end, for each profile the average value between 25 and 30 km of the ratios calculated in the former step were determined. Then for each profile the ratios in the altitude range 12 and 24 km 160 were d ivided by the average value of the ratios between 25 and 30 km from the same profile. The resu ltin g values were considered to be the backscattering ratio (SRo(λ, z)): the ratio between the total (aerosols + molecules) backscattering divided by the molecular backscattering. The normalization procedure assigned the backscattering ratio to be equal to one above 25 km, after assuming the contribution from aerosols was negligible compared to the molecular at those levels. This assumption would lead to an under-estimate of stratospheric aerosol since there would have been aerosol at these altitudes 165 (Russell, et al., 1979).
The SRo(λ, z) from the lidar measurements conducted at Lexington and Fairbanks were reported in tabular format in the Gerald W. Grams PhD Thesis (G-66), and cited in the acknowledgements section of GF-67. It was the unique reference of its existence, the clue that guided us in our search for the lidar measurements.

Algorithm and complementary datasets used in the processing:
Far beyond the mere rescue and deposit on public repositories of datasets, mainly from stratospheric aerosols from past volcanic eruptions from the 60's to the present, the SSiRC Data Rescue Activity is committed, whenever it will possib le, to re-calibrate each dataset and determine its levels of uncertainties (SSiRC, 2020). Because some stratospheric aerosol lidar 175 datasets have already been identified and located, we consider that its recalibration or reprocessin g should be conducted using a standardized a lgorithm, to guarantee the best possible consistence among the different lidar datasets. To contribute to this task we describe below the processing algorithm we used as a first step in that direction.
The lidar backscattering ratio (SR(λ, z)) is defined a s 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 attenuated backscatter 180 (β a A ( λ, z ) ). That definition is associated to the fact that in the retrieval of SR o ( z ) the two-way total transmittance(T T 2 ) correction was neglected (Hosteler et al., 2006): βm(λ, z) is derived using the equation: where Sm = (8π/3)kbw is the molecular extinction to backscatter ratio for the molecular scattering, commonly approximated by 8π/3 (Collins and Russell, 1976) after neglecting the dispersion of the refractive index and the King factor of the air represented by kbw. The volume molecular scattering coefficient, σm(λ, z) is determined by the equation: Where NA = 6.02214×10 23 (1/mol) is Avogadro's number; Ra = 8.314472 (J/K/mol) is the gas constant and Qs(λ) the total 190 molecular scattering cross section per molecule for the standard air. The derived equation for Qs( λ) for standard air is (Hostetler et al., 2006): Then from equation 2 we retrieved β a A ( λ, z ) , using the expression: To derive the aerosol backscatter profiles at 532 nm ( β a ( 532 , z ) ) we used the wavelength exponents (kb(z, t)) for aerosol backscatter in the range of wavelengths between 532 and 694 nm derived for the stratospheric aerosols produced by the https://doi.org/10.5194/essd-2020-246 Applying in addition the corrections for, T m 2 ( 694, z ) and T O3 2 ( 694, z ) the two-way molecular and ozone transmittances at λ 200 = 694 nm to the β a A ( 694 , z ) . The generic definition of the two-way tansmittance is: With the sub index j representing in α (λ, z) the vertical profiles of the scattering by the molecules α (λ, z)), ozone (α O3 (λ, z)) and aerosols (α a (λ, z)).
Finally we derive the aerosols scattering corrected by the total two -way transmittance (α a Ta ( 532, z ) ) applying the correct ion by the two-way aerosol transmittance T a 2 ( z ) 210 α a Ta ( 532, z ) = α a ( 532 ,z ) T a 2 (532 ,z) Because the information available to calculate the T a 2 ( 532, z ) should be determined from total aerosol optical depth (TAOD) measurements from sun photometers we applied a two -step procedure. The first step consists of using the TAOD to calculate a first guess T a 2 ( 532, z ) followed by the calculation of a first guess α a Ta ( 532, z ) * profile. Then the stratospheric AOD (sAOD) is calculated integrating α a Ta ( 532, z ) * between 12 and 24 km. The second step calculates: 215 and the calculation of T a 2 ( 532, z ) is repeated but using tAOD instead of TAOD in equation (11) getting the definitive values of α a Ta ( 532 , z ) .
The algorithms for the solution of the single wavelength lidar equations apply the two -way transmittance correction to the raw lidar returned signal, together with squared distance correction, well before the backscattering ratio is calculated. In 220 our case the available information we have are the backscattering ratios which have been derived without conducting the two way transmittance correction (G-66). That is the reason that correction was included in the retrieval o f β a ( 694, z ) in equation (7). However only the molecular and ozone two way transmittance corrections ( T m 2 ( 694, z ) , T O3 2 (694, z)) were included in this step.
The aerosol two-way transmittance correction, T a 2 ( 532, z ) , was delayed until the final step to derive α a Ta ( 532, z ) . The 225 reason was that the available information on the AOD was at λ = 500 nm and it was the total AOD, includin g the sAOD that we are trying to retrieve. No Ån gstrom exponent contemporary information for the Agung eruption in the eastern US in the range 500 to 694 nm was found. Using the Ån gstrom exponent climatological values covering the cited wavelength range from 1995 to 2019 from the nearest Aerosol Robotic Network (AERONET, 2020) stations we derived , we converted the total AOD at 500 nm to 532 nm only, because it was a very near λ. It was the solution to lack of trusted information 230 with the aim to minimize the error that could be introduced for converting the AOD to 694 nm. There have been abundant accounts about the changes of the physical-chemical properties of aerosols in the eastern US from the sixties until the present (Went, 1960;Husar et al., 1991).

Complementary datasets used: 235
The correction for the attenuation of the lidar signal by the two -way transmission is often considered negligib le and ignored, based on signal to noise ratio considerations or for simplicity as it was the case in the original p rocessin g of the se set of measurements (G-66; GF-67). We were motivated to make that correction by the fact that, during a little more than https://doi.org/10.5194/essd-2020-246 Preprint. Discussion started: 3 November 2020 c Author(s) 2020. CC BY 4.0 License. half a century, the accuracies of the different instruments available for measurements of the stratospheric aerosols from the 240 1963 Mt. Agung eruption have been under, still unsettled, scrutiny and discussions (ex., Deirmendjian, 1965;Dyer, 1971a;Clemesha, 1971;Dyer, 1971b;Deirmendjian, 1971;Stothers, 2001;Timreck et al., 2018). Our goal was to produce the most accurate aerosol extinction and optical depth from the rescued measurements, based on the contemporary state of the art mea surements in the sixties of the XX century. Table 2 summarize the locations of the sites where radiosonde, ozone soundings and atmospheric turbidity measurements 245 were conducted. Also the distances from each individual site to the corresponding lidar site a re provided. Following each individual dataset is described. To derive βm(λ, z), αm(λ, z), αa(λ, z) and T m 2 (z) the required variables are the temperature (Temp(z)), pressure (Pr(z)) and molecular number density (Nd(z)). We used the vertical profiles of Temp(z), Pr(z) and Nd(z) from the 1962 US Standard 255 Atmosphere (U. S. Standard Atmosphere, 1962).
We also made use of the Temp(z) and Pr(z) profiles (deriv ing Nd(z)) from the most complete and nearest sounding station.
In the case of the soundings we took into account the fact that lidar observations were performed at night, typically near Preprint. Discussion started: 3 November 2020 c Author(s) 2020. CC BY 4.0 License.

Datasets used for the estimation of the ozone 2-way transmittance:
We used the NO3(z) from the 1966 US Standard Atmosphere Supplement (COESA, 1967). In addition, we used the seasonal means of NO3(z) between 1963 and 1967 from ozone soundings conducted at L. G. Hanscom Fla., Bedford, MA and Fairbanks, AK for Lexington and Fairbanks respectively (Hering and Borden, 1967). 280 The profile of the ozone absorption coefficient at a given wavelength ( k O3 ( λ, z ) ) is calculated using the profile of ozone cross sections (σ O3 (λ, Temp ( z )) ): at the temperature (Temp(z)), where NO3(z) is the number density of ozone. The σ O3 ( λ, Temp(z) ) at λ = 694 nm is provided by Serdyuchenko et al., (2014) in the temperature range 193 to 293 ºK. We used the average of σ O3 ( λ, Temp(z) ) , 285 9.88e-22 cm 2 molecules -1 , considering that the standard deviation of this averaging profile represents 2.4 % variability of the average value. This set of absorption coefficients have been recommended by the recent status report from the International Ozone Commission from WMO (Orphal et al., 2016).  Flowers et al., (1969). The curve of the monthly mean B belonging to B lue Hill Observatory in the frame plot on the figure were d igit ized (WebPlotDigit ize, 2020). Then TAOD at λ = 500 nm was calculated using the equation (Volz, 1969), result ing from converting the decadic lo garithm used to define B to the neperian logarithm used for AOD: For Fairbanks we did not find contemporary measurements, but there were manually conducted measurements in several places in the Arct ic and Antarctic, including at Fairbanks (64.86°N, 147.85 °W, 133 m), with sun photometers at several wavelengths (Shaw, 1982). Those measurements are reported to be corrected by the molecular scattering and gas absorption (Shaw et al., 1973). The instruments were calibrated at Mauna Loa Observatory using Langley method with root mean square errors (RMSE) of sun photometer voltage output readings (V) of δV V ≈ 10 −3 having a systematic RMSE for 305 AOD = +/-0.002 and the total error estimated as +/-0.004 (Shaw, 1982). At Fairbanks the annual mean AOD = 0.110 from 105 observations at λ = 500 nm is reported on table 1 of the cited reference, but also the AOD annual cycle appears in the lower panel of figure 2, showing the high AOD values in late winter and spring, peaking up to 0.135. We then digitized the mean AOD and its variation range values for July (no data for August appears o n the figure), result ing AOD = 0.082 ± 0.022. Although Shaw (1982) does not provide the information of the year the measurements were conducted this data is 310 cited and cited to have been conducted in 1978 by Freund (1983).
We also used AOD data at 500 nm from the 2 nearest AERONET stations to each site having long-term records. Bonanza Creek, AK, is less than 30 km from the location the lidar measurements were reported to be conducted at Fairbanks. Ångstrom exponents calculated from the TAOD climatological monthly means for the interval 440 to 870 nm ,that we used for Lexington. In the case of Bonanza Creek, we had one "contemporary" value from July 1978 (Shaw, 1982), and we 320 selected the July climatological monthly mean, as the "current" value. After the conversion to 532 nm they were respectively 0.087and 0.242, a nd we used the same value for both July and August 1964 lidar measurements at Fairbanks.
In the case of Lexington, for comparison purposes, we d igitized the average monthly mean TAOD for 26 Eastern US stations from 1972to 1975from Husar et al., (1981. The series of monthly mean TAOD values were converted from 500 nm to 532 nm, with the procedure described above, using the Ån gstrom exponents for the interval 440 to 870 nm. The

Numerical and statistical methods:
For each of the two datasets we calculated the statistical resu lts of the differences ( Δα * US ) between α a ( 532, z ) US calculated using the same β m ( 594 , z ) profile from the 1962 US Standard Atmosphere for all the days and the α a ( 532, z ) * calculated using the β m ( 594 , z ) profiles derived from the daily soundings: 345 Δα a * = α a ( 532, z ) US − α a ( 532, z ) * and the percent differences Δα% a * by the expression: Similarly we defined the differences Δα at2w and the percent differences Δα at2w % between the α a ( 532, z ) * calculated using the β m ( 594 , z ) profiles derived from the daily soundings, and its corrected values α a ( 532, z ) t2w resu lting for accounting 350 for the two-way atmospheric transmittance.

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

Backscattering ratio relative error: 365
As it was explained above, the data we rescue are a reasonable approximation of what we today know as the backscattering ratio described in equation (2). Then we use the equation (19) from Russell (1979)  Where SR ( λ, z ) is the total backscattering ratio; N s is the signal measured; T 2w the two-way transmittance from aerosols; molecules and ozone; β m the molecular backscattering; β m * molecular backscatter at the normalization level; SR ( λ, z ) total backscattering ratio at the normalization level and C FF * 2 the covariance between measured β m and β m * . = 3 % for both sites (Russell et al, 1979). In addition we assumed ( The term δSR min was evaluated according to table (1b) in Russell (1979) for the SRmin = 1.01 and the respective latitudes of both sites. Then for both sites, according to Russell et al. (1979).
385 2.7.2 Aerosol backscattering relative errors: The equation (18) in Russell (1979) The estimated error for the 2-way transmission corrections in Russell et al. (1979) provides the expression: 390 and considering the standard error of determinations of τ a , τ O3 , and τ m are respectively 50, 20 and 10% the following estimates are produced. That is: δτ a = 0.5 τ a , δτ O3 = 0.2 τ O3 and δτ m = 0.1 τ m . However, in our calculus chain of β a only the ozone and molecular two-way transmittances were used.
For this section of the procedure we considered ( δβ m β m ) = 10% because we used radiosonde soundings at both sites 395 (Russell et al, 1979). We neglected the error in computing Qs usin g equation (5)  Next we determined the relative error in β a ( 532, ) associated with the conversion from β a ( 694, ) in equation (7), using the wavelength exponents (kb(z, t)) for aerosol backscatter in the range of wavelengths between 694 and 532 nm (Jäger and

Aerosol extinction relative errors:
In the case of the α a , its relative errors are: 405 The last term in the right side represents the error in the EBC for λ = 532 nm. In the case of the ones we used (Jäger and Deshler, 2002;2003) the error has been estimated in ± 40 % according to Deshler et al., (2003). For α a Ta , the aerosols extinction corrected by the aerosols two-way aerosols transmittance, using the est imates of its relative error described above: 410 Using the cited set of equations and the assumptions described above we evaluated the error for e ach level in each measurement. Figure 2 shows the α a ( 532, z ) cross-sect ions for Lexington. Panel a) α a ( 532, z ) US is calculated using the same β m ( 594, z ) profile from the 1962 US Standard Atmosphere for all the days; b) α a ( 532, z ) * was calculated using the daily β m ( 594 , z ) 420 profiles derived from the sounding at Nantucket, MA. On top of the figures we p lotted the dates the measurements were conducted (red starts at 24.5 km level). In the case of Lexin gton the two data gaps higher than 1 month, March and July to Regarding the magnitudes of α a ( 532, z ) US in f igure 2, they are slight ly higher than the ones from α a ( 532 , z ) * . That is also 435 the case in figure 3 showin g the cross-sections for Fairbanks, with panels similar to figure 2. Th is is quantified in table 3.

Aerosols extinction cross sections and optical depth:
At both sites the mean and maximum values for Δτ a * and Δα * are positive showin g that the magnitudes of α a US and τ a US are in general higher than α a * and τ a * . Also in the table we appreciate that the magnitudes of the mean percent difference increase of both variables is around 1%.
The fact described above disagrees with the possibility G-66 mentions about lower aerosol backscatter from the retrieval 440 they conducted, using the 1962 US Standard Atmosphere, and the more realistic ones using soundings. He arrived to 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 (Herin g and Borden, 1965). He reported N d (z) variability rarely exceeded 5% of the N dUS (z) values at altitudes between 10 and 30 km. 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: The values in the denominators M dUS and M d are the mean values of N dUS (z) and N d (z) between 25 and 30 km 450 respectively, rep licating the procedure used by G-66. In figure 4 the differences ΔN d ( z ) for all the 66 soundings at Nantucket used to calculate N d (z) and the 9 for Fairbanks are plotted. For Lexin gton, on panel a), N dUS (z) values are both negative and positive, but higher values of N dUS (z) dominate. It is confirmed that the relative means and the maximum values of Δα a * between the α a ( 532, z ) US and α a ( 532, z ) * for Lexington in table 3 are the same order of magnitude, 10 -5 km -1 for the relative and absolute means and 10 -4 km-1 for the maximum, larger than for Fairbanks. The values of the 455 relative means Δα * % confirm the higher values when the 1992 US Standard Atmosphere is used, in contradiction with G-66 estimation.
It has been stablished that the errors in lidar retrievals of α a ( 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 3, the magnitudes of the absolute differences between the US 1962 Standard 460 Atmosphere and the soundings at Lexington and Fairbanks for α a ( 532, z ) are in the order of 3%. That magnitude matches the error attributed if case models are used instead of soundings to derive β m ( λ, z ) .
The figure 5 shows τ a * both for Lexington (b lue stars) and for Fairbanks (red diamonds). Also figure 5 shows the monthly mean τ a for the northen hemisphere (Sato et al., 1993). The means for the entire period of measurements available at each https://doi.org/10.5194/essd-2020-246 Lexington are around the magnitude of the mean τ a * at Fairbanks, because of the variability of α a ( 532 , z ) * . Few τ a * values from Lexignton have magnitudes near the values of Sato τ a , the current reference for this period. However, as we will see in the next section a better agreement is found when the measurements are corrected by two way transmittance attenuation.
Taking into account the little d ifference between the results using the US 1962 Standard Atmosphere or the soundings to 470 derive β m ( λ, z ) , the first simpler option could 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.      Figure 6 shows the cross-sections of α a ( 532 , z ) * for uncorrected and corrected two-way transmittance (α a Ta ( 532, z ) ) for

Aerosols extinction cross-sections and optical depth corrected by aerosol two-way transmittance attenuation:
Lexington. The initial values of TAOD were used to obtain a first estimate of α a ( 532 , z ) * 2 . This α a Ta ( 532, z ) * is only 495 used to calcuate sAOD for each day and is subtracted from TAOD to produce the tropospheric corrected value (tAOD) and the calculation is repeated to determine new profiles of the two -way aerosol transmittance and correct α a ( 532, z ) * generating theα a Ta ( 532 , z ) . Panel a) shows the cross-section of uncorrected values of α a ( 532, z ) * , in panel b) the crosssections of α a Ta ( 532 , z ) . The magnitudes of α a Ta ( 532, z ) are higher than α a ( 532, z ) * . The t wo-way transmittance correction is dominated by the aerosols, in particular the tropospheric aerosols. The maximum extinction is at the third maximum, 500 1.071 x 10 -2 km -1 located at 15 km, on March 27 th 1965. Similarly in figure 7, the Fairbanks cross-sect ions for α a ( 532, z ) * and α ( 532, z ) show a notable difference in magnitude. The absolute maximum extinction occurred on August 16 th 1964 at 15 km, with a magnitude of 3.8 x 10 -3 km -1 .  (14) to (17) respectively but for α a ( 532, z ) * vs α a Ta ( 532, z ) and * vs . The magnitude of produced by 505 the two-way transmittance correction is in the order of 10 -3 km -1 for Lexington and 10 -4 km -1 for Fairbanks. They represent an increase of 72 % in the first case and 26 % in the second. 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  Durin g the course of more than two decades after the pioneering stratospheric aerosols measurements with lidar work by Fiocco and Grams (1964) multiple researchers contributed to the development of the processing algorithms to retrieve aerosols optical properties and its errors (Russell et al, 1979, Klett, 1981Klett, 1985, Kovalev, 2015. Those facts explain the limitations that do not allow the retrieval of the full set of optical variables characterizing the stratospheric aerosols 525 from the Fiocco and Grams dataset. However using a Junge size -distribution model, and assuming Mie scattering with refractive index 1.5, they produced 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 use the mean profile they derived to estimate the total of particles/cm 3 , 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 in some of the cited particle concentration estimates at 16 km and in 530 the column are the aerosol extinction at 16 km, 2 x10 -3 km -1 and the aerosol optical depth of 0.015, both at 694 nm (Deirmejian, 1971).
For comparing with the values reported above, we made estimates at of α a ( 694 , z ) from α a ( 532, z ) as well as τ a ( 694 , z ) from τ a ( 532 , z ) usin g the wavelength exponents for aerosols from Mt Pinatubo in the range of wavelengths 532 to 694 nm (Jäger and Deshler, 2002). We made the estimates both for Lexin gton and Fairbanks because no clear assignation of the 535 values cited above to one of the two sites is made in G-66 and GF-67. At the 16 km level the mean value of α a ( 694 , z ) was 10 -3 km -1 for Fairbanks and 2 x 10 -3 km -1 for Lexington matching for both sites the order of magnitude estimated by Deirmejian, (1971).
An additional validation of those results, in particular for τ a Ta ( 532, z ) at Lexington appears in figure 9, where the stratospheric τ a ( 532, z ) for the northern hemisphere from January 1964 to July 1965 has b een plotted (Sato et al., 1993). 540 The magnitude of τ a Ta ( 532 , z ) is the same at Lexington (and also at Fairbanks, figure 8) as the τ a ( 532, z ) from Sato et al., (1993).
From 1963 to mid-1965, in addition to the 1963 Mt Agung, t wo other volcanoes were reported to h ave erupted in the northern hemisphere with Volcanic Explosivity Index (VEI) 3. They were the 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 and. The first 545 was reported to have erupted in April 1963 and its plume reaching 15 km (Decker, 1967). The second remained in in eruption between November 1963 and February 1964, with its plume reaching more than once in November 1963 an altitude around 4.5 km above the tropopa use, located approximately at 10.5 km (Thorairinsson, 1965). They were attributed https://doi.org/10.5194/essd-2020-246 contributing to the replenish ing of aerosols in the mid-latitude lower stratosphere, following the increase of the atmospheric turbidity, determined using twilight measurements (Cronin 1971). 550 Twilight measurements revealed 3 peaks in atmospheric turbidity, between the March 1963 Agung eruption in and the end of 1965 shown in figure 1 from Volz (1970). The first turbidity peak in that figure with the highest magnitude was registered by the end of the 1963, when no lidar measurements were available, but its decaying is seen in the sAOD durin g the first half of 1964 on our figure 5. The second turbidity peak, having approximately the same magnitude than the third, is  Volz (1970) and in the sAOD in our figure 5. 560   Table 5 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 ≤ 5 x 10 -4 km -1 were excluded in the statistics. 580

Relative Errors: 575
Note the increases in the mean relative errors from ( ) to ( ) , 12 % to 48 % at Fairbanks and 13 % to 36 % at Lexington, the higher increases occur durin g the full p rocessin g. It is explained by the factor ( ) in equation (21).
Because the processing algo rithm relies on equation (6) to derive from the squared ratio will be lower than 1 if < 2, increasing as decreases and reaching the value ( ) = 10 4 for = 1.01. On ly with => 2 is the ratio is lower than 1, which in the case of Fairbanks happens at one level on one day. In the case of Lexington, 10% of the are 585 higher than 2.
In table 5, the second high increase in the mean relative error happened in the calculation of ( ) from ( ) . At Fairbanks the increase is 7%, from 54% to 61%. At Lex ington the increase is 20% from 44% to 64%. The error is associated with the magnitudes of the relative errors from ( ) , conducted at this step by the reasons explained above. At Fairbanks the mean value of ( ) is 8% wh ile 44% at Lexington, associated to the expression δτ a = 0.5 τ a . It should be 590 taken into account that the total AOD at both sites are dominated by the magnitude of the tropospheric AOD, which is higher at Lexington. Preprint. Discussion started: 3 November 2020 c Author(s) 2020. CC BY 4.0 License. Table 5: Relative error estimates of the backscattering ratio, aerosol backscatter at 694 nm and 532 nm, aerosol extinction with and without correction for aerosol two -way transmittance a t 532 nm for Lexington and Fairbanks.

Errors for
≤ 5 x 10 -4 km -1 were not included in the statistics. 595 The . As expected at both sites the regions with maximum magnitudes of are associated with the lower relative errors. In figure 10 note that at Lexington, for > 8 x 10 -3 km -1 the relative errors has a magnitude equal or lo wer than 30%. It is also ev ident that relative errors equal or lo wer than 50% dominate both in t ime and altitude. In the case of Fairbanks, figure 11, for > 2 x 10 -3 km -1 the relative error has a magnitude equal or lower than 40%.
Considering the magnitudes of the relative errors for in table 5 it is ev ident that the relative errors are above 605 100%. Those estimated values of the relative errors for together with the ones in table 5 show h igh magnitudes compared with other sets of volcanically perturbed stratospheric aerosols lidar measurements.
As explained above, the highest error introduced in the ( ) at 694 nm estimation could be reduced if the have higher values. In several of the 75 profiles renormalization processing could increase its magnitude. That is possible, because the normalization procedure applied, considered that above 24 km no aerosols were present. Inspection of the plots 610 of vs altitude in figures 14, 15 and 16 in G-66 shows the presence of aerosols between 25 and 30 km and in some of the cases at all of those levels magnitude is above 1, the value representin g no aerosols.
In addition, what will definitely increase the magnitude of , will be the introcution of the two-way transmittance correction in the processing generationg from the raw returned lidar signal.
In the search for the raw lidar data several options are available. Searching for the filmed images of the oscilloscopes used 615 as registers and/or the original punched cards (probably transferred to tapes) both reported in G-66. The last resort would be the digitalization of the from the figures cited above. Then the original signal profiles could be reconstructed inverting the normalization procedure applied to produce the profiles.

630
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 VEI3 discussed in section 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 forcin g datasets are from two different interactive stratospheric aerosol models: two different SO2 emissions scenarios from the UM-UKC A 635 model (Dhomse et al., 2020) and a third simulation from the 2D-AER model (Arfeuille et al. 2014), as applied with in the CMIP6 volcanic aerosol dataset (Luo et al., 2016). The fourth simulations is from an idealised 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).  The progression of volcanic aerosol clouds from major tropical eruptions reaching the stratosphere was established by Dyer et al. (1970;1974) from analysing the extensive synthesis of observations on t he Agung aerosol cloud (Dyer and Hicks, 645 1968), and from knowledge derived from analyses of the global dispersion of radionuclides from Pacific thermonuclear tests in the 1950s (e.g. Machta and List, 1959). The continual slo w upwellin g 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 compare d to mid-latitude eruptions. The Brewer -Dobson circu lation (Brewer, 1949;Dobson, 1956) has a strong seasonal cycle, transporting air 650 preferentially towards the winter pole, causing an increasing mid-latitude sAOD trend during autumn and a decreasing midlatitude sAOD trend during spring (in both hemispheres). Each of the model lines in Figure 12 show this circulation -driven seasonal variation in sAOD, with the transport of the Agung aerosol remaining in the tropical rese rvoir 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 655 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 continu ed to add to the Agung cloud sAOD throughout the autumn of 1964. Furthermore, the 1965 Lexington observations show a 660 continuing increase in sAOD into the springt ime, whereas the models predict the sAOD from Agung would have reduced by a factor of 2 during the first 6 months of 1965. Th e analysis suggests another source of sAOD influential durin g this period (either the two VEI3 volcanic eruptions in 1963/4 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. 665 Figure 13 compares the vertical structure of the 9Tg representation of the Agung aerosol cloud from Dhomse et al. (2020) at 42N, 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 to the original Lexington dataset (Dhomse et al., 2020), with the two-way transmittance https://doi.org/10.5194/essd-2020-246 corrections applied here, the 9Tg simulation is low-b iased compared to the lidar measurements, even in this earlier period, 670 suggest ing a with the 12Tg UM-UKCA simulation likely to compare better (not shown). None of the 4 model-generated Agung forcing datasets can explain the observed increase in extinction during Jan to July 1965, with the sudden peaks in April and June 1965 having a quite different vertical structure to the early 1964 measurements, the sAOD in 1965 having a substantial component from the altitude range 18-20km. This vert ical profile analysis again suggests the ep isodic sAOD enhancements in spring 1965 were from a different source than the 1964 measurements, the altitude of the peak extinction 675 in that first year of the dataset broadly consistent with the UM-UKCA simulations of the Agung cloud. The level of the relative errors are unusually high considering that under high loads of volcanic aerosols in the stratospher e, the signal to noise ratio is high in the returned lidar signal. The analysis of the cont ributions of the variables along the 695 different steps of the processing algorithm, allowed identifying 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, directly linked to low magnitudes of the backscattering ratio. Those lo w 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. The second is the normalization method conducted in the region considered to be empty 700 of aerosols, when in many profiles the signal plots reveal its presence. We suggested alternatives to search for the original signal profile records or to reconstruct the original signal profiles from the plotted backscattering ratio records, includin g the normalization region from 25 to 30 km. The search for original records should include looking for the at least 25 missing profiles from the total of at least 100 Fiocco mentions.
In general the results reported should be considered as the first estimates. We report the comparison of the aerosol 705 extinction values and aerosol optical depths 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 t imeseries to 4 d ifferent model representations of the 1963 Agung aerosol cloud, and illustrate how the model predictions su ggest the sAOD above Lexington from Agung must have decreased from 710 January to July 1965, whereas the 1965 lidar observations show a clear increase in sAOD through the spring into summer.