Retrieval of aerosol backscatter, extinction, and lidar ratio from Raman lidar with optimal estimation

Abstract. Optimal estimation retrieval is a form of nonlinear regression which determines the most probable circumstances that produced a given observation, weighted against any prior knowledge of the system. This paper applies the technique to the estimation of aerosol backscatter and extinction (or lidar ratio) from two-channel Raman lidar observations. It produces results from simulated and real data consistent with existing Raman lidar analyses and additionally returns a more rigorous estimate of its uncertainties while automatically selecting an appropriate resolution without the imposition of artificial constraints. Backscatter is retrieved at the instrument's native resolution with an uncertainty between 2 and 20%. Extinction is less well constrained, retrieved at a resolution of 0.1–1 km depending on the quality of the data. The uncertainty in extinction is > 15%, in part due to the consideration of short 1 min integrations, but is comparable to fair estimates of the error when using the standard Raman lidar technique. The retrieval is then applied to several hours of observation on 19 April 2010 of ash from the Eyjafjallajokull eruption. A depolarising ash layer is found with a lidar ratio of 20–30 sr, much lower values than observed by previous studies. This potentially indicates a growth of the particles after 12–24 h within the planetary boundary layer. A lower concentration of ash within a residual layer exhibited a backscatter of 10 Mm−1 sr−1 and lidar ratio of 40 sr.


Introduction
Aerosols impact the Earth's radiation budget both directly, by reflecting solar radiation back into space (Haywood and Shine, 1995), and indirectly, by altering the properties and distribution of clouds (Lohmann and Feichter, 2005) or reacting with other species (Colbeck, 1998). The lack of knowledge about the global distribution and composition of aerosols is currently the single greatest source of uncertainty in estimates of net radiative forcing and therefore is a factor in the ability to predict the impacts of climate change (IPCC, 2007).
Lidar (light detection and ranging) is an active remote sensing technique for observing the distribution of molecules and particles in the atmosphere as a function of height by means of the light they backscatter from a laser beam. The name intentionally emulates radar as both techniques use the time of flight of a pulsed source to deduce the distance to the scatterer (Fugii and Fukuchi, 2005). Despite its exceptionally high spatial and temporal resolution, lidar is not as widely applied as other techniques in the study of aerosol's effect on climate as a single lidar samples only one location. With the launch of a space-based lidar (Winker et al., 2009) and the development of networks across northern America (Welton et al., 2000), Europe (Pappalardo et al., 2005), and Asia (Sugimoto et al., 2008), the importance and availability of lidar data increases.
The energy observed by a lidar is a function of the extinction and backscattering coefficients -the cross section per unit volume to either attenuate the beam or to scatter light directly back towards the instrument. These coefficients are functions of the microphysical properties of the aerosol present, such as refractive index and size distribution. Deriving such properties directly is possible, but the problem is very poorly constrained. Its solution either requires a greater number of measurements, such as a multi-wavelength system (Müller et al., 1999), or further assumptions about the scatterers. This more complex problem is disregarded here in favour of the better-constrained estimation of extinction and backscatter.
Optimal estimation retrieval is a form of nonlinear regression which determines the most probable circumstances that produced a given observation, weighted against any prior knowledge of the system. For several decades, it has been successfully applied to the analysis of satellite (e.g. Marks and Rodgers, 1993;Li et al., 2008;Watts et al., 2011), radar (Grant et al., 2004), and ground-based radiometer observations (e.g. Guldner and Spankuch, 2001) but has not seen substantial use within the lidar community. This paper applies the technique to the estimation of aerosol extinction and backscatter from two-channel Raman lidar observations. The retrieval processes the entire profile simultaneously, making optimal use of the information available and choosing the most appropriate vertical resolution for the result while fully characterising the covariant error due to measurement noise, model error, and other assumptions. This widely recognised retrieval algorithm, which is less dependant on ad hoc corrections and assumptions while providing rigorous error estimation, provides an alternative to the traditional Raman lidar technique that can make optimal use of the available information.
Section 2 outlines the retrieval algorithm and existing analysis methods. Section 3 evaluates the retrieval's performance against existing techniques with simulated data and considers the error budget. Section 4 applies the algorithm to observations, while Sect. 5 provides some conclusions.

Optimal estimation retrieval
As outlined in Rodgers (2000), optimal estimation solves the inverse problem where y is a column vector describing the measurements; gives the noise on that measurement; and the forward model F (x, b) translates a state of the instrument and atmosphere, summarised by unknown parameters x and known parameters b, into a simulated measurement. Approximating the probability density function (PDF) for all quantities as Gaussian and using Bayes' theorem, the probability that the system has a state x given the measurement y can be written as where the covariance matrix S describes the random experimental error and x a is the a priori, the state expected before the measurement is made. The uncertainty in that expectation is described by the a priori covariance S a . The quantity −2 ln P(x|y) is hereafter referred to as the cost as it measures the goodness of fit for a solution. Good models should have a cost approximately equal to the number of measurements. Hence, the cost will herein be quoted normalised by the length of y.
It can be shown that the iteration converges to the most probable statex, where the Jacobian K i = ∂F /∂x i and i is a scaling constant. General practice, outlined in Fig. 1, is that if the cost increases after an iteration, i is increased by a factor of 10. Otherwise, it is reduced by a factor of 2. Evaluation ceases after -the cost function decreases by much less than the number of measurements; -the cost decreased and the change in the state is much less than the predicted error, where the error covariance matrix of the solution S = (K T i S −1 K i + S −1 a ) −1 ; -the step in state space, is much less than the length of the state vector; -30 iterations, which is considered a failure to converge.
The averaging kernel, describes the extent to which the true and a priori states each contribute to the solution as it can be shown that where a hat indicates the value after convergence. An ideal retrieval would have a kernel equal to the identity. In practice, the rows of A will be peaked functions showing how the information in one retrieved bin is derived from an average of the true values around it. The width of that peak is therefore a measure of the resolution of the retrieval. . Schematic of the optimal estimation retrieval algorithm.
36 Fig. 1. Schematic of the optimal estimation retrieval algorithm.

Existing lidar analyses
The number of photons observed from a height R is expressed by the lidar equation (Measures, 1992) where β(R, λ) is the coefficient for incident laser light, wavelength λ L , to be backscattered at a wavelength λ; α(R, λ) is the extinction coefficient; E L is the number of photons emitted by the laser pulse; and C(R), known as the calibration function, describes the alignment and efficiency of the detection system. As both the extinction and backscatter are unknown, a single profile presents an underconstrained measurement. The atmosphere is assumed to contain only two components such that where β (p) , α (p) denote backscattering and extinction by aerosols and β (m) , α (m) denote scattering by molecules, which is well modelled by Rayleigh scattering.
The dominant return for any lidar is the elastic profile (where λ = λ L ), from which the backscatter is commonly derived by a technique known as onion peeling or the Fernald-Klett method (Klett, 1981;Fernald, 1984). A Raman lidar monitors a second channel containing the Raman scattering from a single species in the atmosphere, such that β becomes a known function of number density. Ansmann et al. (1992) outlined a means to invert Eq. (6) in such circumstances to derive the extinction and backscatter separately.
A few applications of nonlinear regression to lidar have been published, though these retrieve the microphysical properties of aerosol rather than the optical properties outlined above or later in this paper. A retrieval of ice water path and effective radius in cirrus clouds from coincident, space-borne lidar and radar measurements was developed in Delanoe and Hogan (2008), though its results were found to be highly dependant on the microphysical assumptions. Pounder et al. (2012) derived high-quality extinction retrievals from three simultaneous observations with different fields of view using a linearised model of the lidar equation that included multiple scattering while applying Twomey-Tikhonov smoothing rather than an a priori. Marchant et al. (2010) presented an original, if limited, linearised scheme that decomposed scattering over a basis of precomputed aerosols. This was expanded to a retrieval of effective radius in multiwavelength studies via a Kalman filter in Marchant et al. (2012).
A related method known as regularisation has also been used to derive extinction and backscatter from Raman lidar profiles. The introduction of Veselovskii et al. (2002) provides a good review of early attempts and the methodology. Shcherbakov (2007) and Pornsawad et al. (2012) demonstrated that such methods return solely positive extinction and can produce more accurate products than the Ansmann method but use Tikhonov smoothing within the retrieval, which generates significant errors where there are substantial gradients. Though the formulation of that technique resembles that of optimal estimation, they differ in their interpretation. As the Tikhonov matrix H is singular, the smoothing term H T H has no inverse that would correspond to a S a . Hence, the desired smoothing is introduced in a manner that has no physical analogue within the measurement system. Evaluation also requires a set of basis functions to be defined, artificially imposing a structure onto the system.
It is preferable to impose the basis for smoothing solutions through an a priori covariance matrix derived from actual data and the physical processes driving the system, as facilitated by optimal estimation retrieval. The impact of these assumptions can be assessed through the averaging kernel, such that it is clear where the data are the dominant influence on the solution. Other techniques do not provide such a balance and, in fact, rarely discuss the choice of basis functions or their impact on the solution.

Forward model
Lidars frequently use photomultiplier tubes (PMTs) as their detector, which produce a cascade of electrons when struck by a photon. If the rate of photons is less than two per bin, noise can be very effectively removed by applying a threshold to the output to return a count of the number of photons per bin. As count rates increase, multiple pulses are more likely to overlap and be counted only once, such that this mode becomes increasingly nonlinear (Müller, 1973). The most frequent correction for this (Whiteman et al., 1992) assumes that after any count, the detector will be unable to detect another for a "dead time" τ d and that that time is independent of the count rate (known as the non-paralysable correction). The observed profile can then be expressed as where subscript is used to denote a function of range, f (R i ) ≡ f i ; M is the number of laser shots E is averaged over; and τ b is the duration of a bin (such that R i = 1 2 icτ b ). This will apply to both channels, though each may have a distinct value of τ d .
For large count rates, it is also possible to operate the PMT in an analogue mode, which simply averages the output current during each range bin. This is linear over a large dynamic range but suffers additional noise from thermal excitations, variations in pulse height, electrical interference, and other effects. In such circumstances, where a and b are constants and ϕ is measured in volts (whereas it is photon counts in Eq. 9). The correction of Eq. (9) is increasingly unreliable for count rates greater than 10 MHz. When both detection modes are operated simultaneously, it is common practice to "glue" the two signals together using analogue observations near the instrument, where the signal is largest, and photon counting elsewhere (e.g. Newsom et al., 2009). The gluing procedure in essence finds appropriate values for the constants in Eq. (10). Such techniques will not be used in this paper as the two signals will be considered separately. A simple extension of this work would be to include both Eqs. (9) and (10) in the forward model and retrieve using both signals simultaneously.
It is not necessary to retrieve the extinction and backscatter at the native resolution of the instrument. The state vector can be defined on any arbitrary grid, r, and then interpolated onto the instrument's range axis R. For example, a grid with spacing that increases as a function of height could be used such that all retrieved values have uncertainties of a similar magnitude. The use of a coarser grid will also reduce the computational expense of the retrieval. For simplicity, only a regular 33 m grid is considered in this paper, though the general expressions are presented to most accurately represent the model used in the calculations. No significant difference has been found between interpolated solutions and those at full resolution.
Neglecting multiple scattering and assuming α (p) ∝ λ −1 (Ansmann et al., 1992), the number of photons observed from range bin R i will be where superscript denotes functions of wavelength, f (λ X ) ≡ f (X) ; σ R is the cross section for Rayleigh scattering, which has lidar ratio B = 8π/3; N is the atmospheric number density; N i = R i 0 N (R )dR ; and E B is the background count rate, which is estimated from observations as R → ∞. The aerosol optical thickness χ = R 0 α (p) (λ L , R )dR and β are evaluated at λ L , though this dependence is dropped for brevity. The calibration function C is assumed known and is input as a parameter.
A tilde is used to represent variables on the retrieved grid r, which are interpolated using the cubic spline method of Press et al. (1992) onto the measured grid R. The aerosol optical thickness is evaluated on grid r with the trapezium rule, and then interpolated onto R. Note that the extinction is assumed constant through the first bin, such that it acts as a boundary term rather than a physically meaningful value. This avoids various difficulties with observation very near the instrument. Though this problem could be linearised, there will be some error involved in that approximation. Solving the nonlinear problem presented is not an overly intensive calculation and so there is no need to simplify the problem.
Extinction and backscatter are both functions of N and so will be correlated. This should be identified within S a but cannot be easily estimated. Further, the use of correlated variables will emphasise degenerate states of the forward model, which can slow the retrieval's convergence. This can be averted by retrieving the lidar ratio B = α (p) /β (p) instead (Shcherbakov, 2007;R. Hogan, personal communication, 2012), which is independent of N, All elements of x should be positive or there will be a degeneracy in the impact of backscatter and optical depth on the elastic channel, which can impede the retrieval. For the retrieval of β (p) and α (p) , this will be prevented by setting all negative values to zero after evaluating Eq. (3). Though an unsophisticated solution, it is found that values that are zeroed in one iteration generally increase in the next. It is unusual that a pixel retrieves a negative value in the final iteration without also having a large error. For the retrieval of β (p) and B, it was found preferable to instead retrieve ln β (p) while enforcing a lower limit of unity on B. Retrieving the logarithm eliminates negative values from the solution but can contort state space, retarding convergence.
The measurement and state vectors are then The first guess for x in the iteration Eq.
(3) is taken as β (p) = 10 −5 Mm −1 sr −1 and B = 58 sr. These values were chosen as they tend to reduce the number of iterations. Their value does not affect the final result (if the retrieval converges).
One final note must be made of the treatment of measurement error (which is assumed uncorrelated), S . The observed photon counts should be Poisson distributed, such that their variance is equal to their mean. This is widely used to justify approximating the variance of a lidar measurement with the measurement itself. This is not strictly valid as the measurement is only a single sample of a distribution. However, a lidar sums profiles over several seconds or minutes of laser shots during data collection, giving no further measure of their variance.
The optimal estimation scheme requires an unbiased estimate of the variance. Using the measurement itself causes the retrieval to favour observations that coincidentally suffer large, positive noise as they are they appear to be more precise. This effect is most pronounced at low signal levels and introduces a high bias into the retrieval. To alleviate this, the variance will be estimated by the application of a fivebin, sliding-window average through the data. The impact of this will not be explored in detail, though preliminary studies found that even minimal smoothing of the variance vastly reduced biases. A more involved statistical analysis could be applied, considering the correlation of adjacent bins, but is not deemed necessary.
For further details, derivations, and justification of the forward model, please consult Chapter 2 of Povey (2013).

A priori
Arguably the most important component of an optimal estimation scheme is its a priori. Ideally, the a priori would not greatly affect the retrieval, but in practice it constrains which states are deemed to be both physically possible and likely.
In this problem, solutions should be reasonably smooth as aerosols are often well mixed through the planetary boundary layer (PBL), but gradients should not be completely excluded as layering does occur.
The exact composition and optical properties of aerosol are highly variable such that a detailed a priori is unlikely to be representative of the broad range of potential states (Lopatin et al., 2013). Climatologies, from which an a priori would be derived in most applications, rarely exist and most lidar data are derived from elastic instruments and so will be influenced by the assumed lidar ratios. Therefore, the most appropriate a priori would be an order-of-magnitude estimate to weakly constrain the magnitude of the state vector while prescribing vertical and inter-variable correlations.
Some generalised descriptions of representative aerosol types have been explored in the literature. The OPAC model (Optical Properties of Aerosol and Cloud; Hess et al., 1998) defines size distributions, refractive indices, and number densities for a variety of cloud and aerosol particles. These provide the necessary inputs for Mie codes  to calculate the extinction and backscatter. Combinations of these based on expected and observed compositions then produce characteristic aerosol mixtures, such as marine or urban.
For the data to be considered, the continental type should be appropriate -comprising soot with soluble and insoluble aerosols. An ensemble of scattering properties was constructed by randomising the abundance of these components, using the OPAC model values as the mean of a Gaussian distribution with width estimated by 10 % of that mean (the exact value assumed was found to be unimportant). A simple treatment of aspherical particles using the T-matrix code of Dubovik et al. (2006) produces an effectively identical distribution. The resulting distributions of β (p) , α (p) , and B are shown in Fig. 2. The a priori is based on qualitative fits to these, shown in blue when retrieving linearly and in red for a logarithmic retrieval, where the values are given in Table 1. Though the logarithmic retrievals appear to give a better fit to the distributions, the retrieval of ln β (p) and ln α (p) was found to be overly constrained. Though the lidar ratio is theoretically a better description of the state, its distribution is not symmetric and so not necessarily well suited to optimal estimation. A relatively broad a priori distribution has been selected to compensate. These distributions demonstrate an approximately 95 % correlation between β (p) and α (p) , which is included in S a for the linear retrieval. Though its exact value appears to be unimportant, it would be desirable to obtain a more rigorous estimate.  Table 1. The title of each plot gives a scale factor for its vertical axis.
The OPAC model states that the density of non-dust aerosols decreases exponentially with a scale height of 2 km. The prescribed values will therefore decrease similarly in x a . Further, there will almost certainly be some vertical correlation of the measurements due to vertical mixing. The simple model of a Markov process proposed by (2.83) of Rodgers (2000) shall be used with correlations decaying exponentially with separation, where H is a scale height. The suitability of this scale height can be assessed by investigating the covariance of some measure of aerosol scattering. A convenient option is backscatter sondes (NDACC, 1989(NDACC, -2000, which measure the light backscattered from a xenon flashlamp approximately every 30 m during a balloon ascent (Rosen and Kjome, 1991;Rosen et al., 2000). Profiles over 10 years of observations at three sites have been used in Fig. 3 to produce a correlation matrix of backscatter ratio with height. Its rows decay roughly exponentially with height, which when fitted to Eq. (15) give H = 1-2 km in the free troposphere, consistent with the model.  This will not necessarily apply within the PBL, which is only weakly coupled to the free troposphere (Oke, 1987). Several studies of the vertical distribution of aerosol within the PBL have been performed with tethered balloons, though the data could not be readily accessed. These generally find that aerosol concentrations are constant with height (Figs. 10-12, 4, and 2 of Ewell et al. (1989), Greenberg et al. (2009), andFerrero et al. (2010), respectively), but occasionally observe fine structure ( Fig. 6 of Ferrero et al., 2011). Further, a myriad of literature covers observations of aerosol layers tens to hundreds of metres thick (e.g. Di Girolamo et al., 1999;Dacre et al., 2011) or variations within lofted aerosol features (e.g. Althausen et al., 2000;Huang et al., 2010).
A rigorous a priori covariance matrix would represent both the general tendency for aerosols to be well mixed throughout the troposphere and the fine-scale structure that occasionally occurs. The average position of the top of the PBL would be expressed by a significant decrease in correlation between areas above and below it. At the moment, there is insufficient information to quantify these effects with any degree of certainty. As such, a conservative estimate of H = 100 m is used here, which will not make the best use of the available information but does not overconstrain the solution.

Simulations
Simulated data can be easily produced with the forward model, using the NOAA (1976)   PBL extinction profile is modelled by an error function (Steyn et al., 1999) multiplied by an exponential decay above the PBL. Aerosol and cloud layers are modelled by Gaussian peaks (G. Biavati, personal communication, 2011). An analytic model outlined in Povey et al. (2012) is used to generate the calibration function and detector nonlinearity. Once simulated, Poisson noise is added to the profiles.

Sensitivity
The retrieval from six simulated profiles by both proposed configurations is shown in Fig. 4. The two configurations give equivalent results and successfully retrieve the simulated profile in cases a-d. Cases e and f return large costs, such that it is obvious they have failed. In case e, a different nonlinear correction was used in the simulation and causes underestimation of the state where the observed profile has maximal energy. The observation of a cloud in case f is reasonable within the PBL but fails above that. The large scattering within the cloud is outside of the range prescribed by the Sim. a priori and, though it obtains a decent fit to the visible region of the cloud, vertical correlations cause incorrect retrieval beneath it. Successfully fitting cloud and aerosol observations simultaneously requires a forward model and a priori designed for the several orders of magnitude spanned by the state vector. The lidar ratio profiles indicate that there is a decrease in the information content of the measurement above the PBL, where scattering (and therefore the magnitude of the return) is lower. The two configurations react differently to this. The lidar ratio configuration returns a smooth B profile that tends towards its a priori value above the PBL, as would be expected, while the extinction configuration gives a much noisier profile, indicating it is less influenced by the a priori.
A further six simulations containing small-scale fluctuations are presented in Fig. 5. The two configurations behave as before, with the lidar ratio mode returning a smoother profile but losing sensitivity above the PBL. The "layers" of cases g and h are correctly positioned by both modes, if slightly underestimated in magnitude. In case i, the layers are not resolved (see Fig. 6) as the features occupied only one retrieval bin. Doubling the resolution gives equivalent performance to cases g and h but with slightly increased noise and significantly increased processing time. Cases k and l are more difficult retrievals as they present lower signal-to-noise ratio (SNR). They are still consistent with the true profile but with greater errors. Both configurations give a respectable fit to the extinction profile, though they do increasingly underestimate the magnitude of peaks as they become narrower. This decreased sensitivity is clear within the averaging kernels (Fig. 7). Though the backscatter kernels are virtually delta functions in the PBL, the extinction and lidar ratio kernels have widths of 300-1000 m (the effective resolution of those products, which increases with height). The kernels also illustrate the loss of sensitivity in the lidar ratio configuration above the PBL, with the magnitude of both kernels decreasing significantly. In cases k and l, the sensitivity is also lower due to the reduced SNR. The extinction configuration maintains sensitivity throughout the profile, though its extinction kernels are skewed about their centre (which may derive from E (ra) measuring the integral of α, such that bins beneath a level contribute greater information content). Overall, the kernels indicate that the smoother profiles returned by the lidar ratio mode are due to a greater reliance on the a priori.

Retrieval error
The error covariance matrices for case a in both configurations are shown in Fig. 8. They confirm that the linear configuration makes the best use of the available information as that mode behaves identically within and above the PBL, whilst the logarithmic configuration reverts to the a priori covariance in the free troposphere. Where there is information, the  i, show adjacent bins are correlated to ∼ 60 %, whilst other nearby bins are slightly anticorrelated. The intercorrelation of the variables has also evolved, with points above a level being anticorrelated and those below positively correlated. It will need to be confirmed whether real data give similar results.
The diagonals of the covariances, plots a and f, can be used to approximate the error on the products. These give the bounds of Figs. 9 and 10. The lidar ratio a priori uncertainty is too small as it is not consistent with the simulated profile. Since that error is simply the a priori variance, this indicates that the a priori is overly constrictive. The extinction retrieval is better, though it underestimates the error in the PBL. Both fail to appreciate the error caused by assuming a non-paralysable dead-time correction when a different form (Donovan et al., 1993) was simulated.
The above-mentioned figures also compare the retrieval to the Ansmann method. For a fair comparison, the derivative is averaged over 300 m to be equivalent to the effective resolution of the retrieval. They are in good agreement in the PBL and the retrievals exhibit a lesser spread and error than the Ansmann solutions in the free troposphere. The Fernald-Klett method gives equivalent answers when given the correct lidar ratio.

Parameter error
In real retrievals, there will be some error in the model parameters b. This additional uncertainty that can be included in the retrieval by extending the measurement uncertainty to cover all sources of error, where K b = ∂F /∂b,b is the best estimate of the true parameters b, and the last term describes any inability of the forward model to describe the true state. Concentrating on only the parameter error for the moment, this can be implemented by replacing all occurrences of S with This significantly increases the computing cost of the retrieval, as S y must now be inverted in each iteration. A reasonable approximation is to only re-evaluate S y after the last iteration. The full calculation is considered in this section but will be relaxed in Sect. 4 where the quantity of data increases. The Ångström exponent can vary quite significantly but is commonly accepted to lie in the range 0.6-1.4, such that an error of 0.4 is reasonable (Klett, 1985). Radiosondes measure pressure and temperature at a given height with an accuracy of 0.5 hPa, 2 K, and 60 m, from which an error in number density and N of 0.5 % is expected (Kitchen, 1989). The height of the first observed bin, R 0 , can be easily estimated to within 10 m. The standard deviation of the data used to estimate E B can be easily derived.
The remaining parameters are estimated by some calibration procedure (e.g. Wandinger and Ansmann, 2002;Povey et al., 2012). For the purpose of demonstration, Fig. 11 shows the impact of each parameter on the total variance, assuming errors in C of 10 % and τ d of 1 ns; these dominate the total. For the impact of uncertainty in C to be of a similar order to the measurement error, it must be known to within 2 %an unrealistic expectation. However, these errors are unlikely to be correlated such that this term simply increases the total error. The dead time is more troublesome for elastic measurements as it can introduce significant correlations within S y . For the prototype system simulated with τ d = 50 ns and a maximum count rate of 17 MHz, an error greater than 0.1 ns in its estimation significantly reduces the information content available and prevents the retrieval from converging. That is clearly an unrealistic expectation but is a fair representation of the impact that dead time has on the observations for this system. Most laboratory-standard systems will have much smaller dead times, which have a greater tolerance of around 1 ns.
The laser energy E L behaves similarly to C but is considered separately as it can change significantly with time, whilst the calibration function should be fairly consistent. The laser energy may not always be accurately measured and so is retrieved as part of the state vector. As the a priori tends towards zero with height and the molecular component of the scattering is input into the retrieval as a parameter, the most favourable means for the algorithm to fit the signal at the top  of a profile is to vary E L , effectively calibrating the signal against the Rayleigh scattering (as frequently done in Raman analyses). The process could be made more explicit by setting the a priori equal to zero with a small uncertainty at the top of the profile, but this has not been found necessary.

Further errors
The settings of the retrieval that have no bearing on the forward model should not affectx. The initial value of i alters the number of iterations required to converge as it drives the size of each step in state space. A value of 10 5 appears to be optimal in most cases andx appears to be independent of that choice provided it is not too large or small. Similarly, convergence thresholds of 10 −4 on change in cost or step and 10 −1 on error were selected as the highest order for whichx is not affected by the choice. The minimum retrieved height does have a small effect on the retrieval in its first few bins, so r 0 = 100 m was chosen to concentrate these effects within a region where parameter errors will be large regardless. Forward model error is defined in Sect. 3.2.3 of Rodgers (2000) as where G y = ∂x/∂y, the sensitivity of the retrieved state to the measurement, and f is the exact, true profile including any processes the forward model F may not describe. This systematic error is generally difficult to estimate as, if f were known, it would most likely be used as the forward model instead.

Atmos
There are some processes that are clearly not included within the current forward model. Multiple scattering has been neglected as it is mostly important for lidars with a wide footprint, such as space-based system, or for observations within clouds, where this algorithm is already known to perform poorly for other reasons. Though appropriate numerical models of multiple scattering exist (Eloranta, 1998), this is left as an area for future work if retrievals within clouds are desired.
There is a small difference between the bin-averaged backscatter that is sampled and the true backscatter defined by Mie theory, for which the error can be evaluated with Eq. (18). It is greatest in the entrainment layer (or at any other sharp gradient), being at most 1 % of the total error.
The models of the calibration function and detector nonlinearity are idealised versions of the truth. A rough estimate of these contributions can be produced by considering alternative models, such as that of Donovan et al. (1993). For case e, these are over 100 times larger than other errors. This is a circumstance that could benefit from the simultaneous processing of photon-counting and analogue data. The other cases are negligibly affected by the choice of nonlinear correction.
None of these errors describe the discrepancies shown in the free troposphere in Fig. 9 as that is dominated by the a priori uncertainty. In regions where the data are the dominant contribution to the retrieval (i.e. where the area of the averaging kernel is near unity), increasing the a priori variance does not affect the retrieved profiles. Where it is important, the error estimate should clearly be greater to better represent the uncertainty. Hence, the a priori uncertainty in B will be increased to 40 sr. This is effectively a uniform distribution in Fig. 2.

EARLINET intercomparison
The European Aerosol Research Lidar Network (EAR-LINET) produced simulated data with which to compare the performance of the various algorithms used by its members (Pappalardo et al., 2004). These data have been processed with the proposed algorithm, shown in Fig. 12. The true profiles of backscatter, extinction, and the lidar ratio are shown over the top row in black, this retrieval in red, and the local formulation of the Ansmann method is shown with points. The bottom row plots the difference between the retrieved and synthetic profiles, normalised by the estimated error.
The extinction profile is consistent with that simulated throughout the retrieved range and performs equivalently to the EARLINET algorithms (Figs. 4 and 6 of Pappalardo et al., 2004). This retrieval fits the magnitude and shape of all three peaks more accurately than some of the EARLINET algorithms and does not indicate an erroneous peak above 5 km (though the magnitude of noise is similar).
The backscatter product is less satisfactory, with the retrieval overestimating the magnitude below 3 km. This is a matter of calibration. The EARLINET algorithms estimated the calibration constant for the elastic channel using a given value of the aerosol backscatter between 8 and 10 km (described as the Stage II comparison). As the extinction is derived from the gradient of the logarithm of the Raman channel, its calibration constant is not necessary. These are not ideal circumstances for the retrieval, as estimates of both constants are required. These were produced by additionally assuming χ = 0.14, the most commonly observed value at the Chilbolton AERONET station (Woodhouse andAgnew, 2006-2011). The uncertainty of those estimates exceeds 60 %. The retrieval of E L found a suitable value of the Raman calibration constant (as the retrieval returns virtually identical results when both calibrations are multiplied by a random constant), but it has no mechanism to alter the ratio between the two calibrations resulting in the overestimation of the backscatter. The purple line of Fig. 12 shows a retrieval where the elastic calibration has been increased by 10 %, giving a more accurate retrieval. This ratio should not change significantly over time and so would in practice be estimated from infrequent observation where χ (p) is known.
Two additional profiles are plotted in Fig. 12, where the a priori values (both x a and S a ) have been increased and decreased by 20 % (grey and blue, respectively). The solutions cannot be visually distinguished from the original as the changes are less than 10 −7 Mm −1 sr −1 for backscatter and 10 −5 Mm −1 for extinction, demonstrating that the proposed a priori is (as intended) a weak constraint on the solution. The differences can be resolved in the lower row of plots, showing the changes correspond to < 10 % of the total error. Figure 12 also includes an example of the retrieval without the requirement that all state vector elements be positive, shown in green. For this profile, the information content is sufficiently large that there is no difficulty in converging Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |  Fig. 12. Performance of the retrieval (extinction mode) with data simulated by EARLINET. Top -The simulated backscatter, extinction, and lidar ratio (black) compared to that from the Ansmann method (points) and retrieval (red). Also shown are the retrieval where the a priori is increased or decreased by 20 % (grey and blue), where the retrieval has been allowed to take negative values (green), and where the magnitude of the calibration function has been manually selected to improve the comparison. Bottom -The difference between those retrievals and the simulated profile, normalised by the estimated standard deviation (Ansmann not shown).
47 Fig. 12. Performance of the retrieval (extinction mode) with data simulated by EARLINET. (Top) The simulated backscatter, extinction, and lidar ratio (black) compared to that from the Ansmann method (points) and retrieval (red). Also shown are the retrievals where the a priori is increased or decreased by 20 % (grey and blue), where the retrieval has been allowed to take negative values (green), and where the magnitude of the calibration function has been manually selected to improve the comparison. (Bottom) The difference between those retrievals and the simulated profile, normalised by the estimated standard deviation (Ansmann not shown).
(though this has been an issue elsewhere). Retrieved values increase slightly in magnitude throughout the profile, increasing its deviation from the simulated profile. This compensates for the negative values in the profile, such that the integral of the extinction profile differs by only 1 % from its value with the positive lower limit. The limit will therefore be retained.

Individual profiles
The retrieval is now applied to observations by the Chilbolton Ultraviolet Raman lidar (CUV; Agnew, 2003;Agnew andWrench, 2006-2010), which is stationed at the Natural Environment Research Council (NERC) Chilbolton Facility for Atmospheric and Radio Research (CFARR; 51.1445 • N, 1.4270 • W; 84 m a.s.l.; STFC, 2006STFC, -2011. It uses a 355 nm Nd : YAG laser at 350 mJ and 50 Hz for water vapour profiling through the daytime boundary layer on a case-study basis, implementing both photon-counting and analogue data collection. Its observations can be directly compared to those of a Leosphere EZ lidar operated continuously at the same site, which provides volume depolarisation ratio profiles instead of Raman observations. Radiosonde launches are available twice daily from Larkhill, 30 km northwest (UK Meteorological Office, 2006-2011. Six profiles were selected from March 2010 for which the instrument's calibration has been thoroughly investigated using the techniques of Povey et al. (2012). Figure 13 compares the retrieved profiles to those given by the Fernald-Klett and Ansmann methods. A clear atmosphere is assumed between 4 and 5 km using a constant B to give an optical thickness consistent with sun photometer observations and the derivative is evaluated over 150 m. In the PBL, the retrieved backscatter is very similar to that given by the Ansmann ratio and an independent measurement by the EZ lidar. As the SNR decreases, the retrieval tends towards the Fernald-Klett solution. This is a proper response for the retrieval, giving answers similar to existing methods but tending from a two-to one-channel retrieval as the available information decreases. This is also expressed in the backscatter averaging kernels, which widen from 30 to 100 m.  The retrieved extinction is consistent with the Ansmann solution but gives a much smoother solution. The averaging kernels confirm that there is little information available in the free troposphere but also show that the resolution in the PBL is better than that observed in simulations: 100 m. A tendency to find α (p) = 0 at the top of the PBL in cases 1-4 is due to the number density profile. The radiosonde that morning recorded a step decrease in temperature at the top of the PBL, but as that is a low-resolution measurement, linear interpolation overestimates N there. A standard atmosphere does no better. Due to factors such as these, the inclusion of parameter errors significantly reduces the information content in the free troposphere. This can occasionally produce unconstrained solutions (not shown) due to the relatively weak -Larger background levels (being daytime observations) produce large correlations at the top of the profile. The two modes respond differently to this, with the logarithmic configuration reverting to the a priori covariance at the top of the PBL (as observed in simulations) and the linear configuration tending towards complete correlation (representing a more systematic error).

Discussion
-The intercorrelation of extinction and backscatter is different. Bins above a point are still negatively correlated, but those below are more weakly correlated. The exact reason for this is not clear.

Extended periods
Eleven hours of photon counting observations were processed from 2 March 2010. Figure 15 plots retrievals with an error less than 20 % or 30 sr (β (p) and B, respectively). Row a shows the application of the Ansmann method, where the data were averaged over 30 m to give a similar resolution to the retrieval. Row b is the linear retrieval and c the logarithmic. The three backscatter fields are qualitatively similar before 14:00 GMT and after 18:00 GMT. Between these times, the measurement of laser energy has diverged increasingly from reality. As the retrieval has no knowledge of that, it retrieves smaller backscatter to compensate. The Ansmann method is not affected as it considers a ratio of channels. The difference between the Ansmann solution and the retrieval is effectively a constant factor of the failure in the calibration, with the results otherwise being consistent. For example, both methods observe larger backscatter in updrafts than downdrafts (where vertical wind was observed by a Doppler lidar). The retrievals are consistent in their estimates of the lidar ratio and are no worse than the Ansmann method, which is greatly affected by overlap when estimating extinction. The aerosol layer near 1 km at 10:00 GMT gives a lidar ratio of around 30 sr. This is a residual layer lying above a developing mixed layer where lidar ratios are larger (around 50 sr). The low lidar ratio indicates large, likely spherical, particles which are reasonable for an aged residual layer. The results are better in the evening, observing a peak B of 70 sr over a background of 50 sr, indicating the appearance of smaller particles. By this time, convective mixing has collapsed into a persistent updraft, so the increase in depolarisation ratio could indicate that newer, non-spherical particles are being lofted from the surface or are advected over the site. Advected aerosol is more likely considering the brevity of the peak. Figure 16 compares the retrieved χ ∞ during that day to AERONET measurements (Woodhouse andAgnew, 2006-2011). Their agreement is reasonable if not impressive. The retrieval tends to return larger χ than observed by AERONET, though it also contains substantial variability that the latter does not. This is likely due to the inaccurate measurement of laser energy, though this is under investigation. A calibration performed at 10:00 GMT was used throughout this day and that is the only AERONET measurement that was in any way input into the retrieval -the remainder are independent. Regardless, the retrieved values are equivalent to those given by the Ansmann algorithm, indicating that the retrieval is correct for the parameters it has been given. It is the calibration of the system, not the method of retrieval, producing the poor comparison.

Eyjafjallajökull ash
The eruptions of the Eyjafjallajökull volcano in southern Iceland during April and May of 2010 produced the single most significant volcanic ash event over northern Europe in the age of aviation. The closure of airspace cancelled around 100 000 flights, inconveniencing millions of travellers across the globe and resulting in massive losses for airlines and related industries. Owing to the density of personnel and instrumentation within the reach of this plume, it has become one of the most studied atmospheric events in history. The introduction of Johnson et al. (2012) provides a reasonable overview of the literature published to date and more will certainly be published over the years to come.
The CUV was operated, in addition to routine measurements, on 19 April to observe ash within the boundary layer. The data suffer similar difficulties to those already discussed, such as the impact of the calibration function being clear below 600 m. With these limitations in mind, the optimal estimation retrieval (in extinction mode) was applied to these observations, shown for analogue data in Figs. 17 and 18. Plot a presents the volume depolarisation ratio observed by the EZ lidar, while plots c and d show the retrieved backscatter and lidar ratio with errors outlined in Fig. 18. Ash particles have a large depolarisation due to their asphericity, and these measurements indicate the presence of a 400 m thick ash layer within the PBL, highlighted in Figs. 17c and d by plotting the 0.07 depolarisation isoline. It exhibits a low backscatter (< 10 Mm −1 sr −1 ) and lidar ratio (20-35 sr) compared to the remainder of the scene. These are well outside the range of 50-82 sr reported in the literature for similar ash in the free troposphere (Ansmann et al., 2010;Marenco and Hogan, 2011;Hervo et al., 2012), though slightly larger than found if the data is analysed with the Ansmann algorithm (10-25 sr). The retrieval indicates that the properties of the ash have changed significantly after 12-24 h within the PBL. The decreased B implies a growth or shape change of the particles, possibly due to sulfate coating.
A mixed layer forms beneath the aerosol (see the vertical velocity in plot b). There, B = 50-80 sr with minimal depolarisation, which is broadly consistent with urban Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | aerosols (Müller et al., 2007). Backscatter is fairly homogeneous throughout the layer except during a period of updrafts around 13:00 GMT when β (p) decreases from 10 to 6 Mm −1 sr −1 . The absence of a similar change elsewhere in the PBL gives some confidence that this is a real variation rather than a calibration artefact.
A more weakly depolarising aerosol resides in a poorly mixed residual layer above the ash layer. It persists until 14:00 GMT, when they mix. Backscatter and lidar ratios are large at the top of this layer and increase with height. This could be simple stratification within a poorly mixed layer or smaller particles may have concentrated at the top of the layer whilst larger particles have begun to settle.
Finally, a thin layer of aerosol is present above the PBL (labelled in plot c). The EZ lidar did not resolve this, so no measure of the depolarisation is available. Expressing lidar ratios of 40-60 sr with low backscatter, the layer is consistent with aerosol typically observed at CFARR and there is no reason to label it as ash. Figure 19 presents the distribution of retrieved extinction and backscatter for all points for which the error in the volume depolarisation ratio is < 100 %. Lines of constant B are added for reference. Points likely to contain ash are shown in the left plot by filtering for depolarisations greater than 0.07. The residual layer appears as a concentration of points around β (p) = 10 Mm −1 sr −1 and B 40 sr. The mixed layer appears in the right plot as a more continuous distribution between B = 40 and 60 sr. The failure of the retrieval near the surface is evident in a vertical line of points at β (p) = 8 Mm −1 sr −1 . In the free troposphere, β (p) < 1 Mm −1 sr −1 , where poorly constrained retrievals produce a broad distribution in both plots. There are very few observations of the thin ash layer, but their presence is evident in observations near B 20 sr in the left plot not expressed on the right.

Conclusions
An optimal estimation retrieval scheme for aerosol scattering properties from Raman lidar observations was proposed, using the lidar equations as a simple forward model. The a priori state and covariance matrix were based on the properties of aerosol outlined in the OPAC model to provide a weak constraint on the magnitude of scattering whilst assuming them to be vertically correlated over a scale height of 100 m. This is smaller than observed by balloon-borne measurements but ensures that the PBL and free troposphere are not coupled.
The state of the atmosphere can be described at each height by the aerosol backscatter and either the extinction or lidar ratio. These possibilities were assessed by considering their ability to process simulated data. The lidar ratio configuration was found to lose sensitivity in the free troposphere, relying excessively upon its a priori assumptions, as shown by the disappearance of the averaging kernel. If extinction is retrieved instead, it and backscatter should be retrieved linearly with a correlation assumed between them (95 % here, though more investigation of this value is necessary). This configuration maintains sensitivity throughout the profile. This choice is likely influenced by the desire to use a single a priori appropriate for all observations. In future, if aerosol type information was used to better estimate the lidar ratio before analysis, the lidar ratio configuration may perform more favourably.
In the analysis of simulated and real data, the proposed retrieval is consistent with existing analyses. Backscatter was always retrieved at the finest resolution allowed (mostly 33 m, but this remains true at the instrumental limit < 10 m) and with an uncertainty between 2 % in the most ideal circumstances and 20 % in the least. Extinction and the lidar ratio are less well constrained, expressing resolutions of 300-500 m in simulations and 0.1-1 km with real data. Importantly, these are different from the scale of vertical correlations assumed a priori and increase as SNR decreases. The retrieval has selected the most suitable resolution independently, unlike the smoothing filters used in most studies. The uncertainty in extinction retrieved from real data is relatively large (> 15 %) but that is in part due to the short timescales evaluated (1 min in Figs. 15 and 17). The integration time can be increased to reduce these errors to any desired level (at least until atmospheric variability begins to dominate). Regardless, errors are comparable to, if not smaller than, fair estimates of the error resulting from the standard Raman lidar technique of Ansmann et al. (1992) applied to the same data.
The magnitude of the uncertainty was shown to be dominated by the calibration of the instrument -primarily the offset of its vertical axis, the nonlinear response of its detectors, and the calibration function. Ideally, all of these would be determined with dedicated laboratory measurements, especially the offset and nonlinearity, which should change very little over time.
Optimal estimation is a poor technique for exploring unknown and unusual circumstances (such as dense smoke or desert dust plumes) as it can only retrieve states that are consistent with the a priori. (This is generally taken to mean that the retrieved state must be within one standard deviation of the a priori, but it is more accurately considered a balance of the a priori and measurement uncertainties. A state far from the a priori can be returned if the uncertainty in its measurement is sufficiently small.) However, the retrieval will return a cost that indicates when a poor fit of the a priori is encountered. If observations are required there, the magnitude constraint can be altered or practically removed (by setting S a to a very large value), analogous to satellite aerosol retrievals using several aerosol models from which only the result with the lowest cost is reported. It should be emphasised that the proposed a priori is a weak, order-of-magnitude constraint and there is no reason to suspect the retrieval would perform poorly except in the presence of unusually bright or dense aerosols.
The retrieval was applied to several hours of observation on 19 April 2010 of ash from the Eyjafjallajökull eruption. A depolarising ash layer was observed with a lidar ratio of 20-30 sr, much lower than observed in the free troposphere by previous studies and potentially indicating a growth of the particles after 12-24 h within the planetary boundary layer. More dispersed ash within a residual layer exhibited a backscatter of 10 Mm −1 sr −1 and lidar ratio of 40 sr.