Multiple scattering by aerosols as seen from CALIPSO - a Monte-Carlo modelling study.

The impact of multiple scattering (MS) by aerosols on satellite-borne lidar measurements is studied by Monte-Carlo radiative transfer simulations. A total of 48 aerosol scenarios are considered. We find that the frequently used MS correction factor can be parameterized as a function of aerosol size and aerosol optical depth. Its dependencies on vertical distribution and total optical depth can be treated as a random error. We illustrate the use of our parameterization by considering an episode of high sea salt concentrations over the ocean. Neglecting MS, or using a constant value of the MS correction factor, can introduce a negative bias in the computed backscattered power that exceeds the random error in our approach.


Introduction
Multiple scattering poses a challenge to the interpretation of Lidar [1] and Radar measurements [2]. The return signals are commonly interpreted by use of the lidar equation, which is an approximate description of the radiative transfer problem based on neglecting the effect of multiple scattering. The validity of this approximation depends on the geometry of the observing system, the density of the scattering medium, and the size of the scatterers. In general, the field of view (FOV) of the instrument's telescope and the beam divergence need to be sufficiently small. Even when these conditions are met, for space-borne lidar instruments, the beam footprint often has a diameter comparable to the photon mean free path. In such case multiple scattering cannot be neglected (e.g. [3]). Figure 1, diagram a, illustrates the basic principle. A laser pulse, emitted from a satellite-borne instrument, is scattered at a location a distance R away from the instrument, corresponding to an altitude z measured from the bottom of atmosphere (BOA). Part of the energy is scattered in the exact backscattering direction and detected by the instrument's telescope. The laser pulse will be received at a time t = 2R/c after emission, where c is the speed of light. However, part of the emitted energy can also travel along optical paths as in diagram b. If the path length for this process agrees with R within the instrument's vertical resolution ∆R (indicated by the shaded band), then the contribution of this process to the backscattered signal will be correctly attributed to scattering at altitude z. However, this light path is not taken into account when we model the lidar return signal by use of the single-scattering approximation. Thus, a single-scattering model will usually underestimate the backscattered energy.
We further may have to consider optical paths as in diagram c. This illustrates a multiplescattering light path within the same range of pathlengths 2R ± ∆R as in the other two diagrams. But in diagram c the scattering takes place at altitudes much higher than z. Thus the contribution of such processes will be incorrectly attributed to scattering at altitude z. For brevity, let us label the light paths shown in diagram b as regular multiple scattering paths, and those in diagram c as irregular multiple scattering paths. More specifically, we define regular multiple scattering to comprise those paths for which the number of scattering events is at least 2, for which the total path length lies in the range 2R ± ∆R, and for which the lowest altitude along the light path lies in the range z ± ∆R/2, where z = R 0 − R, and R 0 denotes the distance of the instrument from the ground. Irregular multiple scattering comprises those corresponding multiple scattering paths for which the lowest altitude along the light path is higher than z + ∆R/2. Our main concerns in this paper are to (i) establish that the contribution from irregular multiple scattering is negligibly small; and (ii) quantify and parameterize the contribution of regular multiple scattering to the backscattered energy.
Regular multiple scattering will be dominated by second order scattering. All detectable second order scattering paths necessarily lie completely within the detector's FOV. Thus, the contribution of regular multiple scattering to the backscattered energy will increase when increasing the telescope's FOV. Also, the frequency of such paths will be inversely proportional to the angular width of the particles' diffraction peak in the phase function, which, in turn, is inversely proportional to the average size of the particles. Thus, the backscattered energy due to regular multiple scattering increases with increasing particle size.
The effect of multiple scattering by haze and cloud droplets on lidar return signals has been studied since the early 1970s. Both theoretical models [4][5][6][7][8] and experimental studies [9][10][11][12] have been reported. Later, multiple scattering was systematically exploited as an additional source of information. To this end, methods have been developed to measure the effect of multiple scattering by varying the FOV of the receiver [13][14][15][16] and to retrieve the extinction coefficient and effective particle radius [17,18]. This so-called multiple-field-of-view (MFOV) method exploits the dependence of regular multiple scattering on FOV and particle size.
While the majority of studies focus on the backscattered energy in the presence of clouds, there are several studies that discuss the impact of multiple scattering on depolarization (e.g. [19][20][21]), as well as on multiple scattering of lidar signals through aerosol media (e.g. [21,22]). Dedicated studies of space-borne lidar instruments have been reported in [23,24], and in [3,25]. The latter two papers discuss multiple scattering in CALIPSO retrievals for both cloud and aerosol layers.
Various theoretical models have been developed to describe the effects of multiple scattering on lidar return signals (e.g. [26][27][28][29]). A review of different approaches is given in [1]. In the following, we will rely on the use of a Monte-Carlo (MC) model. MC models are easy to implement and to set up for complex instrument specifications, while invoking few approximate assumptions. A disadvantage is that they often are computationally demanding.
We focus on the attenuated backscattering coefficient at a wavelength of 532 nm. Our simulations use the instrument specifications of CALIOP (Cloud-Aerosol Lidar with Orthogonal Polarization) on platform CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations). The main motivation is to account for multiple scattering in an observation operator (forward model) that can be used in a chemical data assimilation system. In general, an observation operator H maps a vector of model variables x to observable parameters H (x). (The vector x contains for each position in the discretized model grid the size distribution and chemical composition of the aerosols.) The observation operator for the attenuated backscattering coefficient is given by Here β m denotes the molecular backscattering coefficient, β a (x) represents the observation operator for the aerosol backscattering coefficient, and the transmission T 2 accounts for the attenuation of the laser pulse through the medium to the scattering location R and back to the source. In the single-scattering approximation, the transmission term is given by where τ(x; R) denotes the observation operator for the optical depth from the source to the point of scattering a distance R away from the source. The factor 2 in the argument of the exponential of the transmission term comes from the fact that the light ray travels twice the distance R on its way through the medium and back to the source. Platt [9] proposed an empirical correction to account for multiple scattering by reducing the optical depth by a factor η. This factor lies usually between 0.5 and 1. Thus, the transmission term becomes In ice clouds the factor η can be dependent on temperature [23]. This is because multiple scattering is dependent on particle size, and ice crystal size in ice clouds is, in turn, temperature dependent.
In this paper, we perform MC simulations for different aerosol scenarios. We model the multiple scattering correction factor η for CALIOP observations at a wavelength of 532 nm under different conditions. In previous publications η has been modelled for use in CALIOP retrievals [3,25]. Here, we want to consider the problem from the point of view of chemical transport modelling and data assimilation (CTM-DA). There are two main reasons for revisiting this topic from a CTM-DA perspective.
• In remote sensing retrieval problems the optical depth of the aerosol layer is unknown a priori. It therefore makes sense to parameterise η as a function of the range ∆r that the pulse has travelled into the aerosol layer, as has been done in [3,25]. By contrast, a CTM provides us with a priori information about the optical depth along the optical path. Therefore, for applications in CTM-DA it seems promising to parameterise η as a function of optical depth τ.
• The CALIOP retrievals use six standard aerosol classes. For each of these classes a separate parameterisation η(∆r) is given [3,25]. By contrast, in a CTM the aerosols are described by their size distribution and chemical composition, which is much more detailed than a classification into six aerosol types. For instance, "sea salt" and "dust" are among the six aerosol classes. However, in a CTM each of these aerosol types can be encountered with highly varying sizes. In view of the strong size dependence of MS in lidar observations, we aim for a parameterisation η(τ, r eff ), in which the effective radius of the aerosols, r eff , enters as a continuous parameter.
Our methods are introduced in Sect. 2. Results are presented and discussed in Sects. 3. and 4., respectively. Concluding remarks are found in Sect. 5.

Monte-Carlo radiative transfer model
We calculate the energy distribution, emitted from the laser source by means of a forward radiative transfer model. Since our intention is to trace individual light rays and scattering events as detailed as possible, a Monte Carlo radiative transfer model is the obvious choice. Our calculations are based on the GRIMALDI Monte Carlo model [30,31]. Relevant parameters for calculating the radiation field are (i) extinction coefficient, (ii) single scattering albedo and (iii) scattering phase function. These parameters are defined for each model voxel. Positions of light-matter interactions are determined by interpreting the transmission as a probability to travel an optical depth τ d (integrated extinction along the path of the light ray) without scattering or absorption. This leads to: where is taken from a series of uncorrelated (pseudo-random) numbers with ∈]0, 1[. The pseudo-random numbers are provided by a routine based on the ran0 random number generator from Numerical Recipes [32]. This random number generator utilizes bit shifting and features a period of 2 31 − 2. This allows for experiments with a large number of scattering processes. Absorption is taken into account by reducing the energy of a test-ray, for each extinction event, according to the local single scattering albedo ω 0 . The scattering phase function is taken as a probability distribution function (PDF). This allows us to determine a new beam direction by mapping another pseudo-random number to a map of scattering angles. This map is generated in 3 steps: First the scattering phase function is inverted. Then it is transformed into a cumulative function. Finally, the PDF is scaled to the range of pseudo-random numbers ]0, 1[. Given a sufficient number of scattering processes, the local scattering phase function is reproduced. The final radiation field is obtained by collecting weights at each location of interest (in this case CALIOP's detector) and scale them with the intensity of the incoming radiation (CALIOP's laser pulse). The error is only dependent on the number of contributing light-rays. For fields of irradiance the tracing and collecting of weights works very well. For the calculation of local radiance, however, it is computationally ineffective, especially for consideration of detectors with small aperture angles. In this case only a small number of light-rays contributes to the signal, and it takes a large number of sample rays to reach a certain level of accuracy.
To use the computation time in a more efficient way, we added the local estimation method [7,33] as an option to GRIMALDI. In the original approach light rays are traced, and if no detector is hit, the information is lost. By contrast, the local estimate approach uses the phase function to sample radiance contributions from each scattering event. The weight or energy of the light ray is modified according to the probability distribution of the scattering phase function, and according to the attenuation along the path between the scattering event and the detector.
We further implemented counters to keep track of the number of scattering events along a light path, the path length, and the minimum altitude along a path. This allows us to keep track of how much of the backscattered power is due to single-scattering, regular multiple-scattering, and irregular multiple scattering events.
Finally, we implemented the option to simulate a LIDAR-like light source and detector in GRIMALDI. The user chooses the instrument and detector specifications. For this study, we used the specifications of CALIOP. More specifically, we set (i) altitude of the instrument 701 km; (ii) receiver telescope 1.0 m diameter; (iii) field of view 130 µrad; (iv) beam divergence 100 µrad; (v) vertical resolution 30 m.

Aerosol scenarios
We compute aerosol optical properties for a wavelength of λ =532 nm. The aerosols are assumed to be homogeneous spheres. The computations have been performed with a standard Lorenz-Mie program developed by Mishchenko et al. [34]. A complex refractive index of m = 1.57 + 0.003i is assumed. This is a generic value that can represent moderately absorbing mineral aerosols (e.g. [35]). Secondary organic aerosols are somewhat less absorbing at visible wavelengths. For instance, toluene has a value of m = 1.55 + 0.0008i [36]. However, secondary aerosols are often internally mixed with small amounts of black carbon (e.g. [37]). Assuming a typical refractive index for soot around m = 1.95 + 0.79i [38], an internal mixture of toluene containing 4 % black carbon would have an effective refractive index of m = 1.56 + 0.003i (using the Maxwell-Garnett mixing rule [39]). Thus, our choice of the refractive index can also be seen as representing the effective refractive index of internal mixtures of secondary organic aerosols with small amounts of black carbon.
We further assume a mono-modal log-normal size distribution given by where r denotes the particle radius, σ g denotes the geometric standard deviation, and r 0 represents the mean radius. Size-averaged optical properties are computed by numerically integrating over 10,000 discrete sizes in the range from 0-50 µm. We consider four different values of r 0 , which are listed in Table 1. These values brace the range from fine mode to coarse mode aerosols. The value for σ g lies in the range of typical values of both accumulation and coarse mode dust aerosols under background conditions, as well as coarse mode particles under heavy winds and dust storms (e.g. [40]). The table also lists the effective radius which is the mean radius of the size distribution weighted with the projected area of the particles. For particles much larger than the wavelength, this is a convenient size measure in the context of aerosol optics, because the extinction and scattering cross sections are proportional to the projected area of the particles. For particles of sizes comparable to the wavelength, it is still a useful proxy for characterizing the average size of an ensemble of particles. For each size distribution, we consider six different values of aerosol optical depth τ * a in the troposphere, namely, τ * a =0.05, 0.2, 0.5, 1.0, 2.0, and 3.0. These values cover the range from very low background concentrations (e.g. [41]) to extremely high aerosol loads during episodic events (e.g. [42,43]). Note that τ * a refers to the aersol optical depth of the entire atmospheric column in the vertical direction, i.e., from the top of atmosphere (TOA) to the bottom of atmosphere (BOA) at the surface. By τ a (z) we denote the vertical aerosol optical depth from TOA to altitude z. Thus, τ * a = τ a (0). For each size distribution and for each optical depth, we consider two vertical distributions of aerosols. In the first case, the aerosol layer extends from 0-3 km altitude; in the second case, the vertical extent is 0-6 km. In either case, we assume that the aerosols are evenly distributed in the vertical direction.

Multiple scattering correction in the lidar equation
Consider an optical path as in Fig. 1a. Part of the emitted energy is backscattered at an altitude z, which lies a distance R away from the lidar instrument. Let us assume that the vertical resolution of the receiver is ∆R. Then the energy that is returned to the detector from points in the interval R − ∆R/2 ≤ R ≤ R + ∆R/2 is given by where E 0 is the laser pulse energy, and A is the area of the receiver telescope. In the singlescattering approximation, the transmission is given by Eq.
(2). The optical depth τ at position R is obtained by integrating the extinction coefficient k ext (dimension L −1 ) according to Thus, the optical depth at r = R and, therefore, the transmission and E(x; R) depend on the aerosol state vector x throughout the column from r = 0 to r = R. The extinction coefficient is the sum of the corresponding extinction coefficients of molecules and aerosols, i.e.
Thus, the aerosol optical depth is given by and similarly for the molecular optical depth τ m . The lidar equation given in (7) in conjunction with (2) only accounts for optical paths as in Fig. 1a. Optical paths such as in Fig. 1b are treated as extinction events in this approximation, even though these rays are detected by the receiver. To account for this error, we follow the approach in [9], which is based on introducing an empirical correction to the transmission term according to Eq. (3) The lower the correction factor η in this equation, the more strongly the return signal is enhanced by multiple scattering.
In GRIMALDI, we keep track of how much of the received power originates from singlescattering paths (Fig. 1a), from regular multiple-scattering paths (Fig. 1b), and from irregular multiple-scattering paths (Fig. 1c). Let us denote these contributions by For computing η, we assume that E ms, irr /E 1 (this will be verified in the next section). Then we can obtain η from the Monte-Carlo simulations as follows: which yields 3. Results Figure 2 shows computational results for three out of the 48 aerosol scenarios we studied. The mean radius and aerosol optical depth are r 0 = 0.2 µm, τ * a =1.0 (left column), r 0 = 2.0 µm, τ * a =1.0 (second column), and r 0 = 2.0 µm, τ * a =0.2 (right column). In all three cases the aerosols are evenly distributed between z =0-3 km altitude. The top row shows that fraction of the pulse energy that is backscattered to the receiver as a function of altitude z (measured from BOA upwards). The profiles show the total backscattered energy (black), the single-scattering contribution (green), the regular multiple-scattering contribution (red), and the irregular multiple scattering contribution (blue). The latter is indistinguishable from the vertical axis, indicating that irregular multiple scattering makes a very small contribution to the total backscattering signal. The relative contributions of regular multiple scattering, E MS, reg /E, and of irregular multiple scattering, E MS, irr /E, are shown in the second and third row, respectively. The fraction of irregular multiple-scattering is always lower than 0.005.
The profiles in the top row display the characteristic pattern of lidar return signals in the presence of aerosols. Above the aerosol layer, very little energy is backscattered. There is a small increase in signal strength with decreasing altitude owing to an exponential increase in air density. As the radiation encounters the aerosol layer at 3 km altitude, there is a sharp increase in the backscattered signal, followed by a decline near the surface due to extinction of the incoming as well as the backscattered radiation through the aerosol layer. In all cases, the simulated single-scattering signal shows very good agreement (not shown) with the lidar equation given in Eq. (7) in conjunction with Eq. (2).

Dependence of multiple scattering on aerosol size
For small particles with r 0 = 0.2 µm (top left) the single-scattering signal (green) agrees relatively well with the total backscattering signal (black). The regular multiple-scattering contribution (red) is small. However, as we increase the particle size to r 0 = 2.0 µm while keeping τ * a fixed (top center), the regular multiple-scattering contribution increases significantly at the expense of the single-scattering signal. Near the surface, the multiple-scattering signal even dominates over the single-scattering signal. This is also reflected in the ratio E MS, reg /E -compare the left and center diagrams in the second row. This effect can be explained by the differences in the aerosol phase functions. Large aerosols have more asymmetric phase functions with pronounced diffraction peaks as compared to small aerosols. Thus, large aerosols scatter more radiation in the near-forward direction than small aerosols, which increases the probability of detecting radiation that has travelled along regular multiple-scattering light paths.

Dependence of multiple scattering on total aerosol optical depth
Next, we reduce τ * a from 1.0 to 0.2 while keeping r 0 fixed -compare the second and third column in Fig. 2. As expected, the total backscattered energy drops by almost a factor of five (note the Fig. 2. Backscattered energy as a function of altitude (top row), fraction of regular (second row) and irregular multiple scattering (third row), and multiple scattering correction factor η as a function of aerosol optical depth τ a . The columns pertain to three different aerosol scenarios as indicated in the column heading. different x-axis scales in the top center and top right panels). Naively, one may expect that the single-and multiple-scattering contributions would drop by the same amount, because we have not changed the aerosol phase function. However, the regular multiple-scattering contribution has been reduced more strongly than the single-scattering signal. The latter lies now close to the total backscattered signal (compare the green and the black curve in the top right diagram). To understand this effect, we need to recall that the angular distribution of the scattered radiation is determined by the total phase function, which is a weighted average of the phase functions of aerosols and molecules. As the aerosol optical depth is reduced, the scattering properties of the medium become more strongly influenced by Rayleigh-scattering gas molecules. Thus, the lower τ * a , the weaker the diffraction peak of the total phase function becomes, resulting in more side-scattering, which results in a higher frequency of optical paths out of the FOV of the detector. This results in a lower frequency of regular multiple scattering events within the detector's FOV. In summary, the results confirm our expectations, namely, that regular multiple scattering contributes significantly to the detected signal for moderate to large optical depths through media containing large particles with strongly forward-peaked phase functions.

Multiple scattering correction factor
The bottom row in Fig. 2 shows the multiple-scattering correction factor η as a function of the aerosol optical depth τ a . For the three scenarios shown here, an aerosol optical depth of τ a = 0 corresponds to altitudes z ≥ 3 km, while τ a = τ * a corresponds to z = 0 km. The triangles in the diagrams represent the values of η computed from the GRIMALDI results by use of Eq. (13). In general, η first decreases sharply with τ a , followed by a moderate increase. This can be understood by inspecting Eq. 13. The term ln E ss (r) E(r) decreases with increasing optical depth τ. If this decrease is stronger than the increase in τ, then η decreases with τ. However, if τ increases more rapidly than ln E ss (r) E(r) decreases, then η increases. Since E ms, irr /E 1, we have E ss (r) E(r) ≈ 1 − E ms, reg /E (see Eq. 11). Thus, E ss /E is the mirror image of E ms, reg /E (second row in Fig. 2). Hence, as we descent into the aerosol layer from above, E ss (r) E(r) first decreases sharply with decreasing altitude, followed by a much slower decrease. As a result, η first decreases sharply with increasing τ a , followed by a mild increase.
Note that η scatters much more in the third column than in the second one. This is because for a lower aerosol optical depth we get fewer extinction events, thus poorer statistics in the Monte Carlo approach. As we will see shortly, in all considered cases this computational uncertainty is much lower than the uncertainty due to different aerosol scenarios.
We can model η as a function of τ a with the following ansatz.
The least-squares fits are represented by the black lines in Fig. 2. Here, each set of data points has been fitted separately. We will see shortly that we can use a similar ansatz to fit η for all 48 aerosol scenarios simultaneously. Figure 3 presents η as a function of τ a and r eff obtained from the 48 different aerosol scenarios. The upper left diagram shows results for r eff =0.4 µm (r 0 =0.2 µm). Different symbols represent scenarios with different total aerosol optical depths τ * a , as indicated in the legend. All scenarios in which the aerosol layer extends from the surface up to z max = 3 km are represented by blue symbols; those in which the aerosol layer extends up to z max = 6 km are shown in red.

Dependence of η on τ * a
For any given value of τ a and z max , η is lower for higher values of τ * a . As mentioned in Sect. 3.2, this effect can be explained by the degree of anisotropy of the phase function. As an example, compare the blue symbols (z max = 3 km) for τ * a =2 (pluses) and τ * a =3 (diamonds) for a given optical depth, say, τ a =1. For the optically thicker aerosol layer (τ * a =3), an optical depth of τ a =1 corresponds to a shorter vertical path ∆z than for the optically thinner layer (τ * a =2). Thus, in the optically thicker case, the aerosols are mixed with a smaller amount of air. As a result, the phase function of the aerosol/air mixture is less influenced by the isotropically scattering air molecules. Hence, the phase function has a more pronounced diffraction peak, which gives rise to a higher probablility of regular multiple scattering in the telescope's FOV. Consequently, the multiple-scattering correction factor η deviates more strongly from unity in the optically thicker case (diamonds) than in the optically thinner case (pluses).

Dependence of η on z max
In a similar way, we can understand the differences in η between aerosol layers of different vertical extent. We take, again, a specific example and focus on the results obtained for r eff =0.4 µm and τ * a =3 (upper left, diamonds). But now we compare the blue and the red symbols. In these two cases, the aerosol loads in the vertical column are identical, but they are differently distributed in the vertical direction. In one case (blue), the aerosols are concentrated in a layer up to z max = 3 km, in the other (red) the same amount of aerosols is mixed with a larger volume of air up to z max = 6 km. Again, the larger the amount of air mixed with the aerosols, the less pronounced the diffraction peak of the phase function becomes, resulting in more side-scattering out of the FOV of the detector. Hence for z max = 6 km (red) η deviates less from unity than for z max = 3 km.

Parameterization of η as a function of τ a and r eff
Although we see a systematic effect of z max and τ * a on η, these effects are minor compared to the dependence of η on the aerosol optical depth through the medium τ a and, most importantly, on the effective radius r eff . For this reason, we propose to model η as a function of τ a and r eff , and to treat all other factors as a source of error. We make the following ansatz to model η: where r eff is given in µm. The non-linear least squares fit is represented by the black lines in the four top diagrams in Fig. 3. The fitting program yields the following values for the free parameters: The coefficients are given with 95 % confidence bounds. The bottom panel in Fig. 3 shows a 3D diagram that summarizes all data points for η. The fitted curve is represented by a surface plot. It is evident that the empirical formula given in Eq. (15) models the data reasonably well. Figure (3) shows a considerable variance among different aerosol scenarios. In an inverse model, it is important to quantify the error we introduce by making use of the fitting expression given in Eq. (15). The left diagram of Fig. 4 shows a surface plot of η(r eff , τ a ) based on Eq. (15); the right diagram shows the standard deviation

Error estimate
Here η MC i denote the multiple scattering corrections factors obtained from the Monte Carlo simulations, η i denote the corresponding values computed with Eq. (15), and N denotes the total number of data points. With the exception of those cases for which r eff > 5 µm and τ a < 0.5, the standard deviation is rather uniform. Thus, for the sake of simplicity, we can average ∆η over two subsets of the total ensemble of data points. We obtain

Discussion
Equations (15) and (16) imply lim Hence, unless we have simultaneously high amounts of coarse particles and moderately to high aerosol loads, we can simply assume a constant value of η = 0.93. How common are the cases in which deviations from this asymptotic value are significant? Clearly, episodic events, such as dust storms and volcanic eruptions, can give rise to such situations. However, even in the absence of such events there can be episodes of high concentrations of coarse sea salt aerosols over the oceans. Figure 5 shows can give rise to high aerosol loads; but these aerosols are usually in the accumulation mode, i.e., they have small values of r eff .
To illustrate the significance of multiple scattering for episodic events of high sea salt concentrations, we select an aerosol profile from the results shown in Fig. 5 over the Nordic Sea, which is indicated in the figure by a black 'X'. The coordinates of this location are 67 • N, 4.5 • E. Figure 6 (left) shows the attenuated backscattering coefficient computed for this aerosol profile by use of the MATCH aerosol optics model [45]. The results have been computed for three different choices of η in the transmission term. The black line shows results obtained by using Eq. (15) to model η for each vertical layer. The shaded region indicates the range of uncertainty, which is estimated by varying η between η min and η max , where η min = max(η −∆η, 0), η max = min(η + ∆η, 1), and ∆η is given by (18). We consider the black curve to be the reference case. The green line shows corresponding results obtained by assuming a constant value of η = 0.93. Finally, the red line has been obtained by neglecting multiple scattering entirely (η = 1). We consider the green and the red lines to represent more or less crude approximations. The right panel presents the magnitude of the relative biases in these two approximations in comparison to the reference case, as well as the relative random error in η for the reference case (black dashed line).
The left diagram illustrates that setting η = 0.93 or η = 1 leads to negative biased in β att of a magnitude that exceeds the random error in the reference case. The right diagram reveals that the assumption of η = 0.93 leads to a relative bias up to 17 %. Taking η = 1 results in a relative bias up to 27 %. By contrast, the random error due to the uncertainty in the parameterization given in (15) only goes up to 12 %. In either approximation we obtain a negative bias for the attenuated backscattering coefficient, because multiple scattering acts as a source term in the radiative transfer equation, which enhances the backscattered energy. The standard deviation ∆η is a measure of the random error introduced by using the parameterisation of η given in Eq. (15), rather than solving the radiative transfer problem exactly. We note that this random error may even be further reduced by slightly modifying the approach followed here. A contributing factor to the error was the varying amount of air that was mixed with the aerosols in the different scenarios. For this reason, one could replace Eq. (3) with thus correcting separately for multiple scattering of aerosols and molecules with corresponding factors η a and η m , respectively. It can be expected that η m would be close to 1, and the standard deviation ∆η a of the aerosol multiple scattering correction factor would be smaller than the corresponding value ∆η we estimated in our approach. However, the main conclusions of our study would not be affected.

Summary and Conclusions
We performed Monte-Carlo radiative transfer simulations of lidar return signals at a wavelength of 532 nm in the presence of aerosols. The specifications of the CALIOP instrument on-board CALIPSO have been used. The Monte-Carlo program kept track of how much of the returned energy is due to single scattering, regular multiple scattering, and irregular multiple scattering. Computations were performed for 48 different aerosol scenarios with varying aerosol size, optical depths, and vertical distributions. It was found that (i) irregular multiple scattering can be neglected; its fractional contribution to the returned energy was always less than 0.005. (ii) The contribution of regular multiple scattering to the detected signal mostly depends on the effective radius of the aerosols and the optical depth. To a lesser extent, it depends on the vertical distribution and on the total optical depth of the aerosol layer. We proposed a parameterization of the multiple-scattering correction factor η as a function of effective radius and optical depth. The dependence on other factors, such as the vertical distribution of the aerosols or the total optical depth, have been treated as a random error, which we quantified. If the optical depth is low, or if the effective radius of the aerosols is not too large, one can use an asymptotic value of η = 0.93. However, we demonstrated that using this asymptotic value outside its range of validity can introduce a negative bias in simulations of lidar return signals. More specifically, we considered an episode of high sea salt concentrations over the ocean. Simulations of the attenuated backscattering coefficient revealed that neglecting multiple scattering, or using a constant value of η = 0.93, resulted in a negative bias, the magnitude of which exceeded the random error in our parameterisation. Although our study covered a significant range of different aerosol scenarios, there are many questions left unanswered. We have not considered multiple aerosol layers, aerosol layers above clouds, nonspherical aerosols, or aerosols with significantly different refractive indices (e.g. pure black carbon). Further, our Monte-Carlo radiative transfer code, in its present form, is based on the scalar approximation; polarization effects are neglected. To what extent polarization effects could modulate the findings of this study can only be answered by use of a fully fledged Monte-Carlo code with polarization. Finally, the investigation was limited to a single wavelength. When considering other wavelengths, there will be two competing effects that may impact the significance of multiple scattering. For example, in the UV the aerosol phase function will have a more pronounced diffraction peak (owing to the higher size parameter), which would enhance the contribution of regular multiple scattering. On the other hand, the contribution of air molecules to the total phase function will increase (owing to the ∼ λ −4 dependency of the Rayleigh scattering cross section); this will reduce the contribution of regular multiple scattering. At infrared wavelengths, these effects will be reversed. It remains an open question which effect will dominate under different conditions.