Development and Application of a Backscatter Lidar Forward Operator for Quantitative Validation of Aerosol Dispersion Models and Future Data Assimilation

A new backscatter-lidar forward operator was developed which is based on the distinct calculation of the aerosols’ backscatter and extinction properties. The forward operator was adapted to the COSMO-ART ash dispersion simulation of the Eyjafjallajökull eruption in 2010. While the particle number concentration was provided as model output variable, the scattering properties of each individual particle type had to be determined by extensive scattering calculations. Sensitivity studies were 5 performed to estimate the uncertainties related to the assumed particle properties. Therefore, scattering calculations for several types of non-spherical particles required the usage of t-matrix routines. Due to the distinct calculation of the backscatter and extinction properties of the models’ volcanic ash size classes, the sensitivity studies could be resolved to each size class individually which is not the case for forward models based on a fixed lidar ratio. Finally, the forward modeled lidar profiles have been compared to ACL measurements both qualitatively and quantitatively while the attenuated backscatter coefficient 10 was chosen as common physical quantity. As the ACL measurements were not calibrated automatically, their calibration had to be performed using CALIPSOs/CALIOP measurements. A slight overestimation of the model predicted volcanic ash number density was observed. By manually reducing the model predicted ash number density, the effect of simple data assimilation methods could be demonstrated. Major issues for future data assimilation of ACL data have been identified. The introduced forward operator offers the flexibility to be adapted to a multitude of model systems and measurement set-ups. 15


Introduction
In Spring 2010, the volcano Icelandic volcano Eyjafjallajökull erupted several times. The emitted ash was found to be harmful for aircraft and due to uncertain information about spatial distribution and concentration of volcanic ash, the European air space was closed for several days (Sandrini et al., 2014). The high economic costs and impact on public transport lead to efforts of DWD (Deutscher Wetterdienst) to improve at monitoring and predicting aerosol cloud movements in the atmosphere. 20 Therefore, DWD decided to start a dedicated project on backscatter lidar forward operators for validating aerosol dispersion models using available remote sensing measurement data.  Chen et al. (2014)). Using these model systems, scientists have analyzed the aerosol influence on, for example, precipitation (Rieger et al., 2014), temperature (Bangert et al., 2012), radiative fluxes (Vogel et al., 2009;Chaboureau et al., 2011).
Lidar (light detection and ranging) is capable of providing information on atmospheric particles with high temporal and 10 spatial resolution. The most basic lidar type is the backscatter lidar which measures the backscattered signal intensity of a volume at a certain range. Comparing the data of such a backscatter lidar that is operated in the UV with simulations of an atmospheric chemistry model, allows, e.g., for the characterization of transport and optical properties of aerosol particles near sources Álvaro M. Valdebenito B et al., 2011). Using ground-based DIAL (differential absorption lidar, Dagan (2008); Späth et al. (2016)) water-vapor can be measured, which can even be combined with backscatter measurements 15 to derive more details of aerosol particle properties (Wulfmeyer and Feingold, 2000). Lidar techniques based on the vibrational and rotational Raman effect, like RRL (rotational Raman lidar) allow for the measurement of trace gas profiles (Whiteman et al., 1992;Turner et al., 2002;Wulfmeyer et al., 2010), as well as profiles of atmospheric temperature, particle backscatter cross section, particle extinction cross section, and particle depolarization properties (Behrendt et al., 2002;Hammann et al., 2015;Radlach et al., 2008). HSRL (high spectral resolution lidar) systems furthermore allow for cloud and particle characterization 20 (Shipley et al., 1983). Multiwavelength lidar systems offer the potential to retrieve the optical, microphysical and chemical properties of aerosols (Mamouri et al., 2012) but these systems are rare and the inversion algorithms are very complex. Profiles of the radial wind speed can be obtained by Doppler-Lidar systems, e.g., Banta et al. (2012).
While the number of sophisticated lidar instruments that provide thermodynamic data  is still low, there are already automated aerosol lidar networks in operation in Europe and Asia Sugimoto et al., 25 2008). For example the German weather service (Deutscher Wetterdienst, DWD) maintained already in 2010 a network of about 36 automated ceilometer lidar (ACL) systems (Flentje et al., 2010a). The data of such a network offers 3D particle information with a high temporal, high vertical and a moderate horizontal resolution. The ACL systems have been used to detect cloud and boundary layer heights (Emeis et al., 2009) but the received signal delivers also information about aerosols.
It is therefore worthwhile to use the ACL network measurements for the validation of particle transport model simulations. 30 Unfortunately, it is not possible to obtain the particle number concentration from an elastic backscatter signal alone without ancillary information and assumptions which are partly critical. The preferable way is thus to use the detailed atmospheric description of the model to simulate lidar profiles for a model-given atmospheric state. Then, the quantitative comparison can be performed between the simulation and observation more accurately. Such a lidar simulator is called lidar forward operator. The lidar forward model of SIČ (2014) has the capability to calculate the atmosphere optical depth (AOD) and the particle extinction coefficient from simulations of the CTM MOCAGE (Chemical-Transport Model, Modèle de Chimie Atmosphérique à Grande Echelle). But similar to the ECMWF forward operator, assumptions on the lidar ratio S lidar are mandatory to calculate the particle backscatter coefficient. At MetOffice, UK, a backscatter lidar forward model is being investigated (Charlton-Perez et al., 2013) which considers clouds, one type of aerosol, and rain simulated by the Unified Model (UM). Another lidar for-10 ward model is being developed for the EURAD-IM (European Air Pollution and Dispersion -Inverse Model, Lange and Elbern (2014)), but no published results have been found in literature. Consequently, all existing lidar forward operators which we found in literature calculate only the extinction coefficient. The backscatter coefficient, however, is not calculated explicitly for given atmospheric scatterer species. On the one hand, this method benefits from the fact that the extinction coefficient is less sensitive to the particle properties than the backscatter coefficient. On the other hand, the precision of this method is limited 15 to the correctness of assumed lidar ratio. The method becomes furthermore unusable once there is a mixture of scatterers for which the effective lidar ratio is not known. To become independent on the assumption of a lidar ratio, we designed a forward operator which is based on the independent calculation of both the extinction and the backscatter coefficients.
Our forward operator can be adapted to any particle-representing atmospheric model system and backscatter lidar system even using multiple wavelengths. The forward operator is optimized to calculate both the attenuated backscatter coefficient and 20 the lidar ratio from model output data with a minimum set of external information. We call this forward model the "backscatter lidar forward operator" (BaLiFOp). This forward operator is planned to be used for particle data assimilation.
Within its priority project KENDA (Kilometer-Scale Ensemble Data Assimilation) the COSMO Consortium has developed an Ensemble Kalman Filter for data assimilation on the convective scale. It is scheduled for introduction into operational use by MeteoSwiss and DWD in 2016. An advantage of the ensemble data assimilation system is that the assimilation can be 25 carried out based on the pure forward operator, and that it is not necessary to calculate derivatives of the forward operator or the adjoint tangential model for carrying out data assimilation. Also, it naturally introduces model increments for all variables where some dynamic covariance is observed from the underlying ensemble model runs. DWD aims to test the assimilation of ACL data into the COSMO-ART model based on the BaLiFOp Lidar forward operator. Here, the current work is a necessary and important basis for further work to improve numerical weather prediction and particle forecasting. 30 In the following we explain the lidar principles and the theoretical background for the backscatter lidar forward operator (Sect. 2). This is followed by an introduction to the case study in Sect. 3. Sensitivity studies of the particles' scattering properties are presented in Sect. 4. Results of the forward operator and a comparison to ACL measurement data are shown in Sect. 5.
Finally, we discuss both the benefits and the requirements of ACL data assimilation systems (Sect. 6).

The Lidar Equation
The lidar principle is based on the emission of UV, visible or IR laser pulses into the atmosphere and the measurement and analysis of the backscatter signals. The received photon number per pulse N rec,λ (z) from range z is described by the following equation for elastic backscatter lidars which detect the backscatter signal at the emitted wavelength Instrument-dependent variables of the lidar equation are the wavelength λ, the laser emitted photon number per pulse N tr,λ , the temporal length of a laser pulse τ , the speed of light c, the efficiency of the receiving system and detectors η λ , the overlap function O(z), and the net area of the receiving telescope A tel . The received signal intensity can be given either as power or in photon counts. Here, we use photon counts per laser pulse unless otherwise noted.

10
The range resolution is related to the temporal resolution of the data acquisition system by τ c 2 = ∆z. Typical ∆z values for ACL systems are a few meters. The overlap function O(z) is zero (no overlap) near ground and becomes 1 (full overlap) above a certain height which is typically 200 m to 1500 m above ground for ACL systems Flentje et al., 2010a). The missing overlap limits the capability to measure in near range but has no effect on heights where full overlap has accomplished. Heights where 0 < O(z) < 1 can be overlap-corrected if the device specific overlap function is known. 15 Processes in the atmosphere are described by the backscatter coefficient β λ (z) and the extinction coefficient α λ (z). The backscatter coefficient β λ (z) describes the scattering strength into the direction of the receiving telescope and depends on the wavelength, the type, shape and size of scatterers, and their respective number-concentrations; β λ (z) is given in units of m −1 sr −1 . The extinction coefficient α λ (z) is a description for the light absorption and light scattering capabilities of objects in a volume; it is given in units of m −1 . 20 Elastic backscatter lidar systems do not allow for the separate measurement of β λ (z) and α λ (z) as two unknowns cannot be determined with one measured variable. For calibrated backscatter lidar systems, it is thus convenient to calculate the attenuated backscatter coefficient γ λ (z) from the measured profiles It is given in units of m −1 sr −1 . The attenuated backscatter coefficient is independent of all instrument-specific parameters 25 except the wavelength. Therefore, it is the best suitable physical quantity for a comparison between ACL measurement and aerosol model using a forward operator as long as no ACL measurements of extinction and backscatter cross section profiles are available for this purpose.

The Backscatter Lidar Forward Operator
According to Eq.
(2), the basic task of the forward operator is to calculate the extinction coefficient α λ (z) and backscatter coefficient β λ (z) based on the atmospheric state simulated by a model. Knowing the extinction coefficient of a vertical column allows for the solution of the two-way-transmission integral and, finally, to determine the attenuated backscatter coefficient γ λ (z). In the following, we describe our method of calculating the attenuated backscatter coefficient from a given set of 5 atmospheric scatterer data.

Scattering of Light by Arbitrary Objects
If there are q s different types of scatterers in an illuminated volume, the total extinction coefficient α λ (z) and the total backscat- where n i (R, z) is the number-size distribution of scatterer type i at range z given in units of m −3 , σ ext,i,λ is the extinction cross-section of scatterer type i given in units of m 2 with radius R, and dσ sca,i,λ dΩ π is the differential backscatter cross section given in units of m 2 sr −1 . 15 For isotropic scattering, the differential backscatter cross section is derived from the scattering cross-section σ sca,i,λ (R) via For non-isotropic scattering, a phase function φ i,λ (θ, R) is used to describe the relative scattering intensity into an angle θ.
In the case of monostatic ACL systems, we thus get θ = π, giving 20 In the following, we treat molecule scattering and particle scattering separately, as the respective calculations depend on different physical theories and algorithms.

Scattering by Molecules
Assuming that a model distinguishes atmospheric gases such as nitrogen, oxygen, argon, and water vapor, the molecule scattering calculation can be performed for each individual gas type and molecule size using the Rayleigh theory (Young, 1981). 25 This becomes relevant if inelastic scattering or molecular absorption is involved. For ACL systems provided that a wavelength is used, which is well outside of molecular absorption lines, the individual gas contribution to the signal does not need to be distinguished.
Consequently, the molecule extinction coefficient α mol,λ (z) and the molecule backscatter coefficient β mol,λ (z) can be calculated with The molecule number density N mol (z) is related to the ideal gas law where p is the model prediction of the atmospheric pressure given in Pascal (Pa), T is temperature given in Kelvin (K), and k 10 is the Boltzmann constant which has a value of about 1.381 × 10 −23 J K −1 .
To calculate the scattering cross-section σ sca,mol,λ and the scattering phase function φ a,λ (θ) of air, empirical equations are available. We used the formulas and look-up tables given by (Buchholtz, 1995). As the empirical equations are only given for wavelengths up to 1000 nm, we simply extrapolated the values to the ACL wavelength in the case study (1064 nm). This method allows us to calculate the scattering properties of air molecules σ sca,mol,λ and dσ sca,mol,i,λ dΩ π directly from model 15 output for a given ACL laser wavelength λ.

Scattering by Particles
The scattering characteristics of spheres with sizes not much smaller or larger than the wavelength are described by Mie's solution of the Maxwell equations (Mie, 1908;Wiscombe, 1980). Methods like the T-matrix (Mishchenko et al., 2002) or the discrete dipole approximation (DDA, Draine and Flatau (1994)) allow for scattering calculations of non-spherical objects; 20 again, with sizes not much smaller or larger than the wavelength. The T-matrix algorithm is a tool for computing scattering by single and compounded particles (Mishchenko et al., 2002). It is faster than DDA but limited to rotationally symmetric objects such as ellipsoids, cylinders or Chebyshev polynomials. DDA, however, can represent arbitrarily shaped objects at the cost of high computational efforts.
As a rough estimate from our studies, the computational costs increase by about one order of magnitude when using the 25 T-matrix instead of the Mie approach and by another two orders of magnitude when using the DDA instead of the T-matrix approach. Another increase in computational time is resulting from larger scatterers, i. e., an increase of the particle size causes an exponential increase of computing time. We therefore utilize Mie scattering algorithms to perform fast calculations although solid particles are in fact non-spherical. The effect of scattering by non-spherical particles is analyzed in a second step using the T-matrix approach for sensitivity studies. We consider these approaches as sufficient for the sensitivity studies performed Mie-scattering related computations were performed using the IDL procedure "mie_single", provided by the Department of Atmospheric, Oceanic and Planetary Physics (AOPP), University of Oxford. Input parameters of the procedure are the real part m and imaginary part m of the refractive index as well as the so-called size parameter X λ (R): where R is the radius of a single particle. The relevant output parameters are the extinction efficiency Q ext,p,λ (R), the scattering 5 efficiency Q sca,p,λ (R), and the backscatter efficiency Q bsc,p,λ (R). These optical efficiencies are defined as ratio between the optical cross section and the physical cross section: As a warning, we like to point out that the procedure changed its definition of the backscatter efficiency: The current (2012) release of mie_single returns the so-called radar backscatter efficiency which is 4π times the backscatter efficiency as we require it within the forward operator. Furthermore, the procedure expects the imaginary part of the refractive index given as 15 negative number. If positive imaginary part values are used, the procedure runs without showing an error but returns wrong results.

Discrete Particle Number Size Distributions
A major problem of discrete size distributions is the high sensitivity of the optical cross sections to the particle size: A slightly different particle radius may lead to quite a large change of the scattering properties. We present in the following a suggestion 20 to overcome this fundamental problem. Due to the fact that naturally occurring particle size distributions are not discrete, averaging the optical cross sections over certain size-intervals seems straightforward. We will show that this approach indeed reduces the problematic and unrealistic sensitivity significantly. If the model represents only one type of particle, i. e. with a constant refractive index but with discrete radii R d , we can define the effective extinction cross section and the effective backscatter cross sections with 25 σ ext,R d ,m,m ,λ = 1 where R da and R d b are size margins for each particle size class d. These integrals are then exchanged by sums in the numerical computation routines giving where n samples is the sampling number and the sampling range R d b − R da is broken down in g subsamples: This calculation of the effective values is performed for every discrete size class d and -if represented by the model -also 10 for every particle type k.
Consequently, the total particle extinction coefficient α par,λ (z) and the total particle backscatter coefficient β par,λ (z) are calculated from: Here, N d,k is the particle number per volume given by the model, σ ext,R d ,m k ,m k ,λ and σ bsc,R d ,m k ,m k ,λ are the effective optical cross sections of particle size class d and particle type class k with the respective real part m k and imaginary part m k of the refractive index.
This simple solution allows for calculating α par,λ (z) and β par,λ (z) by just solving a few multiplications and summations 20 resulting in a minimal demand of computing time.
The forward modeled total extinction coefficient and total backscatter coefficient are the sum of the molecule and the particle extinction and backscatter coefficients: equivalent to Eq. (3) and (4).

Two-Way Transmission
The two-way transmission T λ is calculated by while T λ depends on α λ between the instrument (z = 0) and z. Within the forward operator, the two-way transmission is discretized by: where the number of height levels between ground and the actual height is expressed by n z (z), the profile of the extinction coefficient is α λ (z) and the profile of layer thicknesses is ∆h(z) cos Θ with Θ as zenith angle. In the case of forward modeling, a vertically pointing measurement (Θ = 0), the vertical column of model grid boxes at a specific location is required. By changing the zenith angle Θ, the forward operator can also simulate scanning measurements.

Lidar Ratio
Even though the lidar ratio is not measured directly by current ACL systems, the skill of simulating the lidar ratio for given scatterer types and scatterer mixtures offers a great potential for sensitivity studies but also for applications to research lidar systems which are capable of measuring this parameter such as Raman lidar.
where α meas,par,λ (z) and β meas,par,λ (z) are the measured extinction and backscatter coefficients which represent the numbersize distribution of real scatterers. Assuming a monodisperse number-size distribution (N d,k is a factor), the effective lidar ratio of each size class can be calculated by: 20 where σ ext,d,kλ and σ ext,d,kλ are the extinction and backscatter cross section of a particle size class d and particle type k. This definition is limited to the scope of monodisperse size classes, i.e. the number-size distribution inside a size class is assumed to have a constant value. To calculate the total lidar ratio at distance z, the following equation could be used The 2010 Eyjafjallajökull eruption was extensively analyzed by scientists from many fields of research, resulting in a substantial knowledge base (see ACP special issue "Atmospheric implications of the volcanic eruptions of Eyjafjallajökull, Iceland 2010").
Ash layers were observed from a large set of measurement instruments, allowing for tracking the volcanic ash cloud movement over Europe (Gasteiger et al., 2011a;Zakšek et al., 2013;Mona et al., 2012;Dacre et al., 2013;Waquet et al., 2014). Using images from the geostationary instrument SEVIRI (Spinning Enhanced Visible and Infrared Imager) the spatial extend of the 10 ash cloud and its movement could be tracked and compared to the measurement of ground-based instruments (see Fig. 1). From the synergy of the two measurement systems, signals of the ACL systems could be related to clouds or volcanic ash layers, respectively.
A volcanic eruption case is also beneficial for our analysis due to the fact that the ash emission takes place at one geographic location within well-known time intervals. Other cases are also interesting but may lack a precise definition of the test scenario 15 or respective model prediction capabilities, such as the representation of background aerosol or uncertainties of the background aerosol number density.

The DWD ACL Network
ACL networks are the current data source of choice for an analysis of the vertical structure of aerosols. In 2010, the DWD ACL network consisted of 36 systems CHM15k ACL manufactured by Jenoptik now Lufft (see Fig. 2). A qualitative analysis of the 20 ash cloud over Germany with this data set was made by Flentje et al. (2010a).
We used the NetCDF files with ACL raw data where one file contains the 24 hour measurement of one ACL station. The received photon number per shot is calculated from where beta_raw is the signal-to-noise measurement product, stddev is noise, and base is a daylight correction provided by the 25 ACL software of this model. The calculation of the attenuated backscatter coefficient from ACL data follows Eq. (29) which The emitted photon number per shot N tr,λ was calculated using: where E pulse,λ is the laser pulse energy. E p,λ is the photon energy, calculated according to: with h as Planck's constant having a value of 6.626 07 × 10 −34 J s. The pulse energy of the diode-pumped laser is 8 µJ Unfortunately, the instruments provided no calibrated measurement data at that time. As the true system efficiency η λ and the 10 calibration coefficients are not known, we use the symbol η * as linear calibration factor. From a comparison with the calibrated attenuated backscatter measurements of CALIOP at λ = 1064 nm (Fig. 3), we determined a calibration factor of η * = 0.003: First, we used the CALIOP value of the 1064 nm calibrated attenuated backscatter coefficient at 50.15°lat / 4.81°lon in a height of 2 km ASL (about 5 × 10 −6 m −1 sr −1 ) to estimate the maximum attenuated backscatter coefficient inside the volcanic ash cloud. The second calibration criterion was the forward modeled attenuated backscatter coefficient outside the volcanic 15 ash layer in 3 km ASL which is dominated by molecule scattering and has a value of 1 × 10 −7 m −1 sr −1 , see Fig. 15. This pragmatic approach is of limited precision and only required for uncalibrated ACL measurements. As most ACL networks have been extended by automatic calibration methods after the Eyjafjallajökoll eruption, such a calibration will not be required in future forward operator studies.
As a last step, the high-resolution ACL data was gridded to the model's vertical resolution as well as over 15 min by temporal 20 averaging. This also improved the signal-to-noise ratio of the ACL data. Volcanic ash was represented by six discrete size classes with a spherical shape and aerodynamic diameters of 1 µm, 3 µm, were treated as equal, i. e., they had the same size, shape, and complex index of refraction (monodisperse classes). We therefore calculated effective optical cross sections as explained in Section 2.2.4. The lower and upper size margins R da and R d b were defined as the arithmetic average between two subsequent size classes. The lower margin of the smallest size class was half its nominal diameter. The upper margin of the largest size class was 1.5 times its nominal diameter. The resulting class ranges are shown in Fig. 4.

5
A list of model variables used for the forward operator is given in Table 1. While the p number density is taken directly from model output, the molecule number density of air was calculated from temperature and pressure according to Eq. (9).

Volcanic Ash Properties
A detailed analysis of the emitted ash was performed by Schumann et al. (2011)

15
Electron microscope images from the same study revealed that the volcanic ash particles were sharp-edged with a complex and asymmetric shape. The average asymmetry factor was 1.8 for small particles (<0.5 µm) and 2.0 of larger particles (Schumann et al., 2011). Electron microscope measurements of Rocha-Lima et al. (2014) showed that the asymmetry factor of the volcanic ash fine fraction has a value between 1.2 and 1.8.
The particle growth due to hygroscopic water coating was quantified to about 2 to 5 % at a relative humidity of 90 % 20 (Lathem et al., 2011). A growth of 5 % does not change the scattering properties significantly in relation to the size averaging we perform.
Nevertheless, each measurement contains an error and the ash properties may change during its travel through the atmosphere. It is therefore straightforward to analyze the maximum error due to variations of the volcanic ash properties, namely the particle size, refractive index, and shape. of complex shaped particles but the analysis is limited to size-parameters up to 20.8 (which results in a maximum particle radius of about 3.2 µm at a wavelength of 1064 nm). A study of Kemppinen et al. (2015) is focused on individual ellipsoids but we neither know the ellipsoidal properties nor does COSMO-ART predict changes of the volcanic ash classes. Guessing an ellipsoidal distribution for the model predictions thus may lead to less realistic scattering calculation results than assuming spherical scatterers. Consequently, there is no sufficient scattering description for Eyjafjallajökull ash predictions of COSMO-5 ART available. We therefore treat the volcanic ash as spherical objects with given optical properties (see 3.4) but analyze and discuss the effect of uncertain volcanic ash properties in the following.
It must be noted that these studies are required for nearly any aerosol as most naturally-occurrent aerosols are not perfectly spherical. As we will show below, already slightly aspheric ellipsoids may have very different scattering characteristics compared to perfect spheres.

25
These look-up tables were the base for the refractive index sensitivity study. Fig. 6 show σ ext and σ bsc as well as σ ext and σ bsc against the real and imaginary parts of the complex index of refraction. While the extinction cross section σ ext is more sensitive to the real part than to the imaginary part of the refractive index, the backscatter cross section σ bsc is strongly sensitive to both. These sensitivities are strongly reduced for the effective An overview of the refractive index sensitivity of the effective optical cross sections is given by Fig. 7 which shows the relative errors

Fig. 5 and
and σ bsc,err,p (m, m ) = σ bsc,p (m, m ) − σ bsc,p (m * , m * ) The relative error is the error of the optical cross sections if we assume that the refractive index we defined as reference (m * and m * ) to be true but the real particles have a refractive index of m and m . We can conclude from this analysis that the maximum relative error for the given range of refractive indices is less than 10 % for the extinction cross section but ranges up to 230% for the backscatter cross section.

T-matrix Particle Shape Sensitivity Study
10 T-matrix calculations within this study are based on the FORTRAN code for randomly oriented particles written and provided by Mishchenko and Travis (1998). A detailed description of the method can be found in the work of Mishchenko et al. (2002).
We modified the double-precision version of the T-matrix procedure to perform scattering calculations of multiple particle sizes automatically. In addition to that, the procedure was extended by the calculation of the backscatter cross section σ bsc according to Mishchenko et al. (2002), Eq. (9.10). Then, as a test, both mie_single and our modified T-matrix code were set 15 up to calculate the scattering properties of the same spheres and for the same wavelength. The results of both procedures were indeed identical.
A list of T-matrix options we used for the particle shape sensitivity study is shown in Table 2. The most important particle properties are defined by the variables NP and EPS. NP defines the particle type and has a value of -1 for spheres as well as for ellipsoids. A NP value of -2 is used for cylinders. The variable EPS is an expression for the objects' diameter to length ratio. 20 Consequently, an ellipsoid with EPS=1 is a sphere, prolate objects have an EPS<1 and oblate objects have an EPS>1.
In Figs 8, 9, and 10, the optical cross sections and the lidar ratio of spheres and several aspherical particles are plotted against the equal-volume radius. The aspherical scatterers are 6 ellipsoids with a diameter-to-length-ratio of 0.50, 0.67, 0.75, 1.25, 1.50 and 2.00 as well as 5 types of cylindric particles with a diameter-to-length-ratio of 0.50, 0.80, 1.00, 1.25 and 2.00.
As shown in these plots, the results for a highly asymmetric ellipsoid (EPS: 0.50) is only available up to an equal-volume 25 radius of 3.75 µm. For future research activities in this topic, the quadruple precision version of the t-matrix code should be preferably used.
We found no significant differences of the extinction cross section of spheres and these ellipsoids. By trend, cylindric shaped particles have a higher extinction cross section compared to ellipsoids and the spherical shape has the lowest extinction cross section values over the whole spectrum. Up to a volume-equivalent radius of 0.7 µm, the shape effect is not noticeable. to ellipsoids except for one type of ellipsoid (EPS = 1.25). For cylinders, such interference effects are not observable so the backscatter cross section of the considered cylinders increases monotonically with the equal-volume radius. As a result, the backscatter cross section of spheres is lower than of cylindric particles with the same size if their equal-volume radius is greater than 3.75 µm (for the given wavelength of λ = 1064 nm).
The particle shape effect on the lidar ratio is only weakly pronounced for small particle sizes (less than 0.75 µm), too. For 10 larger particles, the lidar ratio of spheres is generally lower than of the other considered shapes which is in agreement to the higher backscatter cross section observed before. For the fourth size class (equal-volume radii around 5 µm), the previously observed interference effects of the spheres' backscatter cross section leads to extreme values of the lidar ratio (exceeding a lidar ratio of 200 sr). By trend, however, the lidar ratio of spheres in this size class is not that different than the lidar ratio of other considered shapes. But for the other size classes, the lidar ratio of spheres is the lowest of all observations except for 15 cylinders which have a lidar ratio value at the same order of magnitude as spheres.
A summary of the particle shape sensitivity study is shown in Figs 11, 12, and 13, giving the relative differences of the effective optical cross sections and average lidar ratios for different particle shapes. The definition of the relative differences follows Eq. (32) and Eq. (33). Again, the reference is a spherical particle. Positive and negative relative differences indicate that the calculated values for spheres are underestimations and overestimations compared to respective non-spherical particles. 20 The effective extinction cross section of spheres is lower than the effective extinction cross section of other analyzed asymmetric particles. Regarding the effective backscatter cross section, however, we find relative differences of up to 300% and −80%. While the small aspherical particles have a lower effective backscatter cross section compared to spheres, the effective values of the fourth size class are higher for almost all considered aspherical particles. From this analysis, it can be concluded that due to the assumption of sphericity, the backscatter cross section of size classes 1, 2 and 3 are overestimated by about factor 25 1.5 to 5 while the backscatter cross section of the fourth size class is underestimated by factor 2. This allows for quantifying the over-and underestimation of the results for each size class individually which is not possible for forward operators where a fixed lidar ratio is assumed.
The relative difference of the average lidar ratio shows that assuming spheres for the scattering calculations results in much lower lidar ratio values of size classes 2 and 3 but in higher lidar ratio values than would be observed for a mixture of particle 30 shapes. For the fourth size class, the opposite has been observed: assuming spheres leads to an higher values of the lidar ratio of large particles than for non spherical particles of the same size. The lidar ratio of the first class, however, is nearby independent on the particle shape.
Of course, the total effect of different particle shapes for a real aerosol mixture has to be estimated with more extensive scattering calculations. For the current state of the forward operator research, the spherical shape has to be used as reference 35 15 Atmos. Meas. Tech. Discuss., doi:10.5194/amt-2017-142, 2017 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 17 May 2017 c Author(s) 2017. CC-BY 3.0 License. even if the real volcanic ash particles are known to be fractal and complex shaped. One of the reasons for doing so is the fact that there is currently no appropriate shape representation scheme for volcanic ash available. For example from the findings of Rocha-Lima et al. (2014), the average aspect ratio of volcanic ash is known but not a representative particle shape. As we have shown, the backscatter cross section of cylinders also differs tremendously from ellipsoids and from spheres etc. so there is currently a lot of indication that using a given aspherical shape as reference would yield in even higher errors than assuming 5 a spherical shape.
The above results, however, give us valuable insight into the uncertainties of using Mie calculations for non-spherical particles. These findings are important for interpreting the results of the forward operator.

10
A list of the effective extinction cross section and the effective backscatter cross section we determined for atmospheric gas molecules and for the six volcanic ash size classes is shown in Table 3.

Output Variables of the Forward Operator
Using the forward operator allows for plotting each variable of the lidar simulation for analytic purposes (see Fig. 14). These plots of forward-operator output variables are representing the major characteristics of the variables: strong extinction and 15 strong backscattering are usually related. Time and height intervals at which only molecules exist, lead to low values of the extinction coefficient and backscatter coefficient. Due to the decrease of the atmospheric gas number density with height, both extinction and backscatter coefficient decrease with height in an aerosol-free atmosphere. The two-way transmission decreases with height (see Eq. (23)).
In addition, the volcanic ash contribution to the total signal and the total mass density was analyzed for each size class of 20 COSMO-ART within two case studies. The cases were selected to cover extreme situations: Case 1 is the model output from a coordinate inside the volcanic ash layer (Table 4). Case 2 is for a coordinate where the majority of particles are due to size class 4 and 6 (see Table 5).
Regarding case 1, the total backscatter coefficient is dominated by ash size classes 1, 2 and 3 while the signal contribution of classes 4 to 6 is less than 5 % in total. The mass contribution is dominated by classes 3 and 4 while classes 2, 5 and 6 are 25 contributing by 10 % each to the total mass density and class 1 is very low for the total mass density. Regarding case 2, the total backscatter coefficient depends by about 68 % from class 4 and by 30 % from class 6. The mass contribution in case 2 is also dominated by the classes 4 and 6 but, in contrast to the backscatter coefficient, class 6 has a higher contribution to the total mass density than class 4.
General conclusions from this analysis about the relationship between backscattering and mass depending on particle size 30 and wavelength require further investigation. For our application of the forward operator in this study, however, we can conclude 16 Atmos. Meas. Tech. Discuss., doi:10.5194/amt-2017-142, 2017 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 17 May 2017 c Author(s) 2017. CC-BY 3.0 License.
the following: First, the total signal inside the volcanic ash layer (case 1) is predominately dependent on classes 1, 2 and 3 whose backscatter cross sections are also overestimated by the forward operator due to the assumption of sphericity (see Fig.   12). The real values of the attenuated backscatter coefficient may be by factor 2-3 higher. Second, the larger particles of classes 4, 5 and 6 carry a large portion of the mass but contribute only weakly to the total signal. This may be an important information for the selection of future ACL networks as even the systems operating at 1064 nm have a reduced sensitivity for particles 5 within these size classes.

Qualitative Comparison
A qualitative comparison allows for the identification of common and different structures between the measured and simulated lidar profiles. Different ash layer structures can hint, e.g., to errors in the model dynamics, in the source description, or in the sedimentation parametrization. If ash structures are found only in the measured profiles, either the model prediction is wrong 10 or it misses an important aerosol type which is not present in the model. If structures are visible in the forward modeled profiles but missing in the measured profiles, either the ACL signal is too weak because of high extinction in lower heights or the model performed a wrong ash prediction. But it is also possible that the model overestimates the ash concentration so the structures are below the detection limit of the ACL measurements.
For a qualitative comparison between measurement and simulation, we chose the ACL station Deuselbach in West Germany 15 and a time interval from 16 to 17 April 2010 as an example. Here, the ash layer was clearly visible in the measured profiles without being affected by low-level or high-level clouds. We calculated the attenuated backscatter coefficient from the ACL measurement according to Eq. (29), extracted the common time and height intervals and re-sampled the ACL data to the model resolution.
A comparison of ACL measurement and COSMO-ART simulation with an applied forward operator is shown in Fig. 15. Due 20 to the inevitable instrumental noise, the automatic calibration of ACL system and subsequent background subtraction, some data points become negative which is just a statistical effect but causes missing data in the log-scale plots. Volcanic ash plumes are clearly visible on both plots. Looking at the forward operator result, the ash layer begins to cross the ACL station between 06:00 UTC and 12:00 UTC at 16 April 2010. The layer height decreases with time and partially entrains into the planetary boundary layer where it persists even at the end of 17 April 2010. As both model and forward operator only represent volcanic 25 ash and air molecules, the ash layers can be tracked within the planetary boundary layer. This is not possible using ACL measurements alone as the volcanic ash signal is tainted by other aerosol types here. It is, however, rather difficult to determine unambiguously which ash layer structure observed by the ACL instrument can be related to the appropriate structures simulated by the model. Regarding the thin volcanic ash layer which is measured by the ACL instrument in a height between 7 and 9 km ASL at the 16 April 2010, around 06:00 UTC, this feature could be equivalent to the model prediction of ash in a height 30 of 6 km ASL at 7:00 UTC. In this case, the model would have performed a rather precise prediction with only one hour time lag and a two km vertical shift. But it is also possible that the predicted ash entrainment over the ACL station is equivalent to the ash-indicating ACL signals at around 12:00 UTC. In the latter case, the model prediction would be wrong by a time lag of about 6 hours which is insufficient for time-critical applications.
The qualitative comparison is currently limited to coordinates where the major fraction of scatterers are represented by both model and forward operator. There is, however, one scatterer fraction still missing on the present model runs for a comprehensive comparison: Other aerosol types than volcanic ash like anthropogenic emissions, mineral dust, soot, pollen, etc., are not included which leads to differences especially in the planetary boundary layer. We therefore cannot predict yet whether the strong ACL signal in the planetary boundary layer is related to background aerosol or errors of the model. To further investigate 5 this problem, future studies with several types of aerosols incorporated into the model will be helpful.

Quantitative comparison
A major purpose of the backscatter lidar forward operator is the capability to perform also quantitative comparisons of measurement and model output data. Unfortunately, such comparison is of limited validity in this case study due to the unknown ACL calibration as noted in Section 3.2.

10
Outside the volcanic ash layer, the forward operator returns an attenuated backscatter coefficient of 1 × 10 −7 m −1 sr −1 which is equal to the values of the ACL instrument after calibration. This would be expected as both temperature and pressure are rather precisely determinable and the scattering properties of air are represented by the empirical equations which are used for the forward operator. We therefore assume that the selected calibration factor is valid for this scenario.
Regarding the ash layer, however, the forward operator returns stronger signals inside the ash plume as well as a lower factor 10-15 lower than in the same height but without passing through the volcanic ash. The two-way transmission at this height is thus about 7-10 % which seems to be a too low value compared with the observations.
We therefore analyzed the effect on the attenuated backscatter coefficient if the model-predicted ash number densities are being reduced by factors of 10, 20 and 30. If the ash number density is reduced by factor of 20 (Fig. 16), similar maximum values of the attenuated backscatter coefficient are observed inside the ash layer for both the forward operator and the ACL 25 (5.0 × 10 −5 m −1 sr −1 ). This reduction of the number density also results in less extinction, and the two-way transmission has a minimum value of 70 % (plot not shown here) which is more realistic than the minimum value observed for the original dataset (8 %, see Fig. 14).

Conclusions
A backscatter lidar model capable of calculating both the extinction and backscatter coefficients was introduced. Detailed 30 studies concerning the scattering properties of particles and molecules were performed. Instead of assuming a lidar ratio for given particles, this forward operator allows for calculating the scattering properties even for mixtures of different particle types. The shown method allows for an efficient calculation of ACL profiles based on the output of atmospheric chemistry models.
The forward model was applied to the COSMO-ART model but the same approach can be used for any other aerosolchemistry-transport model. Data of a COSMO-ART ash-dispersion simulation was used to run the forward operator and perform both qualitative and quantitative comparison between the output of the forward operator and measurement data of an 5 automated ceilometer lidar (ACL) system.
A major challenge for setting up the forward operator in a given scenario are the scattering calculations, namely the calculation of the effective extinction cross section σ ext,R d ,k,λ and the effective differential backscatter cross section σ bsc,R d ,k,λ of all model-represented particle size and type classes.
The atmospheric gas mixture could be simplified to a uniform mixture of "atmospheric gas". Therefore, empirically deter- Due to uncertain refractive indices of the volcanic ash and the fact that volcanic ash particles show complex shapes, sensitivity studies have been performed to analyze the impact of different particle shapes on the effective extinction and backscatter cross section as well the lidar ratio. While the extinction cross section was only weakly sensitive to variable refractive indices 20 and particle shapes, the backscatter cross section was strongly sensitive to both. However, we found that the sensitivities reduce significantly, when averages over size classes are made. We expect that this very interesting result -described in this manuscript for the first time -will be very helpful for comparisons of modeled with measured backscatter lidar data.
Using these effective optical cross sections reduced the sensitivity of optical cross sections regarding the refractive index as well as the particle shape. But even after averaging, the relative uncertainty of the effective backscatter cross section exceeds 25 280% for uncertain refractive indices. This study also indicates the dependency of the forward operator on precise information about the particle's refractive index. Within a particle shape sensitivity study, we were able to resolve the relative uncertainty of each individual size class for the effective backscatter cross section. Assuming that volcanic ash consists of a mixture of particle shapes, we analyzed the relative differences between the reference particle shape and 11 particle shapes (6 types of ellipsoids and 5 types of cylinders). 30 The forward operator matches the lidar ratio values we find in literature (40 sr to values higher than 100 sr for volcanic ash (Kokkalis et al., 2013;Ansmann et al., 2010;Mortier et al., 2013)). On average, the lidar ratio is 61.17 sr which fits well to the literature findings. Comparing the lidar ratio values of the first two size classes with the lidar ratio values reported by Gasteiger et al. (2011b), a lidar ratio of less than 20 sr seems to be plausible for these particle size to wavelength ratios. The authors found even for irregularly shaped objects a lidar ratio between 5 sr and 20 sr at size parameters between 5 and 15 (equivalent 35 particle diameter at λ = 1064 nm would be 1.6 µm and 4.8 µm, respectively). We therefore assume that the calculation results of the backscatter lidar forward operator are valid and allow for both qualitative and quantitative comparison.
A comparison between ACL measurement and the model predictions used as input for the developed forward operator was shown. Similar structures were observed but some features were referenced to different time and height locations. From our analysis at the ACL station Deuselbach, some ash layer features were predicted quite precisely by the model, for example the Due to unknown calibration coefficients of the ACL system, a calibration constant η * was estimated by comparing the ACL 10 data with calibrated measurements at the same wavelength. Within quantitative comparisons between ACL measurements and the forward operator output, we found that the molecule signal of ACL and forward operator output were of the same order of magnitude which argues that the selected calibration factor was reasonable. Meanwhile, the ACL manufacturers have understood the importance of calibrated backscatter data and implemented technical solutions so that a similar effort as described in this study with data for the year 2010 became obsolete.

15
A comparison the volcanic ash signal led to the conclusion that the model predicted ash concentration was to be too high as the forward modeled attenuated backscatter coefficient within ash layers was 60 times higher and after attenuation 10 times lower than observed by the ACL. If the model-predicted ash concentration is manually reduced by a factor of 20, the forward modeled COSMO-ART predictions and ACL measurements became quantitatively similar. Such a reduction could be part of a simple particle data assimilation system helping to calibrate particle dispersion simulations before in-situ measurements are 20 available. It would of course be beneficial, if even better information on the refractive index and effective particle shape and aspect ratio of volcanic ash particles becomes available in the future.
Furthermore, we analyzed the contribution of each class to the total backscatter coefficient and to the total mass density for two sample cases. Regarding case 1 inside the volcanic ash layer, the classes 1, 2 and 3 were mostly responsible (94.8 %) for the calculated attenuated backscatter coefficient. As these classes contribute most to the forward modeled signal, the value of the 25 lidar ratio would also be expected to be dominated by their contribution, namely a value between 5.23 sr and 58.83 sr (see Table   3). Raman lidar measurements of the Eyjafjallajökull ash resulted lidar ratio values greater than about 40 sr at wavelengths of 355 nmand 532 nm (Großet al., 2012) which is within this range. Thus, the calculated values of both extinction and backscatter cross section as well as the lidar ratio seem to be plausible.
Assuming that the ACL calibration is valid, the model-predicted ash concentration was about 10 times higher than ob- 30 servable. There are, however, some error sources remaining which are: First, there are only molecules and the six volcanic ash classes represented while background aerosol is missing completely. Second, the ACL calibration is of limited precision.
Third, the contribution to the attenuated backscatter coefficient of ash size classes 4, 5 and 6 is relatively low even though these classes carry a large proportion of the mass. This relationship could rely on the ACL's wavelength which probably limits its sensitivity to particles larger than about 10 µm in diameter. Such results strengthen the importance of a joint use of observations 35 20 Atmos. Meas. Tech. Discuss., doi:10.5194/amt-2017-142, 2017 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 17 May 2017 c Author(s) 2017. CC-BY 3.0 License. and model output in combination with data assimilation in order to get the best state of the atmosphere with respect to aerosol distributions and properties.
Conclusively, we recommend further investigation in scattering calculations of non-spherical particles to get more realistic optical cross sections for the forward operator. A decrease of uncertainties related to the forward operator can be achieved by refractive index measurements at the exact ACL wavelength. Refractive index measurements are a basic aspect of the forward 5 operator as the optical cross sections can only be calculated if the aerosols' refractive index is known precisely. The model -and consequently the forward operator -has to represent more aerosol types, especially background aerosols, mineral dust, sea salt and soot as missing extinction near ground cause the forward operator to overestimate the signal from layers behind.
But also qualitatively, more scatterer size classes are required to also represent the fine fraction and very large particles in the atmosphere. One approach for a better representation of the natural size-spectrum of aerosols is the use of continuous number-10 size distributions which are aggregated from multiple distribution functions ("modal" approach). On the one hand, this already includes the size-averaging which is necessary for monodisperse size distributions. But on the other hand, the model delivers exact information about the outer margins, i.e. the number-density of the fine and the extreme coarse fraction which is currently not reproduced by model and forward operator in the Eyjafjallajökull case study.
In the context of international and probably intercontinental ACL networks, the creation of a scattering database for aerosols to be sensitized to data quality and reproducible measurement calibration. Otherwise, the backscatter data is nearby useless for any quantitative comparison or aerosol data assimilation approach. The uncertainties in both modeling and measurements, however, will also require sophisticated data assimilation algorithms not only for typical atmospheric variables but also for aerosol optical properties. Also a very good first guess of model simulations with respect to aerosol particles will be necessary so that more sources, types, and sinks will be to be included.    Table 3. Effective optical cross sections and average lidar ratio of atmospheric gas molecules and six volcanic ash size classes with their respective aerodynamic diameter. The effective optical cross sections were used by the forward operator in the Eyjafjallajökull case study.