A Backscatter Lidar Forward Operator for Particle-Representing Atmospheric Chemistry Models

State-of-the-art atmospheric chemistry models are capable of simulating the transport and evolution of particles and trace gases but there is a lack of reliable methods for model validation and data assimilation. Networks of automated ceilometer lidar systems (ACLs) could be used to fill this gap. These are already used for the detection of clouds and aerosols, providing :::::: provide : a 3D dataset of atmospheric backscatter profiles. However, as the aerosol number concentration cannot be obtained 5 from the ACL data alone, one needs a backscatter-lidar forward model to simulate lidar profiles from the model variables. Such an operator allows then : is ::::::: required :::::: which :::::: allows for a qualitative and quantitative model validation based on ACL data. In this work, we present ::: We :::::::: developed : a new backscatter-lidar operator which contains most of the microphysical properties of aerosol particles, discuss sensitivity studies and compare simulated with measured ACL profiles. A major challenge is :::::: forward ::::::: operator :::::: which :: is ::::: based :: on ::: the ::::::: distinct ::::::::: calculation :: of ::: the :::::::: aerosols’ ::::::::: backscatter :::: and :::::::: extinction ::::::::: properties. :::: The ::::::: forward 10 ::::::: operator ::: was ::::::: adapted ::: to the high sensitivity of the optical cross sections to the particle size and shape: A slightly different particle radius may lead to quite a large change of the scattering properties. As most particle size distributions are continuous in reality, the optical cross sections are averaged over certain size-intervals which also reduces the problematic and unrealistic sensitivity significantly. To calculate the attenuated backscatter coefficient, the size-dependent particle number concentration and the scattering properties of each particle type and size have to be simulated :::::::::::: COSMO-ART ::: ash ::::::::: dispersion ::::::::: simulation :: of ::: the 15 :::::::::::: Eyjafjallajökull :::::::: eruption :: in :::: 2010 ::: as : ::: due :: to ::: its :::::: impact :: on ::: the :::::: public ::::: sector : :::: such :::::: events :::: have ::::::: become :: a ::::: major ::::::::: motivation ::: for :::::: aerosol ::::::::: dispersion :::::::: modeling. While the particle number concentration is a ::::::: provided :: as : model output variable, the scattering properties have :: of :::: each :::::::: individual ::::::: particle :::: type ::: has : to be determined by extensive scattering calculations. As these scattering calculations require assumptions about ::::::::: concerning : the particle refractive indices and shapes :::: index :::: and ::::: shape, sensitivity studies were performed to estimate the uncertainties related to the particle properties as represented by the model system. The strong 20 sensitivity of the scattering characteristics to the particle radius was largely reduced by size-averaging algorithms. We focus on a case study of the eruption of the Islandic volcano Eyjafjallajökull from 20 March 2010 to 24 May 2010. The Consortium for Small-scale Modeling Aerosols and Reactive Trace gases (COSMO-ART) model of DWD (Deutscher Wetterdienst) and KIT (Karlsruhe Institute of Technology) was used during this event for ash-dispersion simulation over Europe. For the forward model ::::::: assumed ::::::: particle ::::::::: properties. ::::::::: Therefore, :::::::: scattering :::::::::: calculations ::: for :::::: several ::::: types :: of :::::::::::: non-spherical ::::::: particles :::::::: required ::: the 25


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. Therefore, DWD decided to start a dedicated project on backscatter lidar forward operators for validating aerosol dispersion models using available remote sensing measurement data.
Aerosol particles in general play an important role in the Earth's weather and climate (Lohmann and Feichter, 2005) as well as for our every day life in terms of air quality (Pöschl, 2005). The role of aerosols in the Earth's atmosphere is commonly separated into direct and indirect effects (Boucher et al., 2013). Direct aerosol effects are due to the interaction with electromagnetic waves and the resulting consequences on the energy balance. The follow-up processes of direct effects are called indirect effects which are, for example, effects on cloud formation, water cycle, and plant growth. In recent years, numerical atmospheric 5 chemistry models have become capable to simulate the evolution, transport and interaction of aerosols in the atmosphere. But observations to validate these simulations are rare.
The simulation of aerosols can be separated in components: sources, sinks and evolution as well as transport and conversion processes. While the transport processes are similar for any aerosol type with only sedimentation being potentially different, the sources and sinks have to be described individually. For example, the emission rate of mineral dust from deserts depends on the wind speed near ground (Vogel et al., 2014a). Processes leading to sinks and conversation are washout, chemical reactions, 20 and particle aggregation. The transport processes mainly depend on the model dynamics (wind speed, wind direction) and sedimentation properties of each individual aerosol type. In consequence, the validation of aerosol transport model simulations with observations allows not only for the validation of the sources and sinks but serves also to validate the model dynamics. If such validation can be done automatically, data assimilation is the next step to continuously converge simulation and reality.
Compared to other atmospheric variables such as pressure and temperature, there are -to the best of our knowledge -no 25 comprehensive measurement systems in use so far validating aerosol simulations operationally.
The only remote sensing technique which is ::::: Lidar ::::: (light :::::::: detection ::: and :::::::: ranging) :: is : capable of providing information on atmospheric particles with high temporal and spatial resolutionis lidar (light detection and ranging). 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., 30 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 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;35 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 (Shipley et al., 1983).
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., 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 10 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.
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 15 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.
There are already several backscatter lidar forward operators available or in development. At ECMWF (European Centre for Medium-Range Weather Forecasts), a lidar forward operator was developed and tested in an ice cloud scenario (Benedetti, 2009), based on assumptions concerning the lidar ratio S lidar . Newer implementations of the ECMWF model allow for sim- 20 ulating the backscatter coefficient measurement of the satellite-mounted CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization) at two wavelengths 532 nm and 1064 nm for spherical sea salt and mineral dust particles (Morcrette et al., 2009 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 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.

35
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 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 5 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 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 10 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.
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.

15
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 20 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 25 photon counts. Here, we use photon counts per laser pulse unless otherwise noted.
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., 2010b). The missing overlap limits the capability to measure in near range but has no effect on heights where full overlap has 30 accomplished. Heights where 0 < O(z) < 1 can be overlap-corrected if the device specific overlap function is known.
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 . 5 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 10 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 15 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 atmospheric scatterer data.

20
If there are q s different types of scatterers in an illuminated volume, the total extinction coefficient α λ (z) and the total backscat- 25 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 .
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 5 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).

10
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.
The molecule number density N mol (z) is related to the ideal gas law 20 where p is the model prediction of the atmospheric pressure given in Pascal (Pa), T is temperature given in Kelvin (K), and k 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 had to extrapolate ::::: simply ::::::::::: extrapolated : the values to the ACL wavelength in the case study

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; again, with sizes not much smaller or larger than the wavelength. The T-matrix algorithm is a tool for computing scattering by 5 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 T-matrix instead of the Mie approach and by another two orders of magnitude when using the DDA instead of the T-matrix 10 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 in this work. 15 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 20 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 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 5 different particle radius may lead to quite a large change of the scattering properties. We present in the following a suggestion 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 10 backscatter cross sections with 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 15 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 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.

5
This simple solution allows for calculating α par,λ (z) and β par,λ (z) by just solving a few multiplications and summations resulting in a minimal demand of computing time.

Lidar Ratio and Two-Way Transmission
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). The lidar ratio S lidar (z) is calculated by 15 Even if 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 that are capable of measuring this parameter like a Raman lidar.
T λ describes the two-way transmission 20 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, 25 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.
Conclusively, our implementation of a backscatter lidar forward operator allows for an efficient calculation of ACL profiles based on the output of atmospheric chemistry models. The biggest challenge for setting up the forward operator in a given scenario are scattering calculations, namely the calculation of the effective extinction cross section σ ext,R d ,k,λ and the effective 5 differential backscatter cross section σ bsc,R d ,k,λ of all represented particle size and type classes. Due to its importance, we dedicate Sect. 4 to sensitivity studies of the optical cross sections regarding the particle size, refractive index, and shape. Prior to that, we introduce the case study for which the forward operator was applied and compared to ACL measurements. 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 15 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 20 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 25 ash cloud over Germany with this data set was made by Flentje et al. (2010b).
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 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 (Flentje et al., 2010c) resulting in an emitted photon number per pulse of about 4.28 × 10 13 . The diameter of the receiving telescope is 100 mm (Flentje et al., 2010c) which results in A tel = 78.54 cm 2 . The vertical resolution ∆z is 15 m over the complete profile. The overlap function O(z) was set to 1 which implies that ranges below about 1500 m cannot be used 10 reliably for comparisons with the forward operator.
Unfortunately, the instruments provided no calibrated measurement data at that time. As the true system efficiency η λ and the 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 15 height of 2 km ASL (about 5 × 10 −6 m −1 sr −1 ) to determine :::::: estimate : the maximum attenuated backscatter coefficient which should be measured by the ACL system inside the volcanic ash cloud. The second calibration criterion was the forward modeled attenuated backscatter coefficient outside the volcanic ash layer in 3 km ASL which is dominated by molecule scattering and has a value of about 1 × 10 −7 m −1 sr −1 , see Fig. 15. This pragmatic approach is of very limited precision and only required for uncalibrated ACL measurements. As most ACL networks have been extended by automatic calibration methods after the 20 Eyjafjallajökoll eruption, this calibration "by hand" ::: 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 averaging. This also improved the signal-to-noise ratio of the ACL data.

Ash Transport Simulation of COSMO-ART
COSMO-ART was set up by the DWD in collaboration with Karlsruhe Institute of Technology (KIT) for an ash-dispersion 25 simulation of the volcanic emissions during the eruptive phase of Eyjafjallajökull in spring 2010 (Vogel et al., 2014b). The model domain had a horizontal resolution of 7 km and 40 height layers. The height layers have a variable thickness ranging from several meters near ground to a layer thickness of about 3 km in 22 km height above ground level. The spatial extend and model set-up is equal to the COSMO-EU domain as it is used at DWD for operational weather forecasting. A more general description of the model run is given by Vogel et al. (2014a). For this study, we used the 78-hour forecast beginning at 15 April 2010, 00:00 UTC, which includes volcanic ash emission data since 14 April 2010, 06:00 UTC.
Volcanic ash was represented by six discrete size classes with a spherical shape and aerodynamic diameters of 1 µm, 3 µm, 5 µm, 10 µm, 15 µm and 30 µm. For each class, a number concentration was predicted by the model. Particles within a class were treated as equal, i. e., they had the same size, shape, and complex index of refraction (monodisperse classes). We therefore 5 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.
A list of model variables used for the forward operator is given in Table 1. While the p number density is taken directly from 10 model output, the molecule number density of standard 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)  Electron microscope images from the same study revealed that the volcanic ash particles were sharp-edged with a complex 20 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 % (Lathem et al., 2011). A growth of 5 % does not change the scattering properties significantly in relation to the size averaging 25 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.

30
The representation of the particles by the model is clearly simplifying, so we must study the effect of these simplifications on the scattering of laser light when developing a forward operator. For a lidar forward model, sensitivities of the backscatter cross section are critical because the received signal intensity is linearly coupled to the backscatter cross section and, consequently, to the attenuated backscatter coefficient.
Prior studies already showed the complexity of non-spherical scattering calculations but there is no universal solution to the problem available. Gasteiger et al. (2011b) use Discrete Dipole Approximation (DDA) to calculate the scattering properties of complex shaped particles but the analysis is limited to size-parameters up to 20.8 (which results in a maximum particle 5 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-ART available. We therefore treat the volcanic ash as spherical objects with given optical properties (see 3.4) but analyze and 10 discuss the effect of uncertain volcanic ash properties in the following.
It must be noted that these studies are required for nearby ::::: 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. The refractive index measurements by Schumann et al. (2011) were not performed for the wavelength of the ACL systems as explained in section 3.4. Therefore, we have to assume the refractive index we use as reference as well as an interval of uncertainty. Schumann et al. (2011) take a refractive index of 1.59 -0.004i for their medium "M" case study and therefore, we also used this value as reference for our study. Our uncertainty intervals of real and imaginary parts were chosen according to 25 the range of measured values at 630 nm and 2000 nm, namely a real part range of 1.54 to 1.64 and an imaginary part range of -0.006 to -0.002. To get an estimate of the overall refractive index sensitivity for such particles, we decided to extend the range of analyzed refractive indices to real parts between 1.49 and 1.69 using increments of 0.001 and to imaginary parts between -0.011 and -0.001 using increments of 0.00005. Consequently, the total element number of one LUT is 4.0 × 10 7 .

Prerequisites
These look-up tables were the base for the subsequent sensitivity studies :::::::: 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 extinction cross section σ ext and the effective backscatter cross section σ bsc .
An overview of the refractive index sensitivity of the effective optical cross sections is given by Fig. 7 which shows the relative errors 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 10 to 230% for the backscatter cross section.

T-matrix Particle Shape Sensitivity Study
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 15 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 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 20 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.
Consequently, an ellipsoid with EPS=1 is a sphere, prolate objects have an EPS<1 and oblate objects have an EPS>1.
For the current state of the forward operator research, the spherical shape has to be used as reference even if the real volcanic ash particlesare known to be fractal and complex shaped. One of the reasons for doing so is the fact that there is currently 30 no appropriate shape representation schemes for volcanic ash available. From the findings of Rocha-Lima et al. (2014), the average aspect ratio of volcanic ash is known but not the particular shape. As we have shown, the backscatter cross section of cylinders differs tremendously from ellipsoids and from spheres etc. so there is currently no indication for using a given aspherical shape as reference because doing so could yield in even higher errors than assuming a spherical shape.

25
The above results, however, give us valuable insight into the uncertainties of using Mie calculations for non-spherical particles. The particle shape sensitivity is reduced by size-and shape-averaging methods but the resulting effective backscatter cross section of the mixture may vary by a factor of ±2. These findings are important for interpreting the results of the forward operator. The assumption of spherical volcanic ash particles results in an overestimation of the backscatter cross section of size classes 1, 2, and 3 by factor 2-3 which results in an equivalent underestimation of the respective lidar ratio values.

Scattering Properties of Volcanic Ash Used Within the Forward Operator
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.

5
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 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 10 with height (see Eq. (24)).
The forward operator yields a lidar ratio for volcanic ash of between 15 sr and 80 sr which differs to lidar ratio values we find in literature (40 sr to 100 sr for volcanic ash (Kokkalis et al., 2013;Ansmann et al., 2010;Mortier et al., 2013)). This difference may be a consequence of assuming spherical particles. But also for spherical particles we would expect the given lidar ratio as the extinction cross section increases exponentially with the particle size while the backscatter cross section is 15 influenced by interference effects (see. Fig. ?? and Fig. ??). This is also represented by the effective optical cross sections which are used by the forward operator (see. Table. 3). Comparing these values with the lidar ratio findings by Gasteiger et al. (2011b), a lidar ratio of less than 20 sr −1 seems to be plausible. The authors found even for irregularly shaped objects a lidar ratio between 5 sr −1 and 20 sr −1 for size parameters between 5 and 15 (equivalent particle diameter at λ = 1064 nm would be 1.6 µm and 4.8 µm, respectively). Due to the different refractive indices and particle shapes, this comparison is no validation 20 for our findings but shows that the result of our analysis is close to the reference calculations.
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 contributing by 10 % each to the total mass density and class 1 is nearby irrelevant ::: very :::: low : for the total mass density.

Qualitative Comparison
A qualitative comparison allows for the identification of common and different structures between the measured and simulated A comparison of ACL measurement and COSMO-ART simulation with an applied forward operator is shown in Fig. 15. Due 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 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 30 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 16th of April :: 16 ::::: April :::: 2010, around 06:00 UTC, this feature could be equivalent to the model prediction of ash in a height 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 5 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 overcome ::::: further :::::::::: investigate this problem, continuous data assimilation of aerosols and the simulation of all aerosol particle sources in the model are required ::::: future :::::: studies :::: with :::::: several ::::: types :: of ::::::: aerosols ::::::::::: incorporated ::: into ::: the :::::: model ::: will ::: be :::::: helpful.

Quantitative comparison 10
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.
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 15 are rather precisely determinable and the scattering properties of air are represented by the empirical equations which are utilized by :::: 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 transmission behind the cloud compared to the ACL measurement. The maximum values of the attenuated backscatter coefficient returned by the forward operator (about 6.0 × 10 −4 m −1 sr −1 ) is 60 times higher than the maximum values reported by 20 the ACL (about 3.0 × 10 −5 m −1 sr −1 ). The minimum attenuated backscatter coefficients calculated by the forward operator (1.0 × 10 −9 m −1 sr −1 ) are about 10 times lower than observable on the ACL plots (1.0 × 10 −8 m −1 sr −1 ), excepting pixels which became negative due to noise. Behind :::::: Above the cloud (in about 12 km ASL), the attenuated backscatter coefficient is by about factor 10-15 lower than in the same height but without passing through the volcanic ash. Due to the definition of the ::: The : two-way transmission (eq. 24), the two-way transmission at this height can be estimated to a value of :: is ::: thus ::::: about : 7-10 25 % which is a rather low value for the already distributed volcanic ash layer ::::: 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 30 (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).
A backscatter lidar model capable of calculating both the extinction and backscatter coefficients was introduced. Detailed studies concerning the scattering properties of particles and molecules were performedwithin this study. Instead of assuming a lidar ratio for given particles, this forward operator allows for calculating the lidar ratio :::::::: scattering :::::::: properties : even for mixtures of different particle typeswhich is an outstanding feature with many fields of application in the future. 5 The forward model was applied to the COSMO-ART model but the same approach can be used for any other aerosolchemistry-transport model. We used data of a COSMO-ART ash-dispersion simulation to run the forward operator and perform both qualitative and quantitative comparison between the output of the forward operator and measurement data of an automated ceilometer lidar (ACL) system.
The atmospheric gas mixture was simplified to a uniform mixture of "atmospheric gas". Therefore, empirically determined 10 scattering equations were used to calculate the optical cross sections of this mixture for the given laser wavelength. From the model-predicted values of temperature and pressure, the molecule number-density and finally the molecule extinction and backscatter coefficients were calculated.
Regarding particle scattering, extensive scattering calculations were performed to create look-up tables of optical cross sections. The range of particle sizes was selected according to the volcanic ash size classes used by COSMO-ART (six monodis- 15 perse classes with diameters of 1 µm, 3 µm, 5 µm, 10 µm, 15 µm and 30 µm). The range of refractive indices were adapted according to in-situ measurements of Schumann et al. (2011). As reference shape, we assumed the volcanic ash to be spherical.
Furthermore, we analyzed the contribution of each class to the total backscatter coefficient and to the total mass density ( Fig. ??) 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, 10 the value of the lidar ratio would also be expected to be dominated by their contribution, namely a value between 4.64 sr −1 and 8.49 sr −1 ::::::::: 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 35 sr −1 (Groß et al., 2012). This discrepancy between the forward modeled and measured lidar ratio values are related to several reasons. First, the measurements are not for the same wavelength, i.e. measurements are only available for :::::: 40 sr at wavelengths of 355 nmand 532 nm but not for the ACLs wavelength (1064 nm). Second, the forward modeled lidar 15 ratio values of classes 1, 2 and 3 are known to be underestimated by factor 2-3 due to the assumption of sphericity (see section 4). Third, the volcanic ash size classes defined by the model are maybe too coarse. Size class 1 represents ash particles with a diameter of 1 µm so lower margin of the size-averaging algorithm of this class was set to its half, namely 0.5 µm in diameter.
This also implies that the majority of the fine fraction is not represented by the model and also not by the forward operator. It should be mentioned here that comparisons between measured and forward modeled lidar ratios require further investigation 20 such as the calculated of the weighted lidar ratio where the contribution to the total lidar ratio of each class depends on the respective contribution to the backscattered signal or on the respective number concentration. Lidar ratio simulations of other scientists, however are exactly in the same scope as the forward modeled lidar ratio values of our study. Gasteiger et al. (2011b), for example found lidar ratio values of less than 20 sr −1 for absorbing and non-absorbing dust particles with a radius close to the laser wavelength (size parameters between 4 and 10. We therefore conclude that the forward modeled values and uncertainty 25 intervals of the :::::::::::::::::::: (Groß et al., 2012) which :: is :::::: within ::: this :::::: range. ::::: Thus, ::: the ::::::::: calculated ::::: values ::: of :::: both ::::::::: extinction ::: and : backscatter cross section are reasonable. Especially the ash size classes 1, 2 and 3 seem to suffer from the assumption of a spherical shape because respective realistic lidar ratio values are at the extreme end of the calculated uncertainty interval. Further studies are therefore required to find a better reference shape than spherical. Scattering calculations for spherical particles are possible with nearby no restrictions related to the wavelength to particle size ratio. For non-spherical shapes, the used t-matrix code we 30 used was limited to particle sizes lower than 10 µm for weakly asymmetric and to particle sizes lower than 3 µm for strongly asymmetric particle shapes. : 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 observable. 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.

35
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 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. 5 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 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 10 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 types ::: size :::::: classes are required to cover all atmospheric features measured by ACL systems :::: 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-size distributions which are aggregated from multiple modes ::::::::: distribution :::::::: functions : ("modal" approach). On the first ::: one hand, this already includes the size-averaging which is 15 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 would be desirable. A central database can increase the development rate and flexibility of current lidar forward operator 20 implementations. The ACL networks themselves are only useful for aerosol research and data assimilation if the calibration is performed automatically and transparent. For future lidar measurement networks, the number of high spectral resolution lidar (HSRL) systems and Raman lidar systems could potentially increase and allow for the assimilation of the extinction coefficient and the backscatter coefficient directly. As many ACL devices are operating a proprietary firmware, the manufacturers have to be sensitized to data quality and reproducible measurement calibration. Otherwise, the backscatter data is nearby useless 25 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.
Time-height cross section of backscatter coefficient, two-way transmission, total extinction coefficient and lidar ratio, calculated by the forward model based on COSMO-ART output at the station Deuselbach (West Germany). The forward model represents clean air molecules and volcanic ash particles (no clouds, rain, fog, background aerosol or other scattering objects).   Fig. 15 with a decreased volcanic ash number density. For this purpose, the ash number density predicted by COSMO-ART was reduced by factor 20 before applying the forward operator.