1. 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 $\Delta 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.

 

Fig. 1. Optical paths for (a) single scattering events, (b) regular multiple scattering, and (c) irregular multiple scattering. The physical length of all three optical paths is $2R\pm \Delta R$, where the resolution of the instrument $\Delta R$ is indicated by the shaded area.

Download Full Size | PDF

We further may have to consider optical paths as in diagram c. This illustrates a multiple-scattering light path within the same range of pathlengths $2R\pm \Delta 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\pm \Delta R$, and for which the lowest altitude along the light path lies in the range $z\pm \Delta 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+\Delta 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–8] and experimental studies [9–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–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–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–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 $\cal H$ maps a vector of model variables $\mathbf {x}$ to observable parameters ${\cal H}(\mathbf {x})$. (The vector $\mathbf {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

Platt [9] proposed an empirical correction to account for multiple scattering by reducing the optical depth by a factor $\eta$. This factor lies usually between 0.5 and 1. Thus, the transmission term becomes

In this paper, we perform MC simulations for different aerosol scenarios. We model the multiple scattering correction factor $\eta$ for CALIOP observations at a wavelength of 532 nm under different conditions. In previous publications $\eta$ 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.

2. Methods

2.1 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 $\tau _d$ (integrated extinction along the path of the light ray) without scattering or absorption. This leads to:

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 $\mu$rad; (iv) beam divergence 100 $\mu$rad; (v) vertical resolution 30 m.

2.2 Aerosol scenarios

We compute aerosol optical properties for a wavelength of $\lambda =$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.003$i 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.0008$i [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.79$i [38], an internal mixture of toluene containing 4 % black carbon would have an effective refractive index of $m=1.56+0.003$i (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

Table 1. Mean radius $r_0$, geometric standard deviation $\sigma _g$, and effective radius $r_{\mathrm {eff}}$ for the lognormal size distributions considered in our aerosol scenarios.

View Table

For each size distribution, we consider six different values of aerosol optical depth $\tau _a^\ast$ in the troposphere, namely, $\tau _a^\ast =$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 $\tau _a^\ast$ 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 $\tau _a(z)$ we denote the vertical aerosol optical depth from TOA to altitude $z$. Thus, $\tau _a^\ast =\tau _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.

2.3 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 $\Delta R$. Then the energy that is returned to the detector from points in the interval $R - \Delta R/2 \le R \le R + \Delta R/2$ is given by

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 $\eta$ 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 single-scattering paths (Fig. 1a), from regular multiple-scattering paths (Fig. 1b), and from irregular multiple-scattering paths (Fig. 1c). Let us denote these contributions by

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$ $\mu$m, $\tau _a^\ast =$1.0 (left column), $r_0=2.0$ $\mu$m, $\tau _a^\ast =$1.0 (second column), and $r_0=2.0$ $\mu$m, $\tau _a^\ast =$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_{\mathrm {MS,\;reg}}/E$, and of irregular multiple scattering, $E_{\mathrm {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).

 

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 $\eta$ as a function of aerosol optical depth $\tau _a$. The columns pertain to three different aerosol scenarios as indicated in the column heading.

Download Full Size | PDF

3.1 Dependence of multiple scattering on aerosol size

For small particles with $r_0=0.2$ $\mu$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$ $\mu$m while keeping $\tau _a^\ast$ 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_{\mathrm {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.

3.2 Dependence of multiple scattering on total aerosol optical depth

Next, we reduce $\tau _a^\ast$ 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 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 $\tau _a^\ast$, 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.

3.3 Multiple scattering correction factor

The bottom row in Fig. 2 shows the multiple-scattering correction factor $\eta$ as a function of the aerosol optical depth $\tau _a$. For the three scenarios shown here, an aerosol optical depth of $\tau _a=0$ corresponds to altitudes $z \ge 3$ km, while $\tau _a=\tau _a^\ast$ corresponds to $z=0$ km. The triangles in the diagrams represent the values of $\eta$ computed from the GRIMALDI results by use of Eq. (13). In general, $\eta$ first decreases sharply with $\tau _a$, followed by a moderate increase. This can be understood by inspecting Eq. 13. The term $\ln \frac {E_\textrm {ss}(r)}{E(r)}$ decreases with increasing optical depth $\tau$. If this decrease is stronger than the increase in $\tau$, then $\eta$ decreases with $\tau$. However, if $\tau$ increases more rapidly than $\ln \frac {E_\textrm {ss}(r)}{E(r)}$ decreases, then $\eta$ increases. Since $E_{\mathrm {ms,\;irr}}/E \ll 1$, we have $\frac {E_\textrm {ss}(r)}{E(r)}\approx 1 - E_{\mathrm {ms,\;reg}}/E$ (see Eq. 11). Thus, $E_{\mathrm {ss}}/E$ is the mirror image of $E_{\mathrm {ms,\;reg}}/E$ (second row in Fig. 2). Hence, as we descent into the aerosol layer from above, $\frac {E_\textrm {ss}(r)}{E(r)}$ first decreases sharply with decreasing altitude, followed by a much slower decrease. As a result, $\eta$ first decreases sharply with increasing $\tau _a$, followed by a mild increase.

Note that $\eta$ 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 $\eta$ as a function of $\tau _a$ with the following ansatz.

Figure 3 presents $\eta$ as a function of $\tau _a$ and $r_{\mathrm {eff}}$ obtained from the 48 different aerosol scenarios. The upper left diagram shows results for $r_{\mathrm {eff}}=$0.4 $\mu$m ($r_0=$0.2 $\mu$m). Different symbols represent scenarios with different total aerosol optical depths $\tau _a^\ast$, as indicated in the legend. All scenarios in which the aerosol layer extends from the surface up to $z_{\mathrm {max}}=3$ km are represented by blue symbols; those in which the aerosol layer extends up to $z_{\mathrm {max}}=6$ km are shown in red.

 

Fig. 3. Fit of the multiple scattering correction factor $\eta$ as a function of effective radius $r_{\mathrm {eff}}$ and aerosol optical depth $\tau _a$ according to Eq. (15). The bottom panel shows the full 3D plot, while the upper four panels show cross sections for fixed values of $r_{\mathrm {eff}}$. The least-square fitting result is represented by the black lines (upper four panels) and by the surface plot (bottom panel).

Download Full Size | PDF

3.4 Dependence of $\eta$ on $\tau _a^\ast$

For any given value of $\tau _a$ and $z_{\mathrm {max}}$, $\eta$ is lower for higher values of $\tau _a^\ast$. 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_{\mathrm {max}}=3$ km) for $\tau _a^\ast =$2 (pluses) and $\tau _a^\ast =$3 (diamonds) for a given optical depth, say, $\tau _a=$1. For the optically thicker aerosol layer ($\tau _a^\ast =$3), an optical depth of $\tau _a=$1 corresponds to a shorter vertical path $\Delta z$ than for the optically thinner layer ($\tau _a^\ast =$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 $\eta$ deviates more strongly from unity in the optically thicker case (diamonds) than in the optically thinner case (pluses).

3.5 Dependence of $\eta$ on $z_{\mathrm {max}}$

In a similar way, we can understand the differences in $\eta$ between aerosol layers of different vertical extent. We take, again, a specific example and focus on the results obtained for $r_{\mathrm {eff}}=$0.4 $\mu$m and $\tau _a^\ast =$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_{\mathrm {max}}=3$ km, in the other (red) the same amount of aerosols is mixed with a larger volume of air up to $z_{\mathrm {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_{\mathrm {max}}=6$ km (red) $\eta$ deviates less from unity than for $z_{\mathrm {max}}=3$ km.

3.6 Parameterization of $\eta$ as a function of $\tau _a$ and $r_{\mathrm {eff}}$

Although we see a systematic effect of $z_{\mathrm {max}}$ and $\tau _a^\ast$ on $\eta$, these effects are minor compared to the dependence of $\eta$ on the aerosol optical depth through the medium $\tau _a$ and, most importantly, on the effective radius $r_{\mathrm {eff}}$. For this reason, we propose to model $\eta$ as a function of $\tau _a$ and $r_{\mathrm {eff}}$, and to treat all other factors as a source of error. We make the following ansatz to model $\eta$:

The bottom panel in Fig. 3 shows a 3D diagram that summarizes all data points for $\eta$. 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.

3.7 Error estimate

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 $\eta (r_{\mathrm {eff}},\tau _a)$ based on Eq. (15); the right diagram shows the standard deviation

 

Fig. 4. Multiple scattering correction factor $\eta$ as a function of $r_{\mathrm {eff}}$ and $\tau _a$ according to Eq. (15) (left), and standard deviation $\Delta \eta$ (right).

Download Full Size | PDF

4. Discussion

Equations (15) and (16) imply

 

Fig. 5. Effective radius 750 m above the ground (left) and aerosol optical depth (right) obtained with the MATCH model on June 13 2006, 1800 UTC. The ’X’ in the Nordic Sea marks the location selected for the test case shown in Fig. 6.

Download Full Size | PDF

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$^\circ$ N, 4.5$^\circ$ 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 $\eta$ in the transmission term. The black line shows results obtained by using Eq. (15) to model $\eta$ for each vertical layer. The shaded region indicates the range of uncertainty, which is estimated by varying $\eta$ between $\eta _{\mathrm {min}}$ and $\eta _{\mathrm {max}}$, where $\eta _{\mathrm {min}}=\mathrm {max}(\eta -\Delta \eta ,\;0)$, $\eta _{\mathrm {max}}=\mathrm {min}(\eta +\Delta \eta ,\;1)$, and $\Delta \eta$ 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 $\eta =0.93$. Finally, the red line has been obtained by neglecting multiple scattering entirely ($\eta =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 $\eta$ for the reference case (black dashed line).

 

Fig. 6. Attenuated backscattering coefficient $\beta _{\mathrm {att}}$ simulated with the MATCH observation operator (left). The results have been obtained by (i) modelling $\eta$ according to Eq. (15) (black, uncertainty indicated by shaded region), (ii) by setting $\eta =$0.93 (green), and (iii) by setting $\eta =$1.00 (red). The right diagram shows the magnitude of the relative biases of cases (ii) and (iii), and the relative random error in case (i).

Download Full Size | PDF

The left diagram illustrates that setting $\eta =0.93$ or $\eta =1$ leads to negative biased in $\beta _{\mathrm {att}}$ of a magnitude that exceeds the random error in the reference case. The right diagram reveals that the assumption of $\eta =0.93$ leads to a relative bias up to 17 %. Taking $\eta =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 $\Delta \eta$ is a measure of the random error introduced by using the parameterisation of $\eta$ 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

5. 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 $\eta$ 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 $\eta =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 $\eta =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 $\sim \lambda ^{-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.