# Multiple-scattering effects on single-wavelength lidar sounding of multi-layered clouds

Valery Shcherbakov, Frédéric Szczap, Guillaume Mioche, Céline Cornet

► To cite this version:

Valery Shcherbakov, Frédéric Szczap, Guillaume Mioche, Céline Cornet. Multiple-scattering effects on single-wavelength lidar sounding of multi-layered clouds. Atmospheric Measurement Techniques, 2024, 17 (9), pp.3011-3028. 10.5194/amt-17-3011-2024. hal-04599479

HAL Id: hal-04599479

https://uca.hal.science/hal-04599479

Submitted on 4 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Multiple-scattering effects on single-wavelength lidar sounding of multi-layered clouds

Valery Shcherbakov, Frédéric Szczap, Guillaume Mioche, Céline Cornet

► To cite this version:

Valery Shcherbakov, Frédéric Szczap, Guillaume Mioche, Céline Cornet. Multiple-scattering effects on single-wavelength lidar sounding of multi-layered clouds. Atmospheric Measurement Techniques, 2024, 17 (9), pp.3011-3028. 10.5194/amt-17-3011-2024. hal-04599479

HAL Id: hal-04599479

https://uca.hal.science/hal-04599479v1

Submitted on 4 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Multiple-scattering effects on single-wavelength lidar sounding of multi-layered clouds

Valery Shcherbakov<fnref n="1" />, Frédéric Szczap<fnref n="1" />, Guillaume Mioche<fnref n="1" />,<fnref n="2" />, and Céline Cornet<fnref n="3" />

**Correspondence:** Valery Shcherbakov (v.shcherbakov@opgc.univ-bpclermont.fr)

Received: 30 May 2023 – Discussion started: 21 June 2023

Revised: 14 March 2024 – Accepted: 25 March 2024 – Published: 17 May 2024

**Abstract.** We performed Monte Carlo simulations of single-wavelength lidar signals from multi-layered clouds with special attention focused on the multiple-scattering (MS) effect in regions of the cloud-free molecular atmosphere (i.e. between layers or outside a cloud system). Despite the fact that the strength of lidar signals from the molecular atmosphere is much lower compared to the in-cloud intervals, studies of MS effects in such regions are of interest from scientific and practical points of view.

The MS effect on lidar signals always decreases with the increasing distance from the cloud far edge. The decrease is the direct consequence of the fact that the forward peak of particle phase functions is much larger than the receiver field of view (RFOV). Therefore, the photons scattered within the forward peak escape the sampling volume formed by the RFOV (i.e. the escape effect). We demonstrated that the escape effect is an inherent part of MS properties within the free atmosphere beyond the cloud far edge.

In the cases of the ground-based lidar, the MS contribution is lower than 5 % within the regions of the cloud-free molecular atmosphere with a distance from the cloud far edge of about 1 km or higher. In the cases of the space-borne lidar, the rate of decrease of the MS contribution is so slow that the threshold of 5 % can hardly be reached. In addition, the effect of non-uniform beam filling is extremely strong. Therefore, practitioners should employ, with proper precautions, lidar data from regions below the cloud base when treating data of a space-borne lidar.

In the case of two-layered cloud, the distance of 1 km is sufficiently large so that the scattered photons emerging from

the first layer do not affect signals from the second layer when we are dealing with the ground-based lidar. In contrast, signals from the near edge of the second cloud layer are severely affected by the photons emerging from the first layer in the case of a space-borne lidar.

We evaluated the Eloranta model (EM) in extreme conditions and showed its good performance in the cases of ground-based and space-borne lidars. At the same time, we revealed the shortcoming that can affect practical applications of the EM. Namely, values of the key parameters – i.e. the ratios of phase functions in the backscatter direction for the nth-order-scattered photon and a singly scattered photon – depend not only on the particle phase function but also on the distance from a lidar to the cloud and the receiver field of view. Those ratios vary within a quite large range, and the MS contribution to lidar signals can be largely overestimated or underestimated if erroneous values of the ratios are assigned to the EM.

## 1 Introduction

It is well recognized that multiple scattering (MS) inevitably affects the data of space-borne lidars (see e.g. Winker, 2003; Young and Vaughan, 2009; Shcherbakov et al., 2022). As for ground-based lidars, the MS relative contribution to lidar signals can exceed 20 % even in the case of cirrus clouds with an extinction coefficient of $ε_p = 0.2 \text{ km}^{-1}$ and can reach 200 % when $ε_p = 1.0 \text{ km}^{-1}$ (Shcherbakov et al., 2022). On the other hand, the MS effect is within measurement errors of ground-

<footnote n="1">¹Université Clermont Auvergne, CNRS, UMR 6016, Laboratoire de Météorologie Physique, 63178 Aubière, France</footnote>

<footnote n="2">²Université Clermont Auvergne, Institut Universitaire de Technologie Clermont
Auvergne - site de Montluçon, 03100 Montluçon, France</footnote>

<footnote n="3">³Université Lille, CNRS, UMR 8518 - LOA - Laboratoire d'Optique Atmosphérique, F-59000 Lille, France</footnote>
<figure><img src="image_2.png" /></figure>

# Multiple-scattering effects on single-wavelength lidar sounding of multi-layered clouds

Valery Shcherbakov<fnref n="1" />, Frédéric Szczap<fnref n="1" />, Guillaume Mioche<fnref n="1" />,<fnref n="2" />, and Céline Cornet<fnref n="3" />

**Correspondence:** Valery Shcherbakov (v.shcherbakov@opgc.univ-bpclermont.fr)

Received: 30 May 2023 – Discussion started: 21 June 2023

Revised: 14 March 2024 – Accepted: 25 March 2024 – Published: 17 May 2024

**Abstract.** We performed Monte Carlo simulations of single-wavelength lidar signals from multi-layered clouds with special attention focused on the multiple-scattering (MS) effect in regions of the cloud-free molecular atmosphere (i.e. between layers or outside a cloud system). Despite the fact that the strength of lidar signals from the molecular atmosphere is much lower compared to the in-cloud intervals, studies of MS effects in such regions are of interest from scientific and practical points of view.

The MS effect on lidar signals always decreases with the increasing distance from the cloud far edge. The decrease is the direct consequence of the fact that the forward peak of particle phase functions is much larger than the receiver field of view (RFOV). Therefore, the photons scattered within the forward peak escape the sampling volume formed by the RFOV (i.e. the escape effect). We demonstrated that the escape effect is an inherent part of MS properties within the free atmosphere beyond the cloud far edge.

In the cases of the ground-based lidar, the MS contribution is lower than 5 % within the regions of the cloud-free molecular atmosphere with a distance from the cloud far edge of about 1 km or higher. In the cases of the space-borne lidar, the rate of decrease of the MS contribution is so slow that the threshold of 5 % can hardly be reached. In addition, the effect of non-uniform beam filling is extremely strong. Therefore, practitioners should employ, with proper precautions, lidar data from regions below the cloud base when treating data of a space-borne lidar.

In the case of two-layered cloud, the distance of 1 km is sufficiently large so that the scattered photons emerging from

the first layer do not affect signals from the second layer when we are dealing with the ground-based lidar. In contrast, signals from the near edge of the second cloud layer are severely affected by the photons emerging from the first layer in the case of a space-borne lidar.

We evaluated the Eloranta model (EM) in extreme conditions and showed its good performance in the cases of ground-based and space-borne lidars. At the same time, we revealed the shortcoming that can affect practical applications of the EM. Namely, values of the key parameters – i.e. the ratios of phase functions in the backscatter direction for the nth-order-scattered photon and a singly scattered photon – depend not only on the particle phase function but also on the distance from a lidar to the cloud and the receiver field of view. Those ratios vary within a quite large range, and the MS contribution to lidar signals can be largely overestimated or underestimated if erroneous values of the ratios are assigned to the EM.

## 1 Introduction

It is well recognized that multiple scattering (MS) inevitably affects the data of space-borne lidars (see e.g. Winker, 2003; Young and Vaughan, 2009; Shcherbakov et al., 2022). As for ground-based lidars, the MS relative contribution to lidar signals can exceed 20 % even in the case of cirrus clouds with an extinction coefficient of $ε_p = 0.2 \text{ km}^{-1}$ and can reach 200 % when $ε_p = 1.0 \text{ km}^{-1}$ (Shcherbakov et al., 2022). On the other hand, the MS effect is within measurement errors of ground-

<footnote n="1">¹Université Clermont Auvergne, CNRS, UMR 6016, Laboratoire de Météorologie Physique, 63178 Aubière, France</footnote>

<footnote n="2">²Université Clermont Auvergne, Institut Universitaire de Technologie Clermont
Auvergne - site de Montluçon, 03100 Montluçon, France</footnote>

<footnote n="3">³Université Lille, CNRS, UMR 8518 - LOA - Laboratoire d'Optique Atmosphérique, F-59000 Lille, France</footnote>
based and airborne lidars – i.e. the single-scattering (SS) ap-
proximation – can be used when the receiver field of view
(RFOV) is quite narrow and/or the extinction coefficient is
quite low (Shcherbakov et al., 2022).

Lidar signals are used as input data to retrieve profiles
of the characteristics of interest – e.g. the extinction coeffi-
cient, the backscatter coefficient, and the lidar ratio (see e.g.
Young and Vaughan, 2009). If the MS contribution cannot
be neglected, a solution to the corresponding inverse prob-
lem has to take into account the MS effect in order to avoid
biased retrievals. To put it differently, a correction for the
effects of multiple scattering has to be applied (Young and
Vaughan, 2009).

Generally, any algorithm to solve an inverse problem is
based on a solution to the corresponding direct, i.e. forward
problem (see e.g. Rodgers, 2000; Neto and Neto, 2013). It
is obvious that retrieval quality is closely related to the ac-
curacy of the direct modelling. To put it in relation to lidar
sounding, the quality of the correction for MS effects di-
rectly depends on the accuracy of lidar signal modelling in
multiple-scattering conditions.

Expressed mathematically, any modelling in MS condi-
tions has to be based on the radiative transfer equation and
corresponding boundary conditions (see e.g. Marchuk et al.,
2013). It is well known that there exist only a few cases when
a solution to the radiative transfer equation can be obtained
as an analytical expression. Thus, numerical methods – e.g.
Monte Carlo (MC) simulations (Marchuk et al., 2013) or ap-
proximate models – are largely applied to obtain solutions.
Referring to lidar sounding, a good review of developed ap-
proximate methods and models can be found in Bissonnette
(2005). Using some cases as examples, good performance of
approximate models was underlined by their authors. At the
same time, we believe that the accuracy level and the appli-
cability bounds of approximate models still need to be rigor-
ously evaluated.

What was suggested in the work by Weinman (1968) was approximating a phase function of cloud particles by a somewhat smoothly varying function of the scattering angle plus a narrow Gaussian peak for small angles. That idea was used to develop models that simulate lidar signals in multiple-scattering conditions, and promising results were obtained (see e.g. Eloranta, 1998; Hogan, 2006, 2008). For example, good agreement with MC simulations of the second, third, and fourth order of scattering was shown for the case of an extinction coefficient of 16.7 km⁻¹ (Eloranta, 1998). Moreover, the models (Eloranta, 1998; Hogan, 2006, 2008) were employed by practitioners to account for MS effects while retrieving clouds’ optical properties on the basis of experimental lidar data (see e.g. Nakoudi et al., 2021; Seifert et al., 2007; Delanoë and Hogan, 2008).

On the other hand, results of MC simulations published
in the literature (see e.g. Flesia and Starkov, 1996; Donovan,
2016; Reverdy et al., 2015; Szczap et al., 2021) evidenced
the following. As it is expected, lidar signals from regions

of the cloud-free molecular atmosphere – namely between
cloud layers and/or beyond the far edge of a cloud system –
are affected by the scattered light emerging from clouds. To
our knowledge, optical processes and cloud characteristics
that govern the effect on lidar signals of the emerging light
(i.e. the MS effect) and its distinctive features have not been
addressed in detail in the literature.

The main objective of this work is to perform Monte Carlo
simulations of single-wavelength lidar signals from multi-
layered clouds with special attention focused on the pecu-
liarities of the MS effect in regions of the cloud-free molec-
ular atmosphere – i.e. between layers or outside a cloud sys-
tem. In our opinion, such a study helps to obtain further in-
sight into the problem that is MS effects. At the same time, it
has practical applications because some inverse-problem al-
gorithms use lidar data taken from the range of the cloud-free
atmosphere beyond the far edge of a cloud – for example, the
two-way transmittance method (see e.g. Young and Vaughan,
2009; Giannakaki et al., 2007; and references therein).

The second objective of this work is to evaluate the perfor-
mance of an approximate model with the focus on the cloud-
free regions. We have chosen for the evaluation the model
developed by Eloranta (1998) because multiple integrals are
in its core (see Eq. 11 in Eloranta, 1998, and Eqs. 3–4 in Elo-
ranta, 2000). A code to compute multiple integrals belongs
to the domain of basic programming tasks. Therefore, it is an
easy matter to develop a code corresponding to the Eloranta
model.

Throughout this work, the majority of the results are
shown and discussed in terms of relative contributions of
multiple scattering (see definitions in Sect. 2.1) – that is,
multiple-to-single-scattering lidar return ratios (Bissonnette
et al., 1995). At the same time, the reader should keep in
mind that the ratios do not show lidar signals but the rela-
tive contribution of MS with respect to lidar signals obtained
under the SS approximation. From the point of view of the
inverse problem, the ratios are of key importance because
they show the level of retrieval errors if the SS approxima-
tion were applied in an inversion algorithm. If different mod-
els, which simulate the direct problem, are compared, the
ratios clearly evidence the differences between the results
(see e.g. Bissonnette et al., 1995). In our opinion, a com-
parison between lidar signals instead of between ratios suf-
fers from a severe disadvantage for the following reasons.
Even when lidar signals are corrected for the offset and in-
strumental factors, they do decrease exponentially because
of the extinction. Consequently, the semi-logarithmic plot is
usually employed to show data, and even important differ-
ences are hardly noticeable in such figures. Moreover, when
differences are seen, it is difficult for the reader to estimate
their importance from the point of view of the inverse prob-
lem.

The MS effect on lidar signals depends on a number of pa-
rameters – namely, the configuration (the distance to a cloud;
the emitter field of view, EFOV; and the RFOV) and the cloud
based and airborne lidars – i.e. the single-scattering (SS) ap-
proximation – can be used when the receiver field of view
(RFOV) is quite narrow and/or the extinction coefficient is
quite low (Shcherbakov et al., 2022).

Lidar signals are used as input data to retrieve profiles
of the characteristics of interest – e.g. the extinction coeffi-
cient, the backscatter coefficient, and the lidar ratio (see e.g.
Young and Vaughan, 2009). If the MS contribution cannot
be neglected, a solution to the corresponding inverse prob-
lem has to take into account the MS effect in order to avoid
biased retrievals. To put it differently, a correction for the
effects of multiple scattering has to be applied (Young and
Vaughan, 2009).

Generally, any algorithm to solve an inverse problem is
based on a solution to the corresponding direct, i.e. forward
problem (see e.g. Rodgers, 2000; Neto and Neto, 2013). It
is obvious that retrieval quality is closely related to the ac-
curacy of the direct modelling. To put it in relation to lidar
sounding, the quality of the correction for MS effects di-
rectly depends on the accuracy of lidar signal modelling in
multiple-scattering conditions.

Expressed mathematically, any modelling in MS condi-
tions has to be based on the radiative transfer equation and
corresponding boundary conditions (see e.g. Marchuk et al.,
2013). It is well known that there exist only a few cases when
a solution to the radiative transfer equation can be obtained
as an analytical expression. Thus, numerical methods – e.g.
Monte Carlo (MC) simulations (Marchuk et al., 2013) or ap-
proximate models – are largely applied to obtain solutions.
Referring to lidar sounding, a good review of developed ap-
proximate methods and models can be found in Bissonnette
(2005). Using some cases as examples, good performance of
approximate models was underlined by their authors. At the
same time, we believe that the accuracy level and the appli-
cability bounds of approximate models still need to be rigor-
ously evaluated.

What was suggested in the work by Weinman (1968) was approximating a phase function of cloud particles by a somewhat smoothly varying function of the scattering angle plus a narrow Gaussian peak for small angles. That idea was used to develop models that simulate lidar signals in multiple-scattering conditions, and promising results were obtained (see e.g. Eloranta, 1998; Hogan, 2006, 2008). For example, good agreement with MC simulations of the second, third, and fourth order of scattering was shown for the case of an extinction coefficient of 16.7 km⁻¹ (Eloranta, 1998). Moreover, the models (Eloranta, 1998; Hogan, 2006, 2008) were employed by practitioners to account for MS effects while retrieving clouds’ optical properties on the basis of experimental lidar data (see e.g. Nakoudi et al., 2021; Seifert et al., 2007; Delanoë and Hogan, 2008).

On the other hand, results of MC simulations published
in the literature (see e.g. Flesia and Starkov, 1996; Donovan,
2016; Reverdy et al., 2015; Szczap et al., 2021) evidenced
the following. As it is expected, lidar signals from regions

of the cloud-free molecular atmosphere – namely between
cloud layers and/or beyond the far edge of a cloud system –
are affected by the scattered light emerging from clouds. To
our knowledge, optical processes and cloud characteristics
that govern the effect on lidar signals of the emerging light
(i.e. the MS effect) and its distinctive features have not been
addressed in detail in the literature.

The main objective of this work is to perform Monte Carlo
simulations of single-wavelength lidar signals from multi-
layered clouds with special attention focused on the pecu-
liarities of the MS effect in regions of the cloud-free molec-
ular atmosphere – i.e. between layers or outside a cloud sys-
tem. In our opinion, such a study helps to obtain further in-
sight into the problem that is MS effects. At the same time, it
has practical applications because some inverse-problem al-
gorithms use lidar data taken from the range of the cloud-free
atmosphere beyond the far edge of a cloud – for example, the
two-way transmittance method (see e.g. Young and Vaughan,
2009; Giannakaki et al., 2007; and references therein).

The second objective of this work is to evaluate the perfor-
mance of an approximate model with the focus on the cloud-
free regions. We have chosen for the evaluation the model
developed by Eloranta (1998) because multiple integrals are
in its core (see Eq. 11 in Eloranta, 1998, and Eqs. 3–4 in Elo-
ranta, 2000). A code to compute multiple integrals belongs
to the domain of basic programming tasks. Therefore, it is an
easy matter to develop a code corresponding to the Eloranta
model.

Throughout this work, the majority of the results are
shown and discussed in terms of relative contributions of
multiple scattering (see definitions in Sect. 2.1) – that is,
multiple-to-single-scattering lidar return ratios (Bissonnette
et al., 1995). At the same time, the reader should keep in
mind that the ratios do not show lidar signals but the rela-
tive contribution of MS with respect to lidar signals obtained
under the SS approximation. From the point of view of the
inverse problem, the ratios are of key importance because
they show the level of retrieval errors if the SS approxima-
tion were applied in an inversion algorithm. If different mod-
els, which simulate the direct problem, are compared, the
ratios clearly evidence the differences between the results
(see e.g. Bissonnette et al., 1995). In our opinion, a com-
parison between lidar signals instead of between ratios suf-
fers from a severe disadvantage for the following reasons.
Even when lidar signals are corrected for the offset and in-
strumental factors, they do decrease exponentially because
of the extinction. Consequently, the semi-logarithmic plot is
usually employed to show data, and even important differ-
ences are hardly noticeable in such figures. Moreover, when
differences are seen, it is difficult for the reader to estimate
their importance from the point of view of the inverse prob-
lem.

The MS effect on lidar signals depends on a number of pa-
rameters – namely, the configuration (the distance to a cloud;
the emitter field of view, EFOV; and the RFOV) and the cloud
optical characteristics (the extinction coefficient, the albedo, the scattering matrix, and so on) (see e.g. Shcherbakov et al., 2022, and references therein). Because the MC method is very time-consuming, it is not suited for taking into account variations in all parameters mentioned. Therefore, our study was restricted to the cases when all cloud layers were within the range between the altitudes of 8 and 11 km. Almost all MC simulations were performed for cloud particles with an extinction coefficient of $ε_p(h) = 1.0 \text{ km}^{-1}$ for the following reasons. On the one hand, technical capacities of contemporary lidars provide us with the possibility of recording signals from the cloud-free atmosphere beyond the far edge of a cloud with an optical thickness of $τ_p = 3.0$. On the other hand, the MS effect cannot be neglected and is clearly seen in a number of cases (Shcherbakov et al., 2022). Our choice of the parameter values was deliberate despite the fact that the altitude range $H \in [8, 11] \text{ km}$ does not correspond to the usual altitudes of warm clouds and that the value $ε_p(h) = 1.0 \text{ km}^{-1}$ is quite small for water clouds and rather high for cirrus clouds. With such a choice, the phase-function impact on multiple scattering is free of the interference of other parameter variations.

For the sake of brevity, the term “sampling volume” is used in this work in reference to the volume bounded by the RFOV of a lidar in 3D space.

Section 2 addresses the methodology of our Monte Carlo simulations in detail. Sections 3 and 4 show our simulation results for a ground-based and a space-borne lidar, respectively. Section 5 is devoted to conclusions. Appendix A is focused on the Eloranta model and its input parameters in homogeneous-cloud conditions; Appendix B is devoted to the simulation uncertainty analysis.

# 2 Methodology

## 2.1 Background

We use the following notations in this work. The function $S_1(h)$ characterizes lidar signals under the SS approximation (corrected for the offset, the instrumental factors, and the two-way ozone transmittance):

$$S_1(h) = [\beta_p(h) + \beta_m(h)] \cdot T^2(h) = [\beta_p(h) + \beta_m(h)] \cdot T_m^2(h) \cdot T_p^2(h), \quad (1)$$

where $h$ is the distance from the lidar; $β_p(h)$ and $β_m(h)$ represent the backscatter contributions from particles and from the atmospheric molecules; $T^2(h) = T_m^2(h) \cdot T_p^2(h)$ is the two-way transmittance from the lidar to the range, $h$; and $T_m^2(h)$ and $T_p^2(h)$ are the molecular and the particulate transmittances, respectively. $T_p^2(h) = 1$ if $h \le h_b$, where $h_b$ is the distance to the cloud near edge; $T_p^2(h) = \exp[-2τ_p(h_b, h)]$ when $h \ge h_b$, where $τ_p(h_b, h) = \int_{h_b}^h ε_p(h') dh'$ is the cloud

optical depth and $ε_p(h)$ is the extinction coefficient of particles.

The notation $S_k(h), k = 2, 3, \dots, n$, represents lidar signals (corrected for the offset, the instrumental factors, and the two-way ozone transmittance) when all scattering events from the first (single scattering) up to the $k$th inclusive are taken into account; $(k+1)$th and higher orders of scattering are neglected. For example, $S_2(h)$ is the double-scattering (DS) approximation (see e.g. Bissonnette, 2005). The notation $S_{MS}(h)$ means that all scattering events are taken into consideration. Of course, only $S_{MS}(h)$ can be obtained from real lidar measurements. $S_k(h)$ carries useful data of simulations computed with the Monte Carlo method or an approximate model.

Below, we use the following characteristics to compare the simulations results:

$$R_{k\text{tol}}(h) = [S_k(h) - S_{k-1}(h)]/S_1(h), \quad (2)$$

$$R_{M\text{Stol}}(h) = [S_{MS}(h) - S_1(h)]/S_1(h). \quad (3)$$

The ratio $R_{k\text{tol}}(h)$ is the relative contribution of the $k$th order of scattering to a lidar signal. For example, $R_{k\text{tol}}(h) = [S_4(h) - S_3(h)]/S_1(h)$ provides the relative contribution of the fourth order of scattering. The ratio $R_{M\text{Stol}}(h)$ is the relative contribution of multiple scattering to a lidar signal; i.e. all orders of scattering were taken into consideration. It is evident that

$$R_{M\text{Stol}}(h) = \sum_{k=2}^{k=n} R_{k\text{tol}}(h). \quad (4)$$

The majority of the results of this work are shown and discussed below in terms of ratios $R_{M\text{Stol}}(h)$ and $R_{2\text{tol}}(h)$.

In order to discuss our results in terms used in the literature (see e.g. Young and Vaughan, 2009; Young et al., 2013; Vaughan et al., 2009; Garnier et al., 2015; and references therein), we employ the following notations and relationships. The “attenuated backscatter”, i.e. lidar signals $S_{MS}(h)$, computed in multiple-scattering conditions (corrected for the offset, the instrumental factors, and the two-way ozone transmittance), is expressed as (see e.g. Young et al., 2013)

$$S_{MS}(h) = [\beta_p(h) + \beta_m(h)] \cdot T_m^2(h) \cdot T_{\text{PA}}^2(h), \quad (5)$$

where $T_{\text{PA}}^2(h) = \exp[-2η(h_b, h)τ_p(h_b, h)]$ is the apparent particulate two-way transmittance; $η(h_b, h)$ is the multiple-scattering function (see Appendix in Shcherbakov et al., 2022, for details) – i.e. a parameterization describing the effect of MS on particulate extinction (Young et al., 2013).

$$S_m(h) = [\beta_m(h)] \cdot T_m^2(h) \quad (6)$$

is a lidar signal computed under the condition of the free molecular atmosphere using available meteorological data (see e.g. Winker et al., 2009).

The scattering ratio, $R(h) = [\beta_p(h) + \beta_m(h)]/\beta_m(h) = [\beta_p(h)/\beta_m(h)] + 1$, is a useful parameter to work with lidar
optical characteristics (the extinction coefficient, the albedo, the scattering matrix, and so on) (see e.g. Shcherbakov et al., 2022, and references therein). Because the MC method is very time-consuming, it is not suited for taking into account variations in all parameters mentioned. Therefore, our study was restricted to the cases when all cloud layers were within the range between the altitudes of 8 and 11 km. Almost all MC simulations were performed for cloud particles with an extinction coefficient of $\epsilon_p(h) = 1.0 \text{ km}^{-1}$ for the following reasons. On the one hand, technical capacities of contemporary lidars provide us with the possibility of recording signals from the cloud-free atmosphere beyond the far edge of a cloud with an optical thickness of $\tau_p = 3.0$. On the other hand, the MS effect cannot be neglected and is clearly seen in a number of cases (Shcherbakov et al., 2022). Our choice of the parameter values was deliberate despite the fact that the altitude range $H \in [8, 11] \text{ km}$ does not correspond to the usual altitudes of warm clouds and that the value $\epsilon_p(h) = 1.0 \text{ km}^{-1}$ is quite small for water clouds and rather high for cirrus clouds. With such a choice, the phase-function impact on multiple scattering is free of the interference of other parameter variations.

For the sake of brevity, the term “sampling volume” is used in this work in reference to the volume bounded by the RFOV of a lidar in 3D space.

Section 2 addresses the methodology of our Monte Carlo simulations in detail. Sections 3 and 4 show our simulation results for a ground-based and a space-borne lidar, respectively. Section 5 is devoted to conclusions. Appendix A is focused on the Eloranta model and its input parameters in homogeneous-cloud conditions; Appendix B is devoted to the simulation uncertainty analysis.

# 2 Methodology

## 2.1 Background

We use the following notations in this work. The function $S_1(h)$ characterizes lidar signals under the SS approximation (corrected for the offset, the instrumental factors, and the two-way ozone transmittance):

$$S_1(h) = [\beta_p(h) + \beta_m(h)] \cdot T^2(h) = [\beta_p(h) + \beta_m(h)] \cdot T_m^2(h) \cdot T_p^2(h), \quad (1)$$

where $h$ is the distance from the lidar; $\beta_p(h)$ and $\beta_m(h)$ represent the backscatter contributions from particles and from the atmospheric molecules; $T^2(h) = T_m^2(h) \cdot T_p^2(h)$ is the two-way transmittance from the lidar to the range, $h$; and $T_m^2(h)$ and $T_p^2(h)$ are the molecular and the particulate transmittances, respectively. $T_p^2(h) = 1$ if $h \le h_b$, where $h_b$ is the distance to the cloud near edge; $T_p^2(h) = \exp[-2\tau_p(h_b, h)]$ when $h \ge h_b$, where $\tau_p(h_b, h) = \int_{h_b}^h \epsilon_p(h') dh'$ is the cloud

optical depth and $\epsilon_p(h)$ is the extinction coefficient of particles.

The notation $S_k(h)$, $k = 2, 3, \dots, n$, represents lidar signals (corrected for the offset, the instrumental factors, and the two-way ozone transmittance) when all scattering events from the first (single scattering) up to the $k$th inclusive are taken into account; $(k+1)$th and higher orders of scattering are neglected. For example, $S_2(h)$ is the double-scattering (DS) approximation (see e.g. Bissonnette, 2005). The notation $S_{MS}(h)$ means that all scattering events are taken into consideration. Of course, only $S_{MS}(h)$ can be obtained from real lidar measurements. $S_k(h)$ carries useful data of simulations computed with the Monte Carlo method or an approximate model.

Below, we use the following characteristics to compare the simulations results:

$$R_{k\text{to}1}(h) = [S_k(h) - S_{k-1}(h)]/S_1(h), \quad (2)$$

$$R_{M\text{Sto}1}(h) = [S_{MS}(h) - S_1(h)]/S_1(h). \quad (3)$$

The ratio $R_{k\text{to}1}(h)$ is the relative contribution of the $k$th order of scattering to a lidar signal. For example, $R_{k\text{to}1}(h) = [S_4(h) - S_3(h)]/S_1(h)$ provides the relative contribution of the fourth order of scattering. The ratio $R_{M\text{Sto}1}(h)$ is the relative contribution of multiple scattering to a lidar signal; i.e. all orders of scattering were taken into consideration. It is evident that

$$R_{M\text{Sto}1}(h) = \sum_{k=2}^{k=n} R_{k\text{to}1}(h). \quad (4)$$

The majority of the results of this work are shown and discussed below in terms of ratios $R_{M\text{Sto}1}(h)$ and $R_{2\text{to}1}(h)$.

In order to discuss our results in terms used in the literature (see e.g. Young and Vaughan, 2009; Young et al., 2013; Vaughan et al., 2009; Garnier et al., 2015; and references therein), we employ the following notations and relationships. The “attenuated backscatter”, i.e. lidar signals $S_{MS}(h)$, computed in multiple-scattering conditions (corrected for the offset, the instrumental factors, and the two-way ozone transmittance), is expressed as (see e.g. Young et al., 2013)

$$S_{MS}(h) = [\beta_p(h) + \beta_m(h)] \cdot T_m^2(h) \cdot T_{\text{pA}}^2(h), \quad (5)$$

where $T_{\text{pA}}^2(h) = \exp[-2\eta(h_b, h)\tau_p(h_b, h)]$ is the apparent particulate two-way transmittance; $\eta(h_b, h)$ is the multiple-scattering function (see Appendix in Shcherbakov et al., 2022, for details) – i.e. a parameterization describing the effect of MS on particulate extinction (Young et al., 2013).

$$S_m(h) = [\beta_m(h)] \cdot T_m^2(h) \quad (6)$$

is a lidar signal computed under the condition of the free molecular atmosphere using available meteorological data (see e.g. Winker et al., 2009).

The scattering ratio, $R(h) = [\beta_p(h) + \beta_m(h)]/\beta_m(h) = [\beta_p(h)/\beta_m(h)] + 1$, is a useful parameter to work with lidar
data; R (h) = 1 for the molecular atmosphere. The attenuated scattering ratio can be defined as (see e.g. Winker et al., 2009)

$$R'_{MS}(h) = S_{MS}(h)/S_{m}(h) = R(h) \cdot T_{pA}^{2}(h). \quad (7)$$

The definitions of a MS lidar signal, S<sub>MS</sub> (h), with Eq. (5); of a SS lidar signal, S<sub>1</sub> (h), with Eq. (1); and of the attenuated scattering ratio with Eq. (7) lead to the following properties for the intervals of the cloud-free molecular atmosphere:

$$R'_{MS}(h) = T_{pA}^{2}, \quad (8)$$

$$S_{1}(h) = [\beta_{m}(h)] \cdot T_{m}^{2}(h) \cdot T_{p}^{2} = S_{m}(h) \cdot T_{p}^{2}, \quad (9)$$

$$S_{MS}(h) = [\beta_{m}(h)] \cdot T_{m}^{2}(h) \cdot T_{pA}^{2} = S_{m}(h) \cdot T_{pA}^{2}, \quad (10)$$

where T<sub>p</sub><sup>2</sup> = T<sub>pA</sub><sup>2</sup> = 1 when h < h<sub>b</sub>; T<sub>p</sub><sup>2</sup> = const and T<sub>pA</sub><sup>2</sup> = const when h > h<sub>end</sub>, where h<sub>end</sub> is the distance from the cloud far edge.

It follows from the relationships above that

$$R_{MS101}(h) = 0, R'_{MS}(h) = 1 \quad (11)$$

when h < h<sub>b</sub> and

$$R_{MS101}(h) = \left( \frac{T_{pA}^2}{T_p^2} \right) - 1 = \text{const, } R'_{MS}(h) = T_{pA}^2 = \text{const} \quad (12)$$

when h > h<sub>end</sub>.

In addition, the following relationship can be useful for the interpretation of Monte Carlo data:

$$R'_{MS}(h) = [R_{MS101}(h) + 1] \cdot T_{p}^{2} = \text{const when } h > h_{\text{end}}. \quad (13)$$

It should be underlined that if measurement errors are neglected, the functions R<sub>MS101</sub> (h) and R'<sub>MS</sub> (h) are expected to be constant when we are dealing within the intervals of the cloud-free molecular atmosphere – i.e. either h < h<sub>b</sub> or h > h<sub>end</sub>. Moreover, the apparent optical thickness η · τ<sub>p</sub> (h<sub>b</sub>, h<sub>end</sub>) of the cloud can be easily computed (see e.g. Garnier et al., 2015):

$$\eta \cdot \tau_{p}(h_{b}, h_{\text{end}}) = -0.5 \cdot \ln \left[ \frac{R'_{MS}(h_{2})}{R'_{MS}(h_{1})} \right], \quad (14)$$

where h<sub>1</sub> and h<sub>2</sub> can be any points satisfying the conditions h<sub>1</sub> < h<sub>b</sub> and h<sub>2</sub> > h<sub>end</sub>.

## 2.2 Simulation software and conditions

Our tool to perform Monte Carlo simulations of lidar signals was the McRALI (Monte Carlo modeling of RAdar and LIidar signals) software developed at the Laboratoire de Météorologie Physique (Szczap et al., 2021; Alkasem et al., 2017). The software employs a forward Monte Carlo (MC) approach, along with the locate estimates method, to simulate the propagation of radiation (see e.g. Marchuk et al., 2013). McRALI is based on 3DMCPOL (three-dimensional

polarized Monte Carlo atmospheric radiative transfer model; Cornet et al., 2010). The polarization state of the radiation is computed using Stokes vectors and scattering matrices of atmospheric compounds. It takes into account molecular scattering. In this work, the properties of the atmosphere were assigned according to the 1976 US standard atmosphere (NOAA, 1976). McRALI is a fully 3D software; that is, the values of the extinction coefficients, the single scattering albedos, and the scattering matrices are assigned in 3D space. Moreover, the mixture of different types of aerosols and/or clouds is allowed. The position of a lidar can be anywhere within or outside of the atmosphere; that is, space-borne, airborne, and ground-based measurement conditions can be simulated. A user can assign a lidar beam direction, a receiver field of view (RFOV), a Stokes vector, and a divergence of the emitted light. The McRALI software was thoroughly tested against data available in the literature (see Appendix in Alkasem et al., 2017). We perform tests when new data are published (under the condition that a paper provides all input data necessary to reproduce simulations). For example, we obtained very good agreement with Fig. 4 of the work by Wang et al. (2021) (including the linear and the circular polarization degree). The good agreement is with data for both ground-based and space-borne lidars.

In this work, MC data were computed so that photons were integrated over the range gate of 20 m; i.e. they correspond to photon counting mode. Such a small value of the range gate was chosen with the aim to study multiple scattering in detail regardless of the fact that it does not correspond to real lidar systems. In other words, the spatial resolution of our MC data is 20 m. The orders of scattering of up to 20 were considered in MC simulations to compute the total multiple scattering. (We have verified that the difference between data obtained with 20 and 10 orders of scattering was not statistically significant for the simulation conditions of this work.)

The results of this work are complementary to the data of the work by Shcherbakov et al. (2022) since most of our MC simulations were performed for the cases of (i) the same water cloud and (ii) the same jet-stream cirrus cloud. The corresponding normalized phase functions, f<sub>W</sub>(θ) and f<sub>JSC</sub>(θ), where θ is the scattering angle, are shown in Fig. 1 by blue and black curves, respectively; their behaviour at forward and backward angles can be seen in the insets. The scattering matrix of water cloud was computed according to the Mie theory for water spheres with the gamma size distribution, with an effective diameter of d<sub>eff</sub> = 18.0 µm (with a standard deviation of 5.3 µm). The SS characteristics of cirrus cloud were computed using the improved geometric optics method (Yang and Liou, 1996); the size distribution of particles was taken to be the gamma distribution, with d<sub>eff</sub> = 56.8 µm (with a standard deviation of 20.1 µm).

To gain a better understanding of the set of parameters that govern the MS effect, we use artificial phase functions. We use the term “chimerical” for them to underline that they are not expected to fit a real phase function but to reproduce
data; R (h) = 1 for the molecular atmosphere. The attenuated scattering ratio can be defined as (see e.g. Winker et al., 2009)

$$R'_{MS}(h) = S_{MS}(h)/S_{m}(h) = R(h) \cdot T_{pA}^{2}(h). \quad (7)$$

The definitions of a MS lidar signal, S<sub>MS</sub>(h), with Eq. (5); of a SS lidar signal, S<sub>1</sub>(h), with Eq. (1); and of the attenuated scattering ratio with Eq. (7) lead to the following properties for the intervals of the cloud-free molecular atmosphere:

$$R'_{MS}(h) = T_{pA}^{2}, \quad (8)$$

$$S_{1}(h) = [\beta_{m}(h)] \cdot T_{m}^{2}(h) \cdot T_{p}^{2} = S_{m}(h) \cdot T_{p}^{2}, \quad (9)$$

$$S_{MS}(h) = [\beta_{m}(h)] \cdot T_{m}^{2}(h) \cdot T_{pA}^{2} = S_{m}(h) \cdot T_{pA}^{2}, \quad (10)$$

where T<sub>p</sub><sup>2</sup> = T<sub>pA</sub><sup>2</sup> = 1 when h < h<sub>b</sub>; T<sub>p</sub><sup>2</sup> = const and T<sub>pA</sub><sup>2</sup> = const when h > h<sub>end</sub>, where h<sub>end</sub> is the distance from the cloud far edge.

It follows from the relationships above that

$$R_{MS101}(h) = 0, R'_{MS}(h) = 1 \quad (11)$$

when h < h<sub>b</sub> and

$$R_{MS101}(h) = \left( \frac{T_{pA}^2}{T_p^2} \right) - 1 = \text{const}, R'_{MS}(h) = T_{pA}^2 = \text{const} \quad (12)$$

when h > h<sub>end</sub>.

In addition, the following relationship can be useful for the interpretation of Monte Carlo data:

$$R'_{MS}(h) = [R_{MS101}(h) + 1] \cdot T_{p}^{2} = \text{const} \text{ when } h > h_{\text{end}}. \quad (13)$$

It should be underlined that if measurement errors are neglected, the functions R<sub>MS101</sub>(h) and R'<sub>MS</sub>(h) are expected to be constant when we are dealing within the intervals of the cloud-free molecular atmosphere – i.e. either h < h<sub>b</sub> or h > h<sub>end</sub>. Moreover, the apparent optical thickness η · τ<sub>p</sub>(h<sub>b</sub>, h<sub>end</sub>) of the cloud can be easily computed (see e.g. Garnier et al., 2015):

$$\eta \cdot \tau_{p}(h_{b}, h_{\text{end}}) = -0.5 \cdot \ln \left[ R'_{MS}(h_{2}) / R'_{MS}(h_{1}) \right], \quad (14)$$

where h<sub>1</sub> and h<sub>2</sub> can be any points satisfying the conditions h<sub>1</sub> < h<sub>b</sub> and h<sub>2</sub> > h<sub>end</sub>.

## 2.2 Simulation software and conditions

Our tool to perform Monte Carlo simulations of lidar signals was the McRALI (Monte Carlo modeling of RAdar and LIidar signals) software developed at the Laboratoire de Météorologie Physique (Szczap et al., 2021; Alkasem et al., 2017). The software employs a forward Monte Carlo (MC) approach, along with the locate estimates method, to simulate the propagation of radiation (see e.g. Marchuk et al., 2013). McRALI is based on 3DMCPOL (three-dimensional

polarized Monte Carlo atmospheric radiative transfer model; Cornet et al., 2010). The polarization state of the radiation is computed using Stokes vectors and scattering matrices of atmospheric compounds. It takes into account molecular scattering. In this work, the properties of the atmosphere were assigned according to the 1976 US standard atmosphere (NOAA, 1976). McRALI is a fully 3D software; that is, the values of the extinction coefficients, the single scattering albedos, and the scattering matrices are assigned in 3D space. Moreover, the mixture of different types of aerosols and/or clouds is allowed. The position of a lidar can be anywhere within or outside of the atmosphere; that is, space-borne, airborne, and ground-based measurement conditions can be simulated. A user can assign a lidar beam direction, a receiver field of view (RFOV), a Stokes vector, and a divergence of the emitted light. The McRALI software was thoroughly tested against data available in the literature (see Appendix in Alkasem et al., 2017). We perform tests when new data are published (under the condition that a paper provides all input data necessary to reproduce simulations). For example, we obtained very good agreement with Fig. 4 of the work by Wang et al. (2021) (including the linear and the circular polarization degree). The good agreement is with data for both ground-based and space-borne lidars.

In this work, MC data were computed so that photons were integrated over the range gate of 20 m; i.e. they correspond to photon counting mode. Such a small value of the range gate was chosen with the aim to study multiple scattering in detail regardless of the fact that it does not correspond to real lidar systems. In other words, the spatial resolution of our MC data is 20 m. The orders of scattering of up to 20 were considered in MC simulations to compute the total multiple scattering. (We have verified that the difference between data obtained with 20 and 10 orders of scattering was not statistically significant for the simulation conditions of this work.)

The results of this work are complementary to the data of the work by Shcherbakov et al. (2022) since most of our MC simulations were performed for the cases of (i) the same water cloud and (ii) the same jet-stream cirrus cloud. The corresponding normalized phase functions, f<sub>W</sub>(θ) and f<sub>JSC</sub>(θ), where θ is the scattering angle, are shown in Fig. 1 by blue and black curves, respectively; their behaviour at forward and backward angles can be seen in the insets. The scattering matrix of water cloud was computed according to the Mie theory for water spheres with the gamma size distribution, with an effective diameter of d<sub>eff</sub> = 18.0 µm (with a standard deviation of 5.3 µm). The SS characteristics of cirrus cloud were computed using the improved geometric optics method (Yang and Liou, 1996); the size distribution of particles was taken to be the gamma distribution, with d<sub>eff</sub> = 56.8 µm (with a standard deviation of 20.1 µm).

To gain a better understanding of the set of parameters that govern the MS effect, we use artificial phase functions. We use the term “chimerical” for them to underline that they are not expected to fit a real phase function but to reproduce
some of its characteristics. A chimerical phase function is a sum of three Gaussian functions, $G_i(\theta)$, $i = 1, 2, 3$: a narrow peak, $G_1(\theta)$, for scattering angles close to 0° (forward direction); a somewhat smoothly varying function, $G_2(\theta)$; and a peak, $G_3(\pi - \theta)$, for angles close to 180°. (We use the notation $\pi - \theta$ in order to underline that we are dealing with the peak in the backward direction.) The first two components are in accordance with the work by Weinman (1968). The third component is inspired by the results of the work (Zhou and Yang, 2015) where rigorous numerical simulations based on solving Maxwell's equations showed that a backscattering peak exists even in cases of randomly oriented ice crystals. We use the formula $G_i(\theta) = a_i \cdot \exp[-(\theta^2)/(\theta_{s,i}^2)]$ to define the Gaussian component, $G_i(\theta)$, where $\theta_{s,i}$ is the 1/e angular half width. The utility of a chimerical phase function consists of the possibility of varying one of the parameters, whereas other characteristics remain unchanged.

The first chimerical phase function, $f_{\text{Ch1}}(\theta) = \sum_{i=1}^{3} G_i(\theta)$, is shown by the red curve in Fig. 1. It was designed to meet the following properties of the normalized phase function, $f_W(\theta)$, of the water cloud. $f_{\text{Ch1}}(0) = f_W(0)$, $f_{\text{Ch1}}(\pi) = f_W(\pi)$, and both phase functions have the same value of the lidar ratio. The width of the Gaussian component $G_1(\theta)$ of $f_{\text{Ch1}}(\theta)$ is assigned so that the function $f_{\text{Ch1}}(\theta)\sin(\theta)$ has the maximum at the same value of $\theta_{\text{max}} = 10.82$ mrad as the function $f_W(\theta)\sin(\theta)$ does. The width of $G_3(\pi - \theta)$ was adjusted so that $f_{\text{Ch1}}(\theta)$ and $f_W(\theta)$ lead to the coincident function $R_{2\text{to}1}(h)$ (see Sect. 3.1). The width of the Gaussian component $G_2(\theta)$ is large; it ensures that a sufficient proportion of photons is scattered sideward. The chimerical phase function $f_{\text{Ch2}}(\theta)$ is shown by the green curve in Fig. 1. It has the same components $G_1(\theta)$ and $G_2(\theta)$ as $f_{\text{Ch1}}(\theta)$. The only difference is $G_3(\pi - \theta)$, which is larger by a factor of 15. Consequently, the curves of $f_{\text{Ch1}}(\theta)$ and $f_{\text{Ch2}}(\theta)$ coincide in Fig. 1 until the scattering angle about 170°. The values of parameters $a_i$ and $\theta_{s,i}$ of the chimerical phase functions are given in Table 1.

Along with MC data, we provide the results, which we obtained using our code based on the Eloranta model (EM) (Eloranta, 1998) (see details in Appendix A). In what follows, the data of our simulations using the EM are referred to as multiple scattering or are labelled MS when up to five orders of scattering were considered. The contribution of the fifth order is sufficiently small in all cases considered below. The addition of the sixth order unreasonably increases the time of computing with the EM.

In this work, each cloud layer is plane-parallel (unless otherwise stated) and homogeneous. All scattering matrices were assigned an angular resolution of 0.01° (about 0.175 mrad); the emitter wavelength $\lambda$ is 0.532 µm.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Normalized phase functions – water cloud: blue lines, jet-stream cirrus: black lines, chimerical phase function $f_{\text{Ch1}}(\theta)$: red lines, chimerical phase function $f_{\text{Ch2}}(\theta)$: green lines.</figcaption></figure>

## 3 Ground-based lidar

The results of this section were obtained for a ground-based lidar, which is at the altitude of $H = 0$ km, and the distance to the cloud base is 8 km. The full receiver field of view (RFOV) is 1.0 mrad (except in Fig. 3c and d); the full emitter field of view (EFOV) is 0.14 mrad. The emitter wavelength, $\lambda$, is 0.532 µm. (We used the characteristics of the lidar system that is in operation at Clermont-Ferrand; Freville et al., 2015.) In order to assure good statistical quality of our Monte Carlo modelling, each signal was simulated with $2 \times 10^{11}$ photons emitted by the lidar (with $4 \times 10^{11}$ photons for the cirrus clouds). Simulations of signals were performed for the orders of scattering of $n = 1$ (single scattering) and $n = 2$ (double scattering) and multiple scattering with an $n$ of 20.

### 3.1 Single-layer cloud

The single-layer cloud is within the altitude range of $H \in [8, 11]$ km; that is, it has an optical thickness of $\tau_p = 3.0$ (with an extinction coefficient of $\epsilon_p = 1.0 \text{ km}^{-1}$). Black and red points in Fig. 2 show the results of our MC simulations reported in terms of the ratios $R_{\text{MSto1}}(d)$ and $R_{\text{2to1}}(d)$, respectively. The parameter $d$ is the distance measured from the cloud base; i.e. the cloud is within the range of $d \in [0, 3]$ km. The interval $d \in (3, 6]$ km is the cloud-free molecular atmosphere – i.e. $H \in (11, 14]$ km. Figure 2a and b correspond to the water and cirrus cloud, respectively. The relative contributions of $R_{\text{MSto1}}(d)$ and $R_{\text{2to1}}(d)$, computed using the EM, are shown in Fig. 2 by the green and blue curves, respectively. The MC in-cloud data – i.e. $d \in [0, 3]$ km – were reported and discussed in the work by Shcherbakov et al. (2022). The focus of interest for this work is the cloud-free molecular atmosphere – that is, the interval of $d \in (3, 6]$ km.
some of its characteristics. A chimerical phase function is a sum of three Gaussian functions, $G_i(\theta)$, $i = 1, 2, 3$: a narrow peak, $G_1(\theta)$, for scattering angles close to 0° (forward direction); a somewhat smoothly varying function, $G_2(\theta)$; and a peak, $G_3(\pi - \theta)$, for angles close to 180°. (We use the notation $\pi - \theta$ in order to underline that we are dealing with the peak in the backward direction.) The first two components are in accordance with the work by Weinman (1968). The third component is inspired by the results of the work (Zhou and Yang, 2015) where rigorous numerical simulations based on solving Maxwell's equations showed that a backscattering peak exists even in cases of randomly oriented ice crystals. We use the formula $G_i(\theta) = a_i \cdot \exp[-(\theta^2)/(\theta_{s,i}^2)]$ to define the Gaussian component, $G_i(\theta)$, where $\theta_{s,i}$ is the 1/e angular half width. The utility of a chimerical phase function consists of the possibility of varying one of the parameters, whereas other characteristics remain unchanged.

The first chimerical phase function, $f_{\text{Ch1}}(\theta) = \sum_{i=1}^{3} G_i(\theta)$, is shown by the red curve in Fig. 1. It was designed to meet the following properties of the normalized phase function, $f_W(\theta)$, of the water cloud. $f_{\text{Ch1}}(0) = f_W(0)$, $f_{\text{Ch1}}(\pi) = f_W(\pi)$, and both phase functions have the same value of the lidar ratio. The width of the Gaussian component $G_1(\theta)$ of $f_{\text{Ch1}}(\theta)$ is assigned so that the function $f_{\text{Ch1}}(\theta)\sin(\theta)$ has the maximum at the same value of $\theta_{\text{max}} = 10.82$ mrad as the function $f_W(\theta)\sin(\theta)$ does. The width of $G_3(\pi - \theta)$ was adjusted so that $f_{\text{Ch1}}(\theta)$ and $f_W(\theta)$ lead to the coincident function $R_{2\text{to}1}(h)$ (see Sect. 3.1). The width of the Gaussian component $G_2(\theta)$ is large; it ensures that a sufficient proportion of photons is scattered sideward. The chimerical phase function $f_{\text{Ch2}}(\theta)$ is shown by the green curve in Fig. 1. It has the same components $G_1(\theta)$ and $G_2(\theta)$ as $f_{\text{Ch1}}(\theta)$. The only difference is $G_3(\pi - \theta)$, which is larger by a factor of 15. Consequently, the curves of $f_{\text{Ch1}}(\theta)$ and $f_{\text{Ch2}}(\theta)$ coincide in Fig. 1 until the scattering angle about 170°. The values of parameters $a_i$ and $\theta_{s,i}$ of the chimerical phase functions are given in Table 1.

Along with MC data, we provide the results, which we obtained using our code based on the Eloranta model (EM) (Eloranta, 1998) (see details in Appendix A). In what follows, the data of our simulations using the EM are referred to as multiple scattering or are labelled MS when up to five orders of scattering were considered. The contribution of the fifth order is sufficiently small in all cases considered below. The addition of the sixth order unreasonably increases the time of computing with the EM.

In this work, each cloud layer is plane-parallel (unless otherwise stated) and homogeneous. All scattering matrices were assigned an angular resolution of 0.01° (about 0.175 mrad); the emitter wavelength $\lambda$ is 0.532 µm.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Normalized phase functions – water cloud: blue lines, jet-stream cirrus: black lines, chimerical phase function $f_{\text{Ch1}}(\theta)$: red lines, chimerical phase function $f_{\text{Ch2}}(\theta)$: green lines.</figcaption></figure>

## 3 Ground-based lidar

The results of this section were obtained for a ground-based lidar, which is at the altitude of $H = 0$ km, and the distance to the cloud base is 8 km. The full receiver field of view (RFOV) is 1.0 mrad (except in Fig. 3c and d); the full emitter field of view (EFOV) is 0.14 mrad. The emitter wavelength, $\lambda$, is 0.532 µm. (We used the characteristics of the lidar system that is in operation at Clermont-Ferrand; Freville et al., 2015.) In order to assure good statistical quality of our Monte Carlo modelling, each signal was simulated with $2 \times 10^{11}$ photons emitted by the lidar (with $4 \times 10^{11}$ photons for the cirrus clouds). Simulations of signals were performed for the orders of scattering of $n = 1$ (single scattering) and $n = 2$ (double scattering) and multiple scattering with an $n$ of 20.

### 3.1 Single-layer cloud

The single-layer cloud is within the altitude range of $H \in [8, 11]$ km; that is, it has an optical thickness of $\tau_p = 3.0$ (with an extinction coefficient of $\epsilon_p = 1.0 \text{ km}^{-1}$). Black and red points in Fig. 2 show the results of our MC simulations reported in terms of the ratios $R_{\text{MSto1}}(d)$ and $R_{\text{2to1}}(d)$, respectively. The parameter $d$ is the distance measured from the cloud base; i.e. the cloud is within the range of $d \in [0, 3]$ km. The interval $d \in (3, 6]$ km is the cloud-free molecular atmosphere – i.e. $H \in (11, 14]$ km. Figure 2a and b correspond to the water and cirrus cloud, respectively. The relative contributions of $R_{\text{MSto1}}(d)$ and $R_{\text{2to1}}(d)$, computed using the EM, are shown in Fig. 2 by the green and blue curves, respectively. The MC in-cloud data – i.e. $d \in [0, 3]$ km – were reported and discussed in the work by Shcherbakov et al. (2022). The focus of interest for this work is the cloud-free molecular atmosphere – that is, the interval of $d \in (3, 6]$ km.
<table><caption>**Table 1.** Values of the input parameters of the chimerical phase functions.</caption><thead><tr><th></th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th><th>θ<sub>s,1</sub><br>(°)</th><th>θ<sub>s,2</sub><br>(°)</th><th>θ<sub>s,3</sub><br>(°)</th></tr></thead><tbody><tr><td>f<sub>Ch1</sub> (θ)</td><td>2.2690</td><td>0.062</td><td>4.086 × 10<sup>-5</sup></td><td>0.88</td><td>200.0</td><td>0.327</td></tr><tr><td>f<sub>Ch2</sub> (θ)</td><td>2.2712</td><td>0.062</td><td>6.26045 × 10<sup>-4</sup></td><td>0.88</td><td>200.0</td><td>5.0</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Monte Carlo simulations of multiple-scattering, $R_{MS101}$ (black points), and double-scattering, $R_{2101}$ (red points), relative contributions to lidar signals. (a) Water cloud. (b) Jet-stream cirrus. Eloranta model simulations are shown by the green (MS) and the blue (DS) curves.</figcaption></figure>

Our MC simulations reveal the following features. The functions $R_{MS101}(d)$ and $R_{2101}(d)$ have undergone a stepwise jump immediately beyond the cloud far edge. It should be kept in mind that lidar signals $S_1(h)$, $S_2(h)$, and $S_{MS}(h)$ sharply decrease immediately beyond the cloud far edge (see e.g. Fig. 7 below). The decrease of $S_1(h)$ is faster; therefore, the ratios $R_{MS101}(d)$ and $R_{2101}(d)$ have undergone a stepwise jump. After the jump, the ratios decrease and tend to zero. It should be noted that the stepwise jump with subsequent decreasing was already reported in the work by Flesia and Starkov (1996), where multiple scattering to a space-borne lidar return from a clear molecular atmosphere obscured by transparent upper-level crystal clouds was assessed by the use of the Monte Carlo technique.

It is seen in Fig. 2a and b that the EM is able to simulate $R_{MS101}(d)$ and $R_{2101}(d)$ within the cloud layer of $d \in [0, 3]$ km well. The simulation results can be considered acceptable in the range of $d \in [3, 6]$ km. That is, our EM data show the stepwise jump of the same amplitude as the MC data immediately beyond the cloud far edge. In contrast, the rate of decrease of $R_{MS101}(d)$ and $R_{2101}(d)$ is lower. Therefore, the EM slightly overestimates the multiple-scattering contribution in the cloud-free molecular atmosphere. We achieved the good agreement with the MC data by adjusting the EM parameters only for the interval of $d \in [0, 3]$ km (see details and the definition of parameters in Appendix A). We have not succeeded in improving the fitting quality for the range of $d \in (3, 6]$ km by varying the fraction $\gamma(h)$ of the energy in the forward peak of the phase function. Thus, all EM data of this work were obtained with $\gamma(h) = 1/2$. The key fact is that the fittings of the water cloud case were

obtained using values of the ratio $P_n(\pi, h)/P_1(\pi, h)$, $n = 2, \dots, 5$, of about 0.5 or lower (see discussion in Appendix A). As for the cloud-free molecular atmosphere, the ratio $P_n(\pi, h)/P_1(\pi, h)$ is “equal to 1.0 due to the broad nature of the molecular phase function near the backscatter direction” (Eloranta, 1998; Whiteman et al., 2001).

### 3.1.1 Stepwise jump and escape effect

It is of importance to underline that the discussion in this section is mainly based on the data obtained under conditions of the double-scattering (DS) approximation. Namely, the ratio $R_{2101}(d)$ of the water cloud case is used as the base to explain the stepwise jump and the escape effect for the following reasons. Both effects are well pronounced in the data for $R_{2101}(d)$ in Fig. 2a; the DS accounts for more than two-thirds of the multiple scattering; and, last but not least, the DS can be understood intuitively.

The properties of $R_{2101}(d)$ – i.e. of double scattering – within the interval $d \in (3, 6]$ km cannot be due to the stretching of pulse length. The definition and a good explanation of the pulse stretching can be found in Miller and Stephens (1999). The simplified explanation could be as follows. In the case of double scattering, a photon goes through two scatterings within the cloud, returns to the receiver, and has the round-trip distance equal to the case of the single scattering from the range of the free atmosphere beyond the cloud far edge. The maximal round-trip distance depends on the configuration geometry – i.e. on the distance from the lidar to the cloud, the cloud depth, the EFOV, and the RFOV. In the case of Fig. 2, the EFOV and the RFOV are narrow, whereas the distance from the lidar to particles of the cloud layer is quite short. The round-trip distance of a double-scattered photon can gain a few metres in such conditions. It means that only the range of $d \in (3, 3.02]$ km of $R_{2101}(d)$ can be somewhat affected by the pulse-length stretching. (See Fig. S1 and the explanation in the Supplement.) In addition, we have to underline the following. The EM uses the assumption that “the multiply scattered photons are scattered from the same slab as the single-scattered photons” (see p. 2466 in Eloranta, 1998). To put it differently, the EM ignores the pulse stretching. Nevertheless, it reproduces, with good accuracy, the stepwise jumps in Fig. 2 on the basis of the phase-function parameters. Therefore, it is safe to assume that the stepwise jump of $R_{MS101}(d)$ and $R_{2101}(d)$ is due to the stepwise jump in phase-function properties at an-
<table><caption>**Table 1.** Values of the input parameters of the chimerical phase functions.</caption><thead><tr><th></th><th>a<sub>1</sub></th><th>a<sub>2</sub></th><th>a<sub>3</sub></th><th>θ<sub>s,1</sub><br>(°)</th><th>θ<sub>s,2</sub><br>(°)</th><th>θ<sub>s,3</sub><br>(°)</th></tr></thead><tbody><tr><td>f<sub>Ch1</sub> (θ)</td><td>2.2690</td><td>0.062</td><td>4.086 × 10<sup>-5</sup></td><td>0.88</td><td>200.0</td><td>0.327</td></tr><tr><td>f<sub>Ch2</sub> (θ)</td><td>2.2712</td><td>0.062</td><td>6.26045 × 10<sup>-4</sup></td><td>0.88</td><td>200.0</td><td>5.0</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Monte Carlo simulations of multiple-scattering, $R_{MS101}$ (black points), and double-scattering, $R_{2101}$ (red points), relative contributions to lidar signals. (a) Water cloud. (b) Jet-stream cirrus. Eloranta model simulations are shown by the green (MS) and the blue (DS) curves.</figcaption></figure>

Our MC simulations reveal the following features. The functions $R_{MS101}(d)$ and $R_{2101}(d)$ have undergone a stepwise jump immediately beyond the cloud far edge. It should be kept in mind that lidar signals $S_1(h)$, $S_2(h)$, and $S_{MS}(h)$ sharply decrease immediately beyond the cloud far edge (see e.g. Fig. 7 below). The decrease of $S_1(h)$ is faster; therefore, the ratios $R_{MS101}(d)$ and $R_{2101}(d)$ have undergone a stepwise jump. After the jump, the ratios decrease and tend to zero. It should be noted that the stepwise jump with subsequent decreasing was already reported in the work by Flesia and Starkov (1996), where multiple scattering to a space-borne lidar return from a clear molecular atmosphere obscured by transparent upper-level crystal clouds was assessed by the use of the Monte Carlo technique.

It is seen in Fig. 2a and b that the EM is able to simulate $R_{MS101}(d)$ and $R_{2101}(d)$ within the cloud layer of $d \in [0, 3]$ km well. The simulation results can be considered acceptable in the range of $d \in [3, 6]$ km. That is, our EM data show the stepwise jump of the same amplitude as the MC data immediately beyond the cloud far edge. In contrast, the rate of decrease of $R_{MS101}(d)$ and $R_{2101}(d)$ is lower. Therefore, the EM slightly overestimates the multiple-scattering contribution in the cloud-free molecular atmosphere. We achieved the good agreement with the MC data by adjusting the EM parameters only for the interval of $d \in [0, 3]$ km (see details and the definition of parameters in Appendix A). We have not succeeded in improving the fitting quality for the range of $d \in (3, 6]$ km by varying the fraction $\gamma(h)$ of the energy in the forward peak of the phase function. Thus, all EM data of this work were obtained with $\gamma(h) = 1/2$. The key fact is that the fittings of the water cloud case were

obtained using values of the ratio $P_n(\pi, h)/P_1(\pi, h)$, $n = 2, \dots, 5$, of about 0.5 or lower (see discussion in Appendix A). As for the cloud-free molecular atmosphere, the ratio $P_n(\pi, h)/P_1(\pi, h)$ is “equal to 1.0 due to the broad nature of the molecular phase function near the backscatter direction” (Eloranta, 1998; Whiteman et al., 2001).

### 3.1.1 Stepwise jump and escape effect

It is of importance to underline that the discussion in this section is mainly based on the data obtained under conditions of the double-scattering (DS) approximation. Namely, the ratio $R_{2101}(d)$ of the water cloud case is used as the base to explain the stepwise jump and the escape effect for the following reasons. Both effects are well pronounced in the data for $R_{2101}(d)$ in Fig. 2a; the DS accounts for more than two-thirds of the multiple scattering; and, last but not least, the DS can be understood intuitively.

The properties of $R_{2101}(d)$ – i.e. of double scattering – within the interval $d \in (3, 6]$ km cannot be due to the stretching of pulse length. The definition and a good explanation of the pulse stretching can be found in Miller and Stephens (1999). The simplified explanation could be as follows. In the case of double scattering, a photon goes through two scatterings within the cloud, returns to the receiver, and has the round-trip distance equal to the case of the single scattering from the range of the free atmosphere beyond the cloud far edge. The maximal round-trip distance depends on the configuration geometry – i.e. on the distance from the lidar to the cloud, the cloud depth, the EFOV, and the RFOV. In the case of Fig. 2, the EFOV and the RFOV are narrow, whereas the distance from the lidar to particles of the cloud layer is quite short. The round-trip distance of a double-scattered photon can gain a few metres in such conditions. It means that only the range of $d \in (3, 3.02]$ km of $R_{2101}(d)$ can be somewhat affected by the pulse-length stretching. (See Fig. S1 and the explanation in the Supplement.) In addition, we have to underline the following. The EM uses the assumption that “the multiply scattered photons are scattered from the same slab as the single-scattered photons” (see p. 2466 in Eloranta, 1998). To put it differently, the EM ignores the pulse stretching. Nevertheless, it reproduces, with good accuracy, the stepwise jumps in Fig. 2 on the basis of the phase-function parameters. Therefore, it is safe to assume that the stepwise jump of $R_{MS101}(d)$ and $R_{2101}(d)$ is due to the stepwise jump in phase-function properties at an-
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Monte Carlo simulations of multiple-scattering, $R_{MSto1}$ (black points), and double-scattering, $R_{2to1}$ (red points), relative contributions to lidar signals. In panels (a) and (b), RFOV = 1 mrad. In panels (c) and (d), RFOV = 110 mrad. (a) Chimerical phase function $f_{Ch1}(\theta)$. Panels (b) and (d) show the chimerical phase function $f_{Ch2}(\theta)$. (c) Water cloud.</figcaption></figure>

gles close to 180° (between the phase function of particles within the cloud and the Rayleigh scattering within the free atmosphere).

That assumption is confirmed by the plots in Fig. 3a and b, where the results of MC simulations are shown as the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$ for the cases of the chimerical phase functions $f_{Ch1}(\theta)$ and $f_{Ch2}(\theta)$, respectively. We adjusted the width of the component $G_3(\pi - \theta)$ of $f_{Ch1}(\theta)$ so that MC simulations with $f_{Ch1}(\theta)$ give the same $R_{2to1}(d)$ within the range of $d \in [0, 3]$ km as $R_{2to1}(d)$ of the water cloud in Fig. 2a. It turns out that the peak in the backward direction of the water cloud is a little bit larger than $G_3(\pi - \theta)$ of $f_{Ch1}(\theta)$ (see Fig. 1). Nevertheless, the MC data in Fig. 3a are coincident with the MC data in Fig. 2a for the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$ within the full range of $d \in [0, 6]$ km. The good agreement is due to the fact that one of the key parameters of multiple scattering is a weighted average of a phase function near the backscatter direction (see Eq. 10 of Eloranta, 1998; see Fig. S2 and the intuitive explanation in the Supplement). That average depends on the properties of the forward and backward peaks of the phase function and on the receiver field of view. This means that we fitted the value of the weighted average when we adjusted the width of the Gaussian function $G_3(\pi - \theta)$.

The chimerical phase functions, $f_{Ch1}(\theta)$ and $f_{Ch2}(\theta)$, have the same values of the parameters $\theta_{s,1}$ and $\theta_{s,2}$, but the width $\theta_{s,3}$ is 15 times higher for $f_{Ch2}(\theta)$. That leads to marked distinctions in the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$ (see Fig. 3b). There is no stepwise jump immediately beyond the cloud far edge. It means that the component $G_3(\pi - \theta)$ of $f_{Ch2}(\theta)$ is large enough to have the weighted average equal to $f_{Ch2}(\pi)$, as in the case of the Rayleigh phase function. If

we use the terms of the Eloranta model (Eloranta, 1998), the ratio $P_2(\pi, h)/P_1(\pi, h) = 1$ for the phase function $f_{Ch2}(\theta)$ – i.e. it has the same value as in the case of the Rayleigh phase function. To put it differently, a higher proportion of photons are scattered by the cloud in the backward direction within the RFOV and contributes to lidar signals in the MS case and under the DS approximation. As a consequence, the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$ of Fig. 3b are much higher than in Figs. 2a and 3a for the in-cloud range of $d \in [0, 3]$ km.

The MC data of Figs. 2a and 3a, b are coincident within the interval of the cloud-free molecular atmosphere, $d \in (3, 6]$ km, if they are superimposed in the same figure. Thus, it is safe to assume that the following two properties are defined by the forward-peak width of the phase function and the RFOV: (i) the proportion of scattered photons emerging from a cloud that are subsequently scattered in the backward direction to the receiver and (ii) the rate of decrease of the relative contribution of multiple scattering outside the cloud. (We recall that we are dealing with ground-based lidar, the cloud optical thickness is not high, and the RFOV and the EFOV are narrow.) That assumption explains the difference in the rates of decrease within the range of $d \in [3, 6]$ km between Fig. 2a and Fig. b. The forward peak of the scattering phase function of the cirrus cloud is much stronger compared to the one of the water cloud. That leads to the slower rate of decrease.

In our opinion, those properties are a direct consequence of the fact that the photons scattered within the forward peak escape the sampling volume. (For the sake of simplicity, we will use the term “escape effect” to refer to that phenomenon.) The following arguments justify that statement. Let us consider the case of the forward scattering. If we refer to Fraunhofer diffraction by large spheres (see e.g. Chap. 8.3 of Van de Hulst, 1981), the first dark ring – i.e. the Airy disc – is at the angle of $\theta_{Airy} \approx 36$ mrad (about 2°) when a particle has a diameter of 18 µm and a wavelength, $\lambda$, of 0.532 µm ($\theta_{Airy} \gg \text{RFOV}/2 = 0.5$ mrad). If we refer to the phase function of the water cloud (see Fig. 1), the intervals of the angles of $[0, 0.5]$ mrad and $[0, 36]$ mrad account for 0.059% and 41.3% of the total scattered light, respectively. In such conditions, it can hardly be expected that a large proportion of small-angle forward-scattered photons always remain within the field of view of the detector; what is more likely is that photons escape the sampling volume.

With the aim of confirming our arguments, we performed MC simulations for the cases of larger RFOVs. All but the RFOV parameters of lidar configuration and the water cloud properties are the same as in the case of Fig. 2a. The full RFOV is 110 mrad. Such a large value ensures that the narrow peak $G_1(\theta)$ (forward direction) of the chimerical phase functions is completely within the angle of $\text{RFOV}/2 = 55$ mrad (about 3°). Black and red points in Fig. 3c and d show the results of our MC simulations reported in terms of the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$, respectively. Figure 3c and d show the cases of the water cloud and the chimerical
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Monte Carlo simulations of multiple-scattering, $R_{MSto1}$ (black points), and double-scattering, $R_{2to1}$ (red points), relative contributions to lidar signals. In panels (a) and (b), RFOV = 1 mrad. In panels (c) and (d), RFOV = 110 mrad. (a) Chimerical phase function $f_{Ch1}(\theta)$. Panels (b) and (d) show the chimerical phase function $f_{Ch2}(\theta)$. (c) Water cloud.</figcaption></figure>

gles close to 180° (between the phase function of particles within the cloud and the Rayleigh scattering within the free atmosphere).

That assumption is confirmed by the plots in Fig. 3a and b, where the results of MC simulations are shown as the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$ for the cases of the chimerical phase functions $f_{Ch1}(\theta)$ and $f_{Ch2}(\theta)$, respectively. We adjusted the width of the component $G_3(\pi - \theta)$ of $f_{Ch1}(\theta)$ so that MC simulations with $f_{Ch1}(\theta)$ give the same $R_{2to1}(d)$ within the range of $d \in [0, 3]$ km as $R_{2to1}(d)$ of the water cloud in Fig. 2a. It turns out that the peak in the backward direction of the water cloud is a little bit larger than $G_3(\pi - \theta)$ of $f_{Ch1}(\theta)$ (see Fig. 1). Nevertheless, the MC data in Fig. 3a are coincident with the MC data in Fig. 2a for the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$ within the full range of $d \in [0, 6]$ km. The good agreement is due to the fact that one of the key parameters of multiple scattering is a weighted average of a phase function near the backscatter direction (see Eq. 10 of Eloranta, 1998; see Fig. S2 and the intuitive explanation in the Supplement). That average depends on the properties of the forward and backward peaks of the phase function and on the receiver field of view. This means that we fitted the value of the weighted average when we adjusted the width of the Gaussian function $G_3(\pi - \theta)$.

The chimerical phase functions, $f_{Ch1}(\theta)$ and $f_{Ch2}(\theta)$, have the same values of the parameters $\theta_{s,1}$ and $\theta_{s,2}$, but the width $\theta_{s,3}$ is 15 times higher for $f_{Ch2}(\theta)$. That leads to marked distinctions in the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$ (see Fig. 3b). There is no stepwise jump immediately beyond the cloud far edge. It means that the component $G_3(\pi - \theta)$ of $f_{Ch2}(\theta)$ is large enough to have the weighted average equal to $f_{Ch2}(\pi)$, as in the case of the Rayleigh phase function. If

we use the terms of the Eloranta model (Eloranta, 1998), the ratio $P_2(\pi, h)/P_1(\pi, h) = 1$ for the phase function $f_{Ch2}(\theta)$ – i.e. it has the same value as in the case of the Rayleigh phase function. To put it differently, a higher proportion of photons are scattered by the cloud in the backward direction within the RFOV and contributes to lidar signals in the MS case and under the DS approximation. As a consequence, the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$ of Fig. 3b are much higher than in Figs. 2a and 3a for the in-cloud range of $d \in [0, 3]$ km.

The MC data of Figs. 2a and 3a, b are coincident within the interval of the cloud-free molecular atmosphere, $d \in (3, 6]$ km, if they are superimposed in the same figure. Thus, it is safe to assume that the following two properties are defined by the forward-peak width of the phase function and the RFOV: (i) the proportion of scattered photons emerging from a cloud that are subsequently scattered in the backward direction to the receiver and (ii) the rate of decrease of the relative contribution of multiple scattering outside the cloud. (We recall that we are dealing with ground-based lidar, the cloud optical thickness is not high, and the RFOV and the EFOV are narrow.) That assumption explains the difference in the rates of decrease within the range of $d \in [3, 6]$ km between Fig. 2a and Fig. b. The forward peak of the scattering phase function of the cirrus cloud is much stronger compared to the one of the water cloud. That leads to the slower rate of decrease.

In our opinion, those properties are a direct consequence of the fact that the photons scattered within the forward peak escape the sampling volume. (For the sake of simplicity, we will use the term “escape effect” to refer to that phenomenon.) The following arguments justify that statement. Let us consider the case of the forward scattering. If we refer to Fraunhofer diffraction by large spheres (see e.g. Chap. 8.3 of Van de Hulst, 1981), the first dark ring – i.e. the Airy disc – is at the angle of $\theta_{Airy} \approx 36$ mrad (about 2°) when a particle has a diameter of 18 µm and a wavelength, $\lambda$, of 0.532 µm ($\theta_{Airy} \gg \text{RFOV}/2 \approx 0.5$ mrad). If we refer to the phase function of the water cloud (see Fig. 1), the intervals of the angles of $[0, 0.5]$ mrad and $[0, 36]$ mrad account for 0.059% and 41.3% of the total scattered light, respectively. In such conditions, it can hardly be expected that a large proportion of small-angle forward-scattered photons always remain within the field of view of the detector; what is more likely is that photons escape the sampling volume.

With the aim of confirming our arguments, we performed MC simulations for the cases of larger RFOVs. All but the RFOV parameters of lidar configuration and the water cloud properties are the same as in the case of Fig. 2a. The full RFOV is 110 mrad. Such a large value ensures that the narrow peak $G_1(\theta)$ (forward direction) of the chimerical phase functions is completely within the angle of $\text{RFOV}/2 \approx 55$ mrad (about 3°). Black and red points in Fig. 3c and d show the results of our MC simulations reported in terms of the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$, respectively. Figure 3c and d show the cases of the water cloud and the chimerical
phase function $f_{Ch2}(\theta)$, respectively. (Note the log scale on the y axis in Fig. 3c and d.)

It is seen that the ratios $R_{2to1}(d)$ have undergone a very strong stepwise jump immediately beyond the cloud far edge. If we consider the configuration geometry – i.e. the angle of RFOV/2 = 55 mrad – and the distance from the lidar to the cloud far edge of 11 km, the pulse stretching is the cause of the jump. The pulse stretching strongly affects $R_{2to1}(d)$ in Fig. 3d within the range of $d \in (3, 3.25]$ km. The ratio $R_{2to1}(d)$ is almost constant within the interval of $d \in [3.25, 6]$ km. (In our opinion, the slight decrease is due to the pulse-stretching effect.) The properties of the ratio $R_{2to1}(d)$ in Fig. 3c are the same but less pronounced. We can hypothesize that photons that were forward-scattered within angles of > 3° play a noticeable role in the case of the water cloud. The results shown in Fig. 3c and d mean that only when the RFOV is sufficiently large do small-angle forward-scattered photons remain within the field of view of the detector, and the escape effect is eliminated. As for the cases of multiple scattering, the ratios $R_{MSto1}(d)$ in Fig. 3c and d are very high due to the large RFOV, and the escape effect is evident.

It is seen that the decrease due to the escape effect is an inherent part of $R_{MSto1}(d)$ properties within the free atmosphere beyond the cloud far edge. That property is in direct contradiction with Eqs. (10)–(14) above, which are the consequences of Eq. (5). Therefore, we can conclude that Eq. (5) is an approximate model of lidar signals under MS conditions because it does not take into account the escape effect.

### 3.1.2 Practical aspects

Knowing that the value of the extinction coefficient of $\epsilon_p = 1.0 \text{ km}^{-1}$ is not typical of cirrus clouds, we provide practitioners with estimations of MS effects for lower values of $\epsilon_p$ of the jet-stream cirrus. Black points in Fig. 4 show the results of our MC simulations reported in terms of the ratio $R_{MSto1}(d)$. The parameter $d$ is the distance measured from the cloud base – i.e. the cloud is within the range of $d \in [0, 3]$ km. The interval of $d \in (3, 6]$ km is the cloud-free molecular atmosphere – i.e. the altitude of $H \in (11, 14]$ km. The MC simulations were performed for the same geometry of planet-parallel homogeneous cloud as in the case of Fig. 2b. The only difference consists of the values of the extinction coefficient, which were $\epsilon_p = 0.06$ and $\epsilon_p = 0.2 \text{ km}^{-1}$ in Fig. 4a and b, respectively. It is seen that the lidar data are affected by MS within the cloud ($d \in [0, 3]$ km) and above the cloud far edge ($d \in (3, 6]$ km). The cases in Fig. 4 are characterized by the quite low values of the extinction coefficient – i.e. the low probability of the interaction of a photon with cloud particles. Moreover, the RFOV – i.e. the sampling volume – and the forward peak of the scattering phase function are narrow. All that led to the rather high dispersion of the MC data in Fig. 4 despite the very high number, $4 \times 10^{11}$, of sampled photons. Nevertheless, it is safe to conclude that the MS contribution decreases with the distance from the cloud

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Monte Carlo simulations of multiple-scattering relative contributions, $R_{MSto1}$, to lidar signals from the jet-stream cirrus; the extinction coefficients are (a) $0.06 \text{ km}^{-1}$ and (b) $0.2 \text{ km}^{-1}$.</figcaption></figure>

<table><caption>**Table 2.** Estimated values of the apparent optical thickness, $\eta \cdot \tau_p(h_b, h_{end})$.</caption><thead><tr><th rowspan="2">MC data case (figure)</th><th rowspan="2"></th><th colspan="2">Apparent optical thickness</th><th>Real optical thickness</th></tr><tr><th>0.2</th><th>3.0</th><th></th></tr></thead><tbody><tr><td rowspan="4">Distance from cloud (km)</td><td>Water cloud (Fig. 2a)</td><td>2.83</td><td>2.99</td><td>3.0</td></tr><tr><td>Jet-stream cirrus (Fig. 2b)</td><td>2.58</td><td>2.94</td><td>3.0</td></tr><tr><td>Jet-stream cirrus (Fig. 4a)</td><td>0.15</td><td>0.18</td><td>0.18</td></tr><tr><td>Jet-stream cirrus (Fig. 4b)</td><td>0.51</td><td>0.60</td><td>0.60</td></tr></tbody></table>

far edge due to the escape effect. The MS contribution is lower than 5% within the regions of the cloud-free molecular atmosphere, with the distance from the cloud far edge of about 1 km or higher when $\epsilon_p \le 0.2 \text{ km}^{-1}$.

Table 2 shows the values of the apparent optical thickness, $\eta \cdot \tau_p(h_b, h_{end})$, of the cloud computed using Eq. (13) with two values of $h_2$ assigned within the interval of $d \in (3, 6]$ km (with an altitude of $H \in (11, 14]$ km). The values of 0.2 and 3.0 km mean the distance from the cloud far edge within the cloud-free molecular atmosphere; they correspond to $d = 3.2 \text{ km}$ ($H = 11.2 \text{ km}$) and $d = 6.0 \text{ km}$ ($H = 14.0 \text{ km}$), respectively. It is seen that the computed value of the apparent optical thickness depends on the chosen value of $h_2$. It is close to the value of the real optical thickness, when $h_2 = 14.0 \text{ km}$, and can be about 15% lower, when $h_2 = 11.2 \text{ km}$, in the cases of the jet-stream cirrus.

## 3.2 Two-layered cloud

The two-layered cloud consists of homogeneous layers situated at the altitudes from 8 to 9 km and from 10 to 11 km. Every layer has the optical thickness of $\tau_p = 1.0$; the total optical thickness is $\tau_p = 2.0$. The black and red points in Fig. 5 show the results of our MC simulations reported in terms of the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$, respectively. The parameter $d$ is the distance measured from the cloud base; i.e. the layers are within the ranges of $d \in [0, 1]$ and $d \in [2, 3]$ km. The intervals of $d \in (1, 2)$ and $d \in (3, 6]$ km are the cloud-free molecular atmosphere. Figure 5a and b correspond to
phase function $f_{Ch2}(\theta)$, respectively. (Note the log scale on the y axis in Fig. 3c and d.)

It is seen that the ratios $R_{2to1}(d)$ have undergone a very strong stepwise jump immediately beyond the cloud far edge. If we consider the configuration geometry – i.e. the angle of RFOV/2 = 55 mrad – and the distance from the lidar to the cloud far edge of 11 km, the pulse stretching is the cause of the jump. The pulse stretching strongly affects $R_{2to1}(d)$ in Fig. 3d within the range of $d \in (3, 3.25]$ km. The ratio $R_{2to1}(d)$ is almost constant within the interval of $d \in [3.25, 6]$ km. (In our opinion, the slight decrease is due to the pulse-stretching effect.) The properties of the ratio $R_{2to1}(d)$ in Fig. 3c are the same but less pronounced. We can hypothesize that photons that were forward-scattered within angles of > 3° play a noticeable role in the case of the water cloud. The results shown in Fig. 3c and d mean that only when the RFOV is sufficiently large do small-angle forward-scattered photons remain within the field of view of the detector, and the escape effect is eliminated. As for the cases of multiple scattering, the ratios $R_{MSto1}(d)$ in Fig. 3c and d are very high due to the large RFOV, and the escape effect is evident.

It is seen that the decrease due to the escape effect is an inherent part of $R_{MSto1}(d)$ properties within the free atmosphere beyond the cloud far edge. That property is in direct contradiction with Eqs. (10)–(14) above, which are the consequences of Eq. (5). Therefore, we can conclude that Eq. (5) is an approximate model of lidar signals under MS conditions because it does not take into account the escape effect.

### 3.1.2 Practical aspects

Knowing that the value of the extinction coefficient of $\epsilon_p = 1.0 \text{ km}^{-1}$ is not typical of cirrus clouds, we provide practitioners with estimations of MS effects for lower values of $\epsilon_p$ of the jet-stream cirrus. Black points in Fig. 4 show the results of our MC simulations reported in terms of the ratio $R_{MSto1}(d)$. The parameter $d$ is the distance measured from the cloud base – i.e. the cloud is within the range of $d \in [0, 3]$ km. The interval of $d \in (3, 6]$ km is the cloud-free molecular atmosphere – i.e. the altitude of $H \in (11, 14]$ km. The MC simulations were performed for the same geometry of planet-parallel homogeneous cloud as in the case of Fig. 2b. The only difference consists of the values of the extinction coefficient, which were $\epsilon_p = 0.06$ and $\epsilon_p = 0.2 \text{ km}^{-1}$ in Fig. 4a and b, respectively. It is seen that the lidar data are affected by MS within the cloud ($d \in [0, 3]$ km) and above the cloud far edge ($d \in (3, 6]$ km). The cases in Fig. 4 are characterized by the quite low values of the extinction coefficient – i.e. the low probability of the interaction of a photon with cloud particles. Moreover, the RFOV – i.e. the sampling volume – and the forward peak of the scattering phase function are narrow. All that led to the rather high dispersion of the MC data in Fig. 4 despite the very high number, $4 \times 10^{11}$, of sampled photons. Nevertheless, it is safe to conclude that the MS contribution decreases with the distance from the cloud

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Monte Carlo simulations of multiple-scattering relative contributions, $R_{MSto1}$, to lidar signals from the jet-stream cirrus; the extinction coefficients are (a) $0.06 \text{ km}^{-1}$ and (b) $0.2 \text{ km}^{-1}$.</figcaption></figure>

<table><caption>**Table 2.** Estimated values of the apparent optical thickness, $\eta \cdot \tau_p(h_b, h_{end})$.</caption><thead><tr><th rowspan="2">MC data case (figure)</th><th rowspan="2"></th><th colspan="2">Apparent optical thickness</th><th>Real optical thickness</th></tr><tr><th>0.2</th><th>3.0</th><th></th></tr></thead><tbody><tr><td rowspan="4">Distance from cloud (km)</td><td>Water cloud (Fig. 2a)</td><td>2.83</td><td>2.99</td><td>3.0</td></tr><tr><td>Jet-stream cirrus (Fig. 2b)</td><td>2.58</td><td>2.94</td><td>3.0</td></tr><tr><td>Jet-stream cirrus (Fig. 4a)</td><td>0.15</td><td>0.18</td><td>0.18</td></tr><tr><td>Jet-stream cirrus (Fig. 4b)</td><td>0.51</td><td>0.60</td><td>0.60</td></tr></tbody></table>

far edge due to the escape effect. The MS contribution is lower than 5% within the regions of the cloud-free molecular atmosphere, with the distance from the cloud far edge of about 1 km or higher when $\epsilon_p \le 0.2 \text{ km}^{-1}$.

Table 2 shows the values of the apparent optical thickness, $\eta \cdot \tau_p(h_b, h_{end})$, of the cloud computed using Eq. (13) with two values of $h_2$ assigned within the interval of $d \in (3, 6]$ km (with an altitude of $H \in (11, 14]$ km). The values of 0.2 and 3.0 km mean the distance from the cloud far edge within the cloud-free molecular atmosphere; they correspond to $d = 3.2 \text{ km}$ ($H = 11.2 \text{ km}$) and $d = 6.0 \text{ km}$ ($H = 14.0 \text{ km}$), respectively. It is seen that the computed value of the apparent optical thickness depends on the chosen value of $h_2$. It is close to the value of the real optical thickness, when $h_2 = 14.0 \text{ km}$, and can be about 15% lower, when $h_2 = 11.2 \text{ km}$, in the cases of the jet-stream cirrus.

## 3.2 Two-layered cloud

The two-layered cloud consists of homogeneous layers situated at the altitudes from 8 to 9 km and from 10 to 11 km. Every layer has the optical thickness of $\tau_p = 1.0$; the total optical thickness is $\tau_p = 2.0$. The black and red points in Fig. 5 show the results of our MC simulations reported in terms of the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$, respectively. The parameter $d$ is the distance measured from the cloud base; i.e. the layers are within the ranges of $d \in [0, 1]$ and $d \in [2, 3]$ km. The intervals of $d \in (1, 2)$ and $d \in (3, 6]$ km are the cloud-free molecular atmosphere. Figure 5a and b correspond to
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Monte Carlo simulations of multiple-scattering, $R_{MStol}$ (black points), and double-scattering, $R_{2to1}$ (red points), relative contributions to lidar signals. (a) Water cloud. (b) Jet-stream cirrus. Eloranta model simulations are shown by the green (MS) and the blue (DS) curves.</figcaption></figure>

the water and cirrus cloud, respectively. The relative contributions $R_{MStol}(d)$ and $R_{2to1}(d)$ computed using the EM are shown in Fig. 5 by green and blue lines, respectively.

The features of the MC data within the intervals of the cloud-free molecular atmosphere are closely similar to those addressed in the previous section. In addition, the following property is noteworthy. The rate of decrease of the MS relative contributions in the case of the water cloud is quite fast; therefore, scattered photons emerging from the first layer almost do not affect the lidar signal from the second layer. (The values of the ratios $R_{MStol}(d = 2.0)$ and $R_{2to1}(d = 2.0)$ are close to zero.) As it was discussed above, the rate of decrease of the MS relative contributions in the case of the cirrus clouds is much slower compared to the water cloud case. Therefore, scattered photons emerging from the first layer affect the lidar signal from the second layer so that the ratios $R_{MStol}(d = 2.0)$ and $R_{2to1}(d = 2.0)$ are about 0.25. In other words, the lidar signal from the near edge of the second layer is about 25 % higher compared to the SS approximation. Obviously, the distance between cloud layers is a key parameter. The values of the ratios $R_{MStol}(d)$ and $R_{2to1}(d)$ at the near edge of the second layer vary if that distance changes.

As for our EM simulations, we used the same values of the ratio $P_n(\pi, h)/P_1(\pi, h)$, $n = 2, \dots, 5$, as in the previous section for each cloud layer. We can conclude one more time that the EM is able to simulate MS contributions within the cloud layers well. As for the intervals of the cloud-free molecular atmosphere, the behaviour of $R_{MStol}(d)$ and $R_{2to1}(d)$ is correct, whereas the MS contribution is slightly overestimated.

# 4 Space-borne lidar

The results of this section were obtained for a space-borne lidar, which is at an altitude of $H = 705$ km. The full receiver field of view (RFOV) is 0.13 mrad; the full emitter field of view (EFOV) is 0.1 mrad. The emitter wavelength, $\lambda$, is 0.532 µm. Those values correspond to the technical characteristics of the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) (Young and Vaughan, 2009). In order to

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Monte Carlo simulations of single-scattering (blue points) and multiple-scattering (black points) lidar signals. (a) Water cloud. (b) Jet-stream cirrus.</figcaption></figure>

assure good statistical quality of our Monte Carlo modelling, each signal was simulated with $2 \times 10^{11}$ photons emitted by the lidar (with $4 \times 10^{11}$ photons for the cirrus clouds). Simulations of signals were performed for the orders of scattering of $n = 1$ (single scattering) and $n = 2$ (double scattering) and multiple scattering with an $n$ of 20.

## 4.1 Single-layer cloud

### 4.1.1 Plane-parallel homogeneous cloud

The single-layer cloud is within the altitude range of $H \in [8, 11]$ km and has an optical thickness of $\tau_p = 3.0$ (with an extinction coefficient of $\epsilon_p(h) = 1.0 \text{ km}^{-1}$). For reasons of direct comparison of our data with those addressed in the literature, we show our simulated signals of the space-borne lidar in Fig. 6. The SS $S_1(h)$ and MS $S_{MS}(h)$ lidar signals are shown by the blue and black points, respectively. Qualitatively, Fig. 6a and b are similar to Fig. B1c and d, respectively, of the work by Reverdy et al. (2015) and to Fig. 1c and d of the work by Donovan (2016). The distinction in the quantitative characteristics is due to the difference in the cloud optical thickness and optical parameters of the particles. The feature that the MS signals below the cloud base tend to the SS signals is well pronounced. Multiple-scattering ratio $R_{MStol}(d)$ is more suitable to quantify that feature compared to lidar signals shown in the semi-logarithmic plot.

Black and red points in Fig. 7 show the results of our MC simulations reported in terms of the ratios $R_{MStol}(d)$ and $R_{2to1}(d)$, respectively. The parameter $d$ is the distance measured from the cloud near edge to the lidar – i.e. the cloud is within the range of $d \in [0, 3]$ km. ($d = 0$ km corresponds to an altitude of $H = 11$ km.) The interval of $d \in (3, 11]$ km (an altitude, $H$, from 8 to 0 km) is the cloud-free molecular atmosphere below the cloud. Figure 7a and b correspond to the water and cirrus cloud, respectively. The relative contri-
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Monte Carlo simulations of multiple-scattering, $R_{MStol}$ (black points), and double-scattering, $R_{2to1}$ (red points), relative contributions to lidar signals. (a) Water cloud. (b) Jet-stream cirrus. Eloranta model simulations are shown by the green (MS) and the blue (DS) curves.</figcaption></figure>

the water and cirrus cloud, respectively. The relative contributions $R_{MStol}(d)$ and $R_{2to1}(d)$ computed using the EM are shown in Fig. 5 by green and blue lines, respectively.

The features of the MC data within the intervals of the cloud-free molecular atmosphere are closely similar to those addressed in the previous section. In addition, the following property is noteworthy. The rate of decrease of the MS relative contributions in the case of the water cloud is quite fast; therefore, scattered photons emerging from the first layer almost do not affect the lidar signal from the second layer. (The values of the ratios $R_{MStol}(d = 2.0)$ and $R_{2to1}(d = 2.0)$ are close to zero.) As it was discussed above, the rate of decrease of the MS relative contributions in the case of the cirrus clouds is much slower compared to the water cloud case. Therefore, scattered photons emerging from the first layer affect the lidar signal from the second layer so that the ratios $R_{MStol}(d = 2.0)$ and $R_{2to1}(d = 2.0)$ are about 0.25. In other words, the lidar signal from the near edge of the second layer is about 25 % higher compared to the SS approximation. Obviously, the distance between cloud layers is a key parameter. The values of the ratios $R_{MStol}(d)$ and $R_{2to1}(d)$ at the near edge of the second layer vary if that distance changes.

As for our EM simulations, we used the same values of the ratio $P_n(\pi, h)/P_1(\pi, h)$, $n = 2, \dots, 5$, as in the previous section for each cloud layer. We can conclude one more time that the EM is able to simulate MS contributions within the cloud layers well. As for the intervals of the cloud-free molecular atmosphere, the behaviour of $R_{MStol}(d)$ and $R_{2to1}(d)$ is correct, whereas the MS contribution is slightly overestimated.

# 4 Space-borne lidar

The results of this section were obtained for a space-borne lidar, which is at an altitude of $H = 705$ km. The full receiver field of view (RFOV) is 0.13 mrad; the full emitter field of view (EFOV) is 0.1 mrad. The emitter wavelength, $\lambda$, is 0.532 µm. Those values correspond to the technical characteristics of the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) (Young and Vaughan, 2009). In order to

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Monte Carlo simulations of single-scattering (blue points) and multiple-scattering (black points) lidar signals. (a) Water cloud. (b) Jet-stream cirrus.</figcaption></figure>

assure good statistical quality of our Monte Carlo modelling, each signal was simulated with $2 \times 10^{11}$ photons emitted by the lidar (with $4 \times 10^{11}$ photons for the cirrus clouds). Simulations of signals were performed for the orders of scattering of $n = 1$ (single scattering) and $n = 2$ (double scattering) and multiple scattering with an $n$ of 20.

## 4.1 Single-layer cloud

### 4.1.1 Plane-parallel homogeneous cloud

The single-layer cloud is within the altitude range of $H \in [8, 11]$ km and has an optical thickness of $\tau_p = 3.0$ (with an extinction coefficient of $\epsilon_p(h) = 1.0 \text{ km}^{-1}$). For reasons of direct comparison of our data with those addressed in the literature, we show our simulated signals of the space-borne lidar in Fig. 6. The SS $S_1(h)$ and MS $S_{MS}(h)$ lidar signals are shown by the blue and black points, respectively. Qualitatively, Fig. 6a and b are similar to Fig. B1c and d, respectively, of the work by Reverdy et al. (2015) and to Fig. 1c and d of the work by Donovan (2016). The distinction in the quantitative characteristics is due to the difference in the cloud optical thickness and optical parameters of the particles. The feature that the MS signals below the cloud base tend to the SS signals is well pronounced. Multiple-scattering ratio $R_{MStol}(d)$ is more suitable to quantify that feature compared to lidar signals shown in the semi-logarithmic plot.

Black and red points in Fig. 7 show the results of our MC simulations reported in terms of the ratios $R_{MStol}(d)$ and $R_{2to1}(d)$, respectively. The parameter $d$ is the distance measured from the cloud near edge to the lidar – i.e. the cloud is within the range of $d \in [0, 3]$ km. ($d = 0$ km corresponds to an altitude of $H = 11$ km.) The interval of $d \in (3, 11]$ km (an altitude, $H$, from 8 to 0 km) is the cloud-free molecular atmosphere below the cloud. Figure 7a and b correspond to the water and cirrus cloud, respectively. The relative contri-
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Monte Carlo simulations of multiple-scattering, $R_{MSto1}$ (black points), and double-scattering, $R_{2to1}$ (red points), relative contributions to lidar signals. (a) Water cloud. (b) Jet-stream cirrus. Eloranta model simulations are shown by the green (MS) and the blue (DS) curves.</figcaption></figure>

butions $R_{MSto1}(d)$ and $R_{2to1}(d)$ computed using the EM are shown in Fig. 7 by the green and blue curves, respectively. The MC in-cloud data – i.e. $d \in [0, 3]$ km – were reported and discussed in the work by Shcherbakov et al. (2022). The focus of interest for this work is the cloud-free molecular atmosphere – that is, the interval of $d \in (3, 11]$ km.

The features of the MC data within the range of the cloud-free molecular atmosphere have a lot in common with those addressed in Sect. 3.1. In addition, the following properties are noteworthy. In the case of the space-borne lidar, the difference between values of $R_{MSto1}(d)$ and $R_{2to1}(d)$ is much larger compared to Fig. 2a and b; that is, the third and higher orders of scattering dominate. The escape effect does affect the lidar signals. At the same time, the rate of decrease of $R_{MSto1}(d)$ and $R_{2to1}(d)$ is much slower. The distance of 8 km from the cloud is not sufficient for reaching the conditions of the SS approximation. That result is of importance for practitioners who work with data of space-borne lidars because the estimated value of the apparent optical thickness depends on the distance from the cloud far edge chosen as the reference point $h_2$ (see details in Sect. 4.1.3).

As in the cases of the ground-based lidar, the EM is able to simulate MS contributions within the cloud layers well. We employed the same approach to obtain values of the ratio $P_n(\pi, h)/P_1(\pi, h)$, $n = 2, \dots, 5$, as in Sect. 3.1 (see Appendix A for details). As in the cases of the ground-based lidar, the fittings of the water cloud case were obtained using values of the ratios $P_n(\pi, h)/P_1(\pi, h)$, $n = 2, \dots, 5$, about 0.5 or lower. The behaviours of $R_{MSto1}(d)$ and $R_{2to1}(d)$ are correct within the range of the cloud-free molecular atmosphere, whereas the MS contribution is slightly overestimated. It is seen that the EM is not able to reproduce the amplitude of the stepwise jump of $R_{MSto1}(d)$ immediately beyond the cloud far edge. Thus, the stepwise jump in those cases is not only due to the stepwise jump in phase-function properties for angles close to 180°. We can suggest that the range of $d \in (3, 3.1]$ km is somewhat affected by the pulse stretching.

### 4.1.2 Non-uniform beam filling

Another problem associated with a reference range taken below the cloud base consists of the following. Space-borne lidars such as CALIOP have a quite large laser footprint so that a cloud field can be horizontally heterogeneous within the footprint. That raises the well-known problem of the non-uniform beam filling (NUBF). Some effects of the NUBF on lidar data within cloud fields were addressed in other works (Szczap et al., 2021; Alkasem et al., 2017). The lidar signals from the cloud-free atmosphere below the cloud base are affected as well.

In order to understand some basic properties under 3D multiple-scattering conditions, we performed simulations for the case of a very simple 3D field. (All data of this subsection were obtained using the McRALI software – i.e. the Monte Carlo method.) A homogeneous cloud covers half of the field, when viewed from the top. The centre of the CALIOP laser footprint is exactly on the cloud border. In other words, half of the laser beam passes through the cloud, and the other half goes through the molecular atmosphere. The thickness of the cloud is 3 km, and the cloud top altitude is 11 km. The extinction coefficient of the particles of the cloud is $2 \text{ km}^{-1}$. To put it differently, we have chosen the value of the extinction coefficient so that the amount of particles within the volume bounded by the EFOV in the 3D case is exactly the same as in the case of the plane-parallel homogeneous cloud. The same is true if we consider the sampling volume. Therefore, the optical thickness of the 3D cloud, when averaged over the EFOV, is $3 = (0 + 6)/2$.

The NUBF effect is so high in such conditions that it has to be shown in terms of lidar signals. Figure 8 is complementary to Fig. 7; i.e. it shows the results based on the same MC simulations but as the lidar signals. Figure 8a and b correspond to the water and cirrus cloud, respectively. The black and red lines are lidar signals (corrected for the offset and instrumental factors) obtained with the Monte Carlo method in MS conditions and the SS approximation, respectively. We recall that the black and red lines were obtained for the case of the homogeneous plane-parallel cloud with an optical thickness of $\tau_p = 3.0$.

The green lines in Fig. 8 are the lidar signals computed in the case of the 3D cloud field in MS conditions. Within the in-cloud altitude range of $H \in [8, 11]$ km, they are lower than the MS signals from the corresponding homogeneous cloud (the black lines), which is in total agreement with the theory (see Chap. 3.1 in Alkasem et al., 2017). As for the range below the cloud base, it is seen that the NUBF makes the situation much more aggravated. The lidar signals (the green curves) are around 200 times higher compared to the SS approximation for the homogeneous plane-parallel cloud (the red curves). That is the direct consequence of the fact that in the 3D case, half of the laser beam passes through the molecular atmosphere, whereas the other half passes through the cloud with $2\tau_p = 6$.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Monte Carlo simulations of multiple-scattering, $R_{MSto1}$ (black points), and double-scattering, $R_{2to1}$ (red points), relative contributions to lidar signals. (a) Water cloud. (b) Jet-stream cirrus. Eloranta model simulations are shown by the green (MS) and the blue (DS) curves.</figcaption></figure>

butions $R_{MSto1}(d)$ and $R_{2to1}(d)$ computed using the EM are shown in Fig. 7 by the green and blue curves, respectively. The MC in-cloud data – i.e. $d \in [0, 3]$ km – were reported and discussed in the work by Shcherbakov et al. (2022). The focus of interest for this work is the cloud-free molecular atmosphere – that is, the interval of $d \in (3, 11]$ km.

The features of the MC data within the range of the cloud-free molecular atmosphere have a lot in common with those addressed in Sect. 3.1. In addition, the following properties are noteworthy. In the case of the space-borne lidar, the difference between values of $R_{MSto1}(d)$ and $R_{2to1}(d)$ is much larger compared to Fig. 2a and b; that is, the third and higher orders of scattering dominate. The escape effect does affect the lidar signals. At the same time, the rate of decrease of $R_{MSto1}(d)$ and $R_{2to1}(d)$ is much slower. The distance of 8 km from the cloud is not sufficient for reaching the conditions of the SS approximation. That result is of importance for practitioners who work with data of space-borne lidars because the estimated value of the apparent optical thickness depends on the distance from the cloud far edge chosen as the reference point $h_2$ (see details in Sect. 4.1.3).

As in the cases of the ground-based lidar, the EM is able to simulate MS contributions within the cloud layers well. We employed the same approach to obtain values of the ratio $P_n(\pi, h)/P_1(\pi, h)$, $n = 2, \dots, 5$, as in Sect. 3.1 (see Appendix A for details). As in the cases of the ground-based lidar, the fittings of the water cloud case were obtained using values of the ratios $P_n(\pi, h)/P_1(\pi, h)$, $n = 2, \dots, 5$, about 0.5 or lower. The behaviours of $R_{MSto1}(d)$ and $R_{2to1}(d)$ are correct within the range of the cloud-free molecular atmosphere, whereas the MS contribution is slightly overestimated. It is seen that the EM is not able to reproduce the amplitude of the stepwise jump of $R_{MSto1}(d)$ immediately beyond the cloud far edge. Thus, the stepwise jump in those cases is not only due to the stepwise jump in phase-function properties for angles close to 180°. We can suggest that the range of $d \in (3, 3.1]$ km is somewhat affected by the pulse stretching.

### 4.1.2 Non-uniform beam filling

Another problem associated with a reference range taken below the cloud base consists of the following. Space-borne lidars such as CALIOP have a quite large laser footprint so that a cloud field can be horizontally heterogeneous within the footprint. That raises the well-known problem of the non-uniform beam filling (NUBF). Some effects of the NUBF on lidar data within cloud fields were addressed in other works (Szczap et al., 2021; Alkasem et al., 2017). The lidar signals from the cloud-free atmosphere below the cloud base are affected as well.

In order to understand some basic properties under 3D multiple-scattering conditions, we performed simulations for the case of a very simple 3D field. (All data of this subsection were obtained using the McRALI software – i.e. the Monte Carlo method.) A homogeneous cloud covers half of the field, when viewed from the top. The centre of the CALIOP laser footprint is exactly on the cloud border. In other words, half of the laser beam passes through the cloud, and the other half goes through the molecular atmosphere. The thickness of the cloud is 3 km, and the cloud top altitude is 11 km. The extinction coefficient of the particles of the cloud is $2 \text{ km}^{-1}$. To put it differently, we have chosen the value of the extinction coefficient so that the amount of particles within the volume bounded by the EFOV in the 3D case is exactly the same as in the case of the plane-parallel homogeneous cloud. The same is true if we consider the sampling volume. Therefore, the optical thickness of the 3D cloud, when averaged over the EFOV, is $3 = (0 + 6)/2$.

The NUBF effect is so high in such conditions that it has to be shown in terms of lidar signals. Figure 8 is complementary to Fig. 7; i.e. it shows the results based on the same MC simulations but as the lidar signals. Figure 8a and b correspond to the water and cirrus cloud, respectively. The black and red lines are lidar signals (corrected for the offset and instrumental factors) obtained with the Monte Carlo method in MS conditions and the SS approximation, respectively. We recall that the black and red lines were obtained for the case of the homogeneous plane-parallel cloud with an optical thickness of $\tau_p = 3.0$.

The green lines in Fig. 8 are the lidar signals computed in the case of the 3D cloud field in MS conditions. Within the in-cloud altitude range of $H \in [8, 11]$ km, they are lower than the MS signals from the corresponding homogeneous cloud (the black lines), which is in total agreement with the theory (see Chap. 3.1 in Alkasem et al., 2017). As for the range below the cloud base, it is seen that the NUBF makes the situation much more aggravated. The lidar signals (the green curves) are around 200 times higher compared to the SS approximation for the homogeneous plane-parallel cloud (the red curves). That is the direct consequence of the fact that in the 3D case, half of the laser beam passes through the molecular atmosphere, whereas the other half passes through the cloud with $2\tau_p = 6$.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Monte Carlo simulations of lidar signals. Homogeneous plane-parallel cloud: black curves correspond to multiple-scattering and red curves correspond to single-scattering conditions; green curves correspond to multiple scattering in the 3D cloud field. (a) Water cloud. (b) Jet-stream cirrus.</figcaption></figure>

### 4.1.3 Practical aspects

Knowing that the value of the extinction coefficient of $ε_p = 1.0 \text{ km}^{-1}$ is not typical of cirrus clouds, we provide practitioners with estimations of the escape and NUBF effects for lower values of $ε_p$ of the jet-stream cirrus. The black points in Fig. 9 show the values of the MS relative contribution, $R_{MSTo1}(d)$ (see Eq. 3), in the case of the homogeneous plane-parallel cloud. The MC simulations were performed for the same geometry of plane-parallel homogeneous cloud as in the case of Fig. 7b. The only difference consists of the values of the extinction coefficient, which were $ε_p = 0.06$ and $ε_p = 0.2 \text{ km}^{-1}$ in Fig. 9a and b, respectively. It is seen that the lidar data are affected by MS within the cloud ($d \in [0, 3] \text{ km}$) and below the cloud base ($d \in (3, 11] \text{ km}$).

NUBF effects are shown by the red points in Fig. 9. As previously, (i) a homogeneous cloud covers half of the field, when viewed from the top; (ii) the centre of the CALIOP laser footprint is exactly on the cloud border; and (iii) the value of $ε_p$ is a factor of 2 higher than the extinction coefficient of the corresponding plane-parallel cloud. The values of the extinction coefficient of the 3D cloud, when averaged over the EFOV, are $0.06 = (0+0.12)/2$ and $0.2 = (0+0.4)/2$ in the cases of Fig. 9a and b, respectively. The relative difference

$$R_{\text{NUBFto1}}(h) = [S_{\text{NUBF}}(h) - S_1(h)]/S_1(h) \quad (15)$$

is taken to be a measure of the NUBF effects, where $S_{\text{NUBF}}(h)$ is the MS lidar signal simulated using the Monte Carlo method in the case of the 3D cloud field and $S_1(h)$ is the lidar signal computed in the case of the correspond-

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Monte Carlo simulations of multiple-scattering effects (jet-stream cirrus). Black points are the relative contributions, $R_{MSTo1}$, to lidar signals from a homogeneous plane-parallel cloud; red points are the relative contributions, $R_{NUBFto1}$, of the MS and the NUBF to lidar signals from a 3D cloud. The extinction coefficients are (a) $0.06 \text{ km}^{-1}$ and (b) $0.2 \text{ km}^{-1}$.</figcaption></figure>

<table><caption>**Table 3.** Estimated values of the apparent optical thickness, $\eta \cdot \tau_p(h_b, h_{\text{end}})$. The real optical thickness of 3D clouds is the value averaged over the EFOV.</caption><thead><tr><th>MC data case (figure)</th><th>Apparent optical thickness</th><th>Real optical thickness</th></tr></thead><tbody><tr><td>Distance from cloud (km)</td><td>0.2</td><td>7.98</td></tr><tr><td>Water cloud (Fig. 7a)</td><td>1.91</td><td>2.81</td><td>3.0</td></tr><tr><td>Jet-stream cirrus (Fig. 7b)</td><td>1.60</td><td>2.36</td><td>3.0</td></tr><tr><td>Water cloud (3D) (Fig. 8a)</td><td>0.32</td><td>0.35</td><td>3.0</td></tr><tr><td>Jet-stream cirrus (3D) (Fig. 8b)</td><td>0.31</td><td>0.34</td><td>3.0</td></tr><tr><td>Jet-stream cirrus (Fig. 9a)</td><td>0.10</td><td>0.14</td><td>0.18</td></tr><tr><td>Jet-stream cirrus (3D) (Fig. 9a)</td><td>0.09</td><td>0.12</td><td>0.18</td></tr><tr><td>Jet-stream cirrus (Fig. 9b)</td><td>0.34</td><td>0.48</td><td>0.60</td></tr><tr><td>Jet-stream cirrus (3D) (Fig. 9b)</td><td>0.21</td><td>0.26</td><td>0.60</td></tr></tbody></table>

ing plane-parallel cloud under the SS approximation. In other words, the parameter $R_{\text{NUBFto1}}$ is devoted to showing the relative contribution of the MS and the NUBF taken together.

At another time, it is seen that the MS effect is somewhat weakened by the NUBF within the cloud range ($d \in [0, 3] \text{ km}$). When the optical thickness within the cirrus cloud is quite low, MS lidar signals from a 3D cloud are fairly close to the SS lidar signals from the corresponding plane-parallel homogeneous cloud (see Figs. 8b and 9). That property is not valid in the case of the water cloud (see Fig. 8a). It can be hypothesized that more photons escape through the cloud side when the forward-scattering peak of the particles is larger.

As for lidar signals within the range of the free atmosphere below the cloud base ($d \in (3, 11] \text{ km}$), they are always affected either by the multiple scattering or by the non-uniform beam filling. There is no interval where the MS effect is lower than 5%. Also, the NUBF effect on lidar signals is always higher than the MS effect in the case of plane-parallel cloud. We recall that our comparison is made on the condition that the number of particles within the volume bounded by the EFOV in the 3D case is exactly the same as in the case of the plane-parallel homogeneous cloud.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Monte Carlo simulations of lidar signals. Homogeneous plane-parallel cloud: black curves correspond to multiple-scattering and red curves correspond to single-scattering conditions; green curves correspond to multiple scattering in the 3D cloud field. (a) Water cloud. (b) Jet-stream cirrus.</figcaption></figure>

### 4.1.3 Practical aspects

Knowing that the value of the extinction coefficient of $ε_p = 1.0 \text{ km}^{-1}$ is not typical of cirrus clouds, we provide practitioners with estimations of the escape and NUBF effects for lower values of $ε_p$ of the jet-stream cirrus. The black points in Fig. 9 show the values of the MS relative contribution, $R_{MSTo1}(d)$ (see Eq. 3), in the case of the homogeneous plane-parallel cloud. The MC simulations were performed for the same geometry of plane-parallel homogeneous cloud as in the case of Fig. 7b. The only difference consists of the values of the extinction coefficient, which were $ε_p = 0.06$ and $ε_p = 0.2 \text{ km}^{-1}$ in Fig. 9a and b, respectively. It is seen that the lidar data are affected by MS within the cloud ($d \in [0, 3] \text{ km}$) and below the cloud base ($d \in (3, 11] \text{ km}$).

NUBF effects are shown by the red points in Fig. 9. As previously, (i) a homogeneous cloud covers half of the field, when viewed from the top; (ii) the centre of the CALIOP laser footprint is exactly on the cloud border; and (iii) the value of $ε_p$ is a factor of 2 higher than the extinction coefficient of the corresponding plane-parallel cloud. The values of the extinction coefficient of the 3D cloud, when averaged over the EFOV, are $0.06 = (0+0.12)/2$ and $0.2 = (0+0.4)/2$ in the cases of Fig. 9a and b, respectively. The relative difference

$$R_{\text{NUBFto1}}(h) = [S_{\text{NUBF}}(h) - S_{1}(h)]/S_{1}(h) \quad (15)$$

is taken to be a measure of the NUBF effects, where $S_{\text{NUBF}}(h)$ is the MS lidar signal simulated using the Monte Carlo method in the case of the 3D cloud field and $S_1(h)$ is the lidar signal computed in the case of the correspond-

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Monte Carlo simulations of multiple-scattering effects (jet-stream cirrus). Black points are the relative contributions, $R_{MSTo1}$, to lidar signals from a homogeneous plane-parallel cloud; red points are the relative contributions, $R_{NUBFto1}$, of the MS and the NUBF to lidar signals from a 3D cloud. The extinction coefficients are (a) $0.06 \text{ km}^{-1}$ and (b) $0.2 \text{ km}^{-1}$.</figcaption></figure>

<table><caption>**Table 3.** Estimated values of the apparent optical thickness, $\eta \cdot \tau_p(h_b, h_{\text{end}})$. The real optical thickness of 3D clouds is the value averaged over the EFOV.</caption><thead><tr><th>MC data case (figure)</th><th>Apparent optical thickness</th><th>Real optical thickness</th></tr></thead><tbody><tr><td>Distance from cloud (km)</td><td>0.2</td><td>7.98</td></tr><tr><td>Water cloud (Fig. 7a)</td><td>1.91</td><td>2.81</td><td>3.0</td></tr><tr><td>Jet-stream cirrus (Fig. 7b)</td><td>1.60</td><td>2.36</td><td>3.0</td></tr><tr><td>Water cloud (3D) (Fig. 8a)</td><td>0.32</td><td>0.35</td><td>3.0</td></tr><tr><td>Jet-stream cirrus (3D) (Fig. 8b)</td><td>0.31</td><td>0.34</td><td>3.0</td></tr><tr><td>Jet-stream cirrus (Fig. 9a)</td><td>0.10</td><td>0.14</td><td>0.18</td></tr><tr><td>Jet-stream cirrus (3D) (Fig. 9a)</td><td>0.09</td><td>0.12</td><td>0.18</td></tr><tr><td>Jet-stream cirrus (Fig. 9b)</td><td>0.34</td><td>0.48</td><td>0.60</td></tr><tr><td>Jet-stream cirrus (3D) (Fig. 9b)</td><td>0.21</td><td>0.26</td><td>0.60</td></tr></tbody></table>

ing plane-parallel cloud under the SS approximation. In other words, the parameter $R_{\text{NUBFto1}}$ is devoted to showing the relative contribution of the MS and the NUBF taken together.

At another time, it is seen that the MS effect is somewhat weakened by the NUBF within the cloud range ($d \in [0, 3] \text{ km}$). When the optical thickness within the cirrus cloud is quite low, MS lidar signals from a 3D cloud are fairly close to the SS lidar signals from the corresponding plane-parallel homogeneous cloud (see Figs. 8b and 9). That property is not valid in the case of the water cloud (see Fig. 8a). It can be hypothesized that more photons escape through the cloud side when the forward-scattering peak of the particles is larger.

As for lidar signals within the range of the free atmosphere below the cloud base ($d \in (3, 11] \text{ km}$), they are always affected either by the multiple scattering or by the non-uniform beam filling. There is no interval where the MS effect is lower than 5%. Also, the NUBF effect on lidar signals is always higher than the MS effect in the case of plane-parallel cloud. We recall that our comparison is made on the condition that the number of particles within the volume bounded by the EFOV in the 3D case is exactly the same as in the case of the plane-parallel homogeneous cloud.
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Monte Carlo simulations of multiple-scattering, $R_{MSto1}$ (black points), and double-scattering, $R_{2to1}$ (red points), relative contributions to lidar signals. (a) Water cloud. (b) Jet-stream cirrus. Eloranta model simulations are shown by the green (MS) and the blue (DS) curves.</figcaption></figure>

The escape effect is noteworthy as well. It is seen that $R_{MSto1}(d)$ and $R_{NUBFto1}(d)$ decrease with the increasing distance from the cloud base in all four cases of Fig. 9. Examples of its consequences in terms of the apparent optical thickness of the cloud, $\eta \cdot \tau_p(h_b, h_{end})$, are given in Table 3. The data were computed using Eq. (13), with two values of $h_2$ assigned within the interval of $d \in (3, 11]$ km (with an altitude of $H \in (8, 0]$ km). The values 0.2 and 7.98 km mean the distance from the cloud far edge within the cloud-free molecular atmosphere; they correspond to $d = 3.2$ ($H = 7.8$) and $d = 10.98$ ($H = 0.02$) km, respectively. It is seen that the estimated values of $\eta \cdot \tau_p(h_b, h_{end})$ are always lower than the values of the real optical thickness; they are much lower in the cases of the 3D cloud (NUBF) compared to the plane-parallel cloud, and they strongly depend on the assigned value of $h_2$.

As in the case of the ground-based lidar, we can conclude that Eq. (5) does not take into account the escape effect. Therefore, it is only an approximate model of MS signals from space-based lidars.

## 4.2 Two-layered cloud

The two-layered cloud is the same as in Sect. 3.2 – that is, two homogeneous layers are at altitudes from 8 to 9 km and from 10 to 11 km. Every layer has an optical thickness of $\tau_p = 1.0$; the total optical thickness is $\tau_p = 2.0$. The black and red points in Fig. 10 show the results of our MC simulations reported in terms of the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$, respectively. The parameter $d$ is the distance measured from the cloud near edge – i.e. $d = 0$ km corresponds to an altitude of $H = 11$ km. The cloud layers are within the ranges of $d \in [0, 1]$ and $d \in [2, 3]$ km. The intervals of $d \in (1, 2)$ and $d \in (3, 11]$ km (with altitudes, $H$, from 10 to 9 km and from 8 to 0 km) are the cloud-free molecular atmosphere. Figure 10a and b correspond to the water and cirrus cloud, respectively.

The main properties of our MC simulations in Fig. 10 are in agreement with the results of the work by Flesia and Starkov (1996); the distinctions are due to differences in the phase functions and the configuration characteristics. In ad-

dition, the vertical resolution in our work is 7.5 times better, and, therefore, the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$ are seen in more detail.

The general behaviour of MC data around $d = 1$ km and within the range of $d \ge 3$ km in Fig. 10 is the same as it is within the range of $d \ge 3$ km in Fig. 7. Thus, we focus attention only on the interval around $d = 2$ km – that is, on the near edge of the second cloud layer. The difference between the cases of the ground-based and the space-borne lidars is better seen from comparison of Figs. 5a and 10a, i.e. the water clouds. In the case of the space-borne lidar, the sampling volume is so large (due to the large receiver footprint) that the majority of forward-scattered photons remain within it. This leads to the value $R_{MSto1}(d = 2.0$ km), which is almost the same as $R_{MSto1}(d = 1.0$ km) (see Fig. 10a). In contrast, the sampling volume, i.e. the footprint, is narrow in the case of the ground-based lidar. The majority of forward-scattered photons escape it when they are going through the cloud-free molecular atmosphere within the interval $d \in (1, 2)$ km. This leads to the value $R_{MSto1}(d = 2.0$ km) close to zero (see Fig. 5a). Generally, the same features are observed in Figs. 5b and 10b. They are less pronounced because the forward-scattering peak of the phase function of cirrus cloud is much narrower.

The relative contributions $R_{MSto1}(d)$ and $R_{2to1}(d)$ computed using the EM are shown in Fig. 10 by green and blue curves, respectively. The simulations were performed with the same values of the ratio $P_n(\pi, h)/P_1(\pi, h)$, $n = 2, \dots, 5$, as in Sect. 4.1.1. Generally, the EM curves follow the MC results well. It is seen another time that the EM slightly overestimates the MS contribution in the cloud-free intervals and it is not able to reproduce the effect related to the pulse stretching.

## 5 Conclusions

We performed Monte Carlo simulations of single-wavelength lidar signals from multi-layered clouds with special attention focused on features of the multiple-scattering (MS) effect in regions of the cloud-free molecular atmosphere, i.e. between layers or outside a cloud system. Despite the fact that the strength of lidar signals from molecular atmosphere is much lower compared to the in-cloud intervals, studies of MS effects in such regions are of interest from scientific and practical points of view. The results of this work are shown and discussed in terms of relative contributions of multiple scattering, that is, multiple-to-single-scattering lidar return ratios $R_{MSto1}(d)$. That provides the possibility of accentuating the visibility of MS effects.

The scattered photons emerging from a cloud do affect lidar signals received from the intervals of the cloud-free molecular atmosphere. The MS effect is rather high within the region that is close to the far edge of a cloud. Those high values of the ratios $R_{MSto1}(d)$ are due to the features
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Monte Carlo simulations of multiple-scattering, $R_{MSto1}$ (black points), and double-scattering, $R_{2to1}$ (red points), relative contributions to lidar signals. (a) Water cloud. (b) Jet-stream cirrus. Eloranta model simulations are shown by the green (MS) and the blue (DS) curves.</figcaption></figure>

The escape effect is noteworthy as well. It is seen that $R_{MSto1}(d)$ and $R_{NUBFto1}(d)$ decrease with the increasing distance from the cloud base in all four cases of Fig. 9. Examples of its consequences in terms of the apparent optical thickness of the cloud, $\eta \cdot \tau_p(h_b, h_{end})$, are given in Table 3. The data were computed using Eq. (13), with two values of $h_2$ assigned within the interval of $d \in (3, 11]$ km (with an altitude of $H \in (8, 0]$ km). The values 0.2 and 7.98 km mean the distance from the cloud far edge within the cloud-free molecular atmosphere; they correspond to $d = 3.2$ ($H = 7.8$) and $d = 10.98$ ($H = 0.02$) km, respectively. It is seen that the estimated values of $\eta \cdot \tau_p(h_b, h_{end})$ are always lower than the values of the real optical thickness; they are much lower in the cases of the 3D cloud (NUBF) compared to the plane-parallel cloud, and they strongly depend on the assigned value of $h_2$.

As in the case of the ground-based lidar, we can conclude that Eq. (5) does not take into account the escape effect. Therefore, it is only an approximate model of MS signals from space-based lidars.

## 4.2 Two-layered cloud

The two-layered cloud is the same as in Sect. 3.2 – that is, two homogeneous layers are at altitudes from 8 to 9 km and from 10 to 11 km. Every layer has an optical thickness of $\tau_p = 1.0$; the total optical thickness is $\tau_p = 2.0$. The black and red points in Fig. 10 show the results of our MC simulations reported in terms of the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$, respectively. The parameter $d$ is the distance measured from the cloud near edge – i.e. $d = 0$ km corresponds to an altitude of $H = 11$ km. The cloud layers are within the ranges of $d \in [0, 1]$ and $d \in [2, 3]$ km. The intervals of $d \in (1, 2)$ and $d \in (3, 11]$ km (with altitudes, $H$, from 10 to 9 km and from 8 to 0 km) are the cloud-free molecular atmosphere. Figure 10a and b correspond to the water and cirrus cloud, respectively.

The main properties of our MC simulations in Fig. 10 are in agreement with the results of the work by Flesia and Starkov (1996); the distinctions are due to differences in the phase functions and the configuration characteristics. In ad-

dition, the vertical resolution in our work is 7.5 times better, and, therefore, the ratios $R_{MSto1}(d)$ and $R_{2to1}(d)$ are seen in more detail.

The general behaviour of MC data around $d = 1$ km and within the range of $d \ge 3$ km in Fig. 10 is the same as it is within the range of $d \ge 3$ km in Fig. 7. Thus, we focus attention only on the interval around $d = 2$ km – that is, on the near edge of the second cloud layer. The difference between the cases of the ground-based and the space-borne lidars is better seen from comparison of Figs. 5a and 10a, i.e. the water clouds. In the case of the space-borne lidar, the sampling volume is so large (due to the large receiver footprint) that the majority of forward-scattered photons remain within it. This leads to the value $R_{MSto1}(d = 2.0$ km), which is almost the same as $R_{MSto1}(d = 1.0$ km) (see Fig. 10a). In contrast, the sampling volume, i.e. the footprint, is narrow in the case of the ground-based lidar. The majority of forward-scattered photons escape it when they are going through the cloud-free molecular atmosphere within the interval $d \in (1, 2)$ km. This leads to the value $R_{MSto1}(d = 2.0$ km) close to zero (see Fig. 5a). Generally, the same features are observed in Figs. 5b and 10b. They are less pronounced because the forward-scattering peak of the phase function of cirrus cloud is much narrower.

The relative contributions $R_{MSto1}(d)$ and $R_{2to1}(d)$ computed using the EM are shown in Fig. 10 by green and blue curves, respectively. The simulations were performed with the same values of the ratio $P_n(\pi, h)/P_1(\pi, h)$, $n = 2, \dots, 5$, as in Sect. 4.1.1. Generally, the EM curves follow the MC results well. It is seen another time that the EM slightly overestimates the MS contribution in the cloud-free intervals and it is not able to reproduce the effect related to the pulse stretching.

## 5 Conclusions

We performed Monte Carlo simulations of single-wavelength lidar signals from multi-layered clouds with special attention focused on features of the multiple-scattering (MS) effect in regions of the cloud-free molecular atmosphere, i.e. between layers or outside a cloud system. Despite the fact that the strength of lidar signals from molecular atmosphere is much lower compared to the in-cloud intervals, studies of MS effects in such regions are of interest from scientific and practical points of view. The results of this work are shown and discussed in terms of relative contributions of multiple scattering, that is, multiple-to-single-scattering lidar return ratios $R_{MSto1}(d)$. That provides the possibility of accentuating the visibility of MS effects.

The scattered photons emerging from a cloud do affect lidar signals received from the intervals of the cloud-free molecular atmosphere. The MS effect is rather high within the region that is close to the far edge of a cloud. Those high values of the ratios $R_{MSto1}(d)$ are due to the features
of the molecular backscattering, i.e. to the fact that the ratios $P_n(\pi, h)/P_1(\pi, h)$ are “equal to 1.0 due to the broad nature of the molecular phase function near the backscatter direction” (Whiteman et al., 2001; Eloranta, 1998). In the cases of the space-borne lidar, the additional MS contribution is due to the pulse stretching.

The MS effect on lidar signals decreases with the increasing distance from the cloud far edge; i.e. the ratio $R_{MSto1}(d)$ tends to zero. The decrease is the direct consequence of the fact that the forward peak of particle phase functions is much larger than the receiver field of view. Therefore, the photons scattered within the forward peak escape the sampling volume formed by the RFOV – i.e. the escape effect. The escape effect is an inherent part of MS properties within the free atmosphere beyond the cloud far edge. That property is in direct contradiction with Eq. (5). Consequently, Eq. (5) is an approximate model of lidar signals under MS conditions.

The two-way transmittance method (see e.g. Young and Vaughan, 2009; Giannakaki et al., 2007) is based on Eq. (5) and used to deduce the values of the apparent optical thickness of clouds. In view of the results of this work, it is advisable at least to choose the reference point always at the same distance from the cloud far edge when estimating the apparent optical thickness of clouds.

In the cases of the ground-based lidar, the MS contribution is lower than 5% within the regions of the cloud-free molecular atmosphere, with the distance from the cloud far edge of about 1 km or higher. Therefore, it is safe to say that practitioners can use those regions as a reference and estimate the real optical thickness of clouds (see e.g. Giannakaki et al., 2007, and references therein) if the EFOV and the RFOV are not very large. In the cases of the space-borne lidar, the rate of decrease of the MS contribution is so slow that the threshold of 5% can hardly be reached.

Using an example of a very simple 3D field, we demonstrated that the effect of non-uniform beam filling can be extremely strong in the case of a space-borne lidar. It is so strong that, in our opinion, practitioners should employ, with proper precautions, lidar data from regions below the cloud base, when treating data from a space-borne lidar. At the same time, it should be underlined that the effects of the NUBF need further study, which will provide statistically significant results.

In the case of two-layered cloud, the distance of 1 km is sufficiently large so that the scattered photons that emerge from the first layer do not affect signals from the second layer when we are dealing with the ground-based lidar. In contrast, signals from the near edge of the second cloud layer are severely affected by the photons emerging from the first layer in the case of a space-borne lidar.

We evaluate the Eloranta model (EM) (Eloranta, 1998) in extreme conditions and show its good performance in the cases of ground-based and space-borne (CALIOP) lidars. When the extinction coefficient is about $1.0 \text{ km}^{-1}$ or lower, and the EFOV and the RFOV are quite narrow, five orders of

scattering are sufficient for obtaining satisfying accuracy of the simulations. At the same time, we revealed the shortcoming that affects practical applications of the EM. Namely, values of the key parameters, i.e. of the ratio $P_n(\pi, h)/P_1(\pi, h)$ of phase functions in the backscatter direction for the nth-order-scattered photon and a singly scattered photon, depend not only on the particle phase function but also on the distance from a lidar to the cloud and the receiver field of view. Values of the ratio $P_n(\pi, h)/P_1(\pi, h)$ vary within a quite large range. Therefore, the multiple-scattering contribution to lidar signals can be largely overestimated or underestimated if erroneous values of the ratios are assigned to the EM. That problem can be circumvented using Monte Carlo simulations or the empirical model (Shcherbakov et al., 2022) to calibrate the ratio $P_n(\pi, h)/P_1(\pi, h)$.

## Appendix A: Input parameters of the Eloranta model in homogeneous-cloud conditions

### A1 Eloranta model

We have chosen to evaluate the Eloranta model (EM) (Eloranta, 1998) due to its following attractive features. The input parameters have a clear physical meaning. The contribution of each of the *n* orders of scattering can be computed separately. The corresponding code can be developed without much difficulty because multiple integrals are the core of the EM. The good performance of the Eloranta model (EM) (Eloranta, 1998) in homogeneous-cloud conditions has already been reported in the literature (see e.g. Eloranta, 1998; Donovan and Van Lammeren, 2001). In our opinion, our results above demonstrate the reliability of the EM in extreme conditions – i.e. when the extinction coefficient has undergone a stepwise jump. The objective of this appendix is to reveal some significant features of the EM’s input parameters that are related to the particle phase function.

We developed two versions of the code based on the EM. The first version corresponds to Eq. (8) of Eloranta (1998). In order to avoid ambiguity, we rewrite that equation in a way that all functions, including $\varepsilon(x)$, $\gamma(x)$, and $\Theta_s(x)$, are assigned in the coordinate system where the lidar is at $h=0$.
of the molecular backscattering, i.e. to the fact that the ratios $P_n(\pi, h)/P_1(\pi, h)$ are “equal to 1.0 due to the broad nature of the molecular phase function near the backscatter direction” (Whiteman et al., 2001; Eloranta, 1998). In the cases of the space-borne lidar, the additional MS contribution is due to the pulse stretching.

The MS effect on lidar signals decreases with the increasing distance from the cloud far edge; i.e. the ratio $R_{MSto1}(d)$ tends to zero. The decrease is the direct consequence of the fact that the forward peak of particle phase functions is much larger than the receiver field of view. Therefore, the photons scattered within the forward peak escape the sampling volume formed by the RFOV – i.e. the escape effect. The escape effect is an inherent part of MS properties within the free atmosphere beyond the cloud far edge. That property is in direct contradiction with Eq. (5). Consequently, Eq. (5) is an approximate model of lidar signals under MS conditions.

The two-way transmittance method (see e.g. Young and Vaughan, 2009; Giannakaki et al., 2007) is based on Eq. (5) and used to deduce the values of the apparent optical thickness of clouds. In view of the results of this work, it is advisable at least to choose the reference point always at the same distance from the cloud far edge when estimating the apparent optical thickness of clouds.

In the cases of the ground-based lidar, the MS contribution is lower than 5% within the regions of the cloud-free molecular atmosphere, with the distance from the cloud far edge of about 1 km or higher. Therefore, it is safe to say that practitioners can use those regions as a reference and estimate the real optical thickness of clouds (see e.g. Giannakaki et al., 2007, and references therein) if the EFOV and the RFOV are not very large. In the cases of the space-borne lidar, the rate of decrease of the MS contribution is so slow that the threshold of 5% can hardly be reached.

Using an example of a very simple 3D field, we demonstrated that the effect of non-uniform beam filling can be extremely strong in the case of a space-borne lidar. It is so strong that, in our opinion, practitioners should employ, with proper precautions, lidar data from regions below the cloud base, when treating data from a space-borne lidar. At the same time, it should be underlined that the effects of the NUBF need further study, which will provide statistically significant results.

In the case of two-layered cloud, the distance of 1 km is sufficiently large so that the scattered photons that emerge from the first layer do not affect signals from the second layer when we are dealing with the ground-based lidar. In contrast, signals from the near edge of the second cloud layer are severely affected by the photons emerging from the first layer in the case of a space-borne lidar.

We evaluate the Eloranta model (EM) (Eloranta, 1998) in extreme conditions and show its good performance in the cases of ground-based and space-borne (CALIOP) lidars. When the extinction coefficient is about $1.0 \text{ km}^{-1}$ or lower, and the EFOV and the RFOV are quite narrow, five orders of

scattering are sufficient for obtaining satisfying accuracy of the simulations. At the same time, we revealed the shortcoming that affects practical applications of the EM. Namely, values of the key parameters, i.e. of the ratio $P_n(\pi, h)/P_1(\pi, h)$ of phase functions in the backscatter direction for the nth-order-scattered photon and a singly scattered photon, depend not only on the particle phase function but also on the distance from a lidar to the cloud and the receiver field of view. Values of the ratio $P_n(\pi, h)/P_1(\pi, h)$ vary within a quite large range. Therefore, the multiple-scattering contribution to lidar signals can be largely overestimated or underestimated if erroneous values of the ratios are assigned to the EM. That problem can be circumvented using Monte Carlo simulations or the empirical model (Shcherbakov et al., 2022) to calibrate the ratio $P_n(\pi, h)/P_1(\pi, h)$.

## Appendix A: Input parameters of the Eloranta model in homogeneous-cloud conditions

### A1 Eloranta model

We have chosen to evaluate the Eloranta model (EM) (Eloranta, 1998) due to its following attractive features. The input parameters have a clear physical meaning. The contribution of each of the *n* orders of scattering can be computed separately. The corresponding code can be developed without much difficulty because multiple integrals are the core of the EM. The good performance of the Eloranta model (EM) (Eloranta, 1998) in homogeneous-cloud conditions has already been reported in the literature (see e.g. Eloranta, 1998; Donovan and Van Lammeren, 2001). In our opinion, our results above demonstrate the reliability of the EM in extreme conditions – i.e. when the extinction coefficient has undergone a stepwise jump. The objective of this appendix is to reveal some significant features of the EM’s input parameters that are related to the particle phase function.

We developed two versions of the code based on the EM. The first version corresponds to Eq. (8) of Eloranta (1998). In order to avoid ambiguity, we rewrite that equation in a way that all functions, including $\varepsilon(x)$, $\gamma(x)$, and $\Theta_s(x)$, are assigned in the coordinate system where the lidar is at $h=0$.
and *h* is the distance from the lidar:

$$
\begin{align}
R_{\text{ntol}}(h) &= \frac{P_n(\pi, h)}{P_1(\pi, h)} \left[ 1 - \exp\left(-\frac{\rho_l^2}{\rho_l^2}\right) \right]^{-1} \cdot \left\{ \int_{-d}^{d} \right. \nonumber \\
&\qquad \left. \gamma(h - |x_1|) \cdot \varepsilon(h - |x_1|) \int_{x_1}^{d} \gamma(h - |x_2|) \cdot \varepsilon(h - |x_2|) \right. \nonumber \\
&\qquad \left. \dots \int_{x_{n-2}}^{d} \gamma(h - |x_{n-1}|) \cdot \varepsilon(h - |x_{n-1}|) \cdot [1.0 - \exp\left(\frac{\rho_l^2 h^2}{x_1^2 \Theta_s^2 (h - |x_1|) + x_2^2 \Theta_s^2 (h - |x_2|) + \dots + x_{n-1}^2 \Theta_s^2 (h - |x_{n-1}|) + \rho_l^2 h^2}\right)] \right\} \nonumber \\
&\qquad \times dx_{n-1} \dots dx_2 dx_1 \tag{A1}
\end{align}
$$

Here, $n \ge 2$ and $\varepsilon(h) = \varepsilon_p(h) + \varepsilon_m(h)$, where $\varepsilon_p(h)$ and $\varepsilon_m(h)$ are the particles and the molecular extinction coefficients, respectively. The function $\gamma(h)$ is the fraction of the energy in the forward peak of the phase function, $h_b$ is the distance from the lidar to the cloud near edge, and $d = h - h_b$. If calculations are performed only within a cloud, $d$ is just the cloud penetration depth. As usual, the notation $|x_i|$ means the absolute value of $x_i$. $\rho_l$ is the half angle of the receiver field of view, and $\rho_l$ is the half angle of the emitter divergence. $\Theta_s$ is the $1/e$ diffraction peak angular half width (Whiteman et al., 2001); in other words, $\Theta_s$ is the parameter that characterizes a Gaussian approximation for the forward-scattered peak. $\Theta_s$ can be estimated using the relationship (Eloranta, 1998; Hogan, 2006)

$$ \Theta_s = \lambda / (\pi r_G), \qquad (A2) $$

where $r_G$ is the equivalent-area radius of the particles size distribution and $\lambda$ is the wavelength.

Anther input parameter of the EM is the ratio $P_n(\pi, h)/P_1(\pi, h)$ of phase functions in the backscatter direction for the nth-order-scattered photon and a singly scattered photon (Whiteman et al., 2001). It is assumed that the backscattered phase function $P_n(\pi, h)$ for nth-order scattering is independent of the angle near 180° with a value which is a weighted average near the backscatter direction (Eloranta, 1998). The following is underlined in the work by Eloranta (1998): “for observations it will be necessary to use assumed values. For typical phase functions, $P_n(\pi, h)/P_1(\pi, h)$ is between 0.5 and 1”.

The second version of the EM code corresponds to Eq. (11) of Eloranta (1998), where the constant value of $\gamma(x) = 1/2$ was assumed using a reference to diffraction theory. We reproduce the equation for the nth order of scattering with the reformulation outlined in Eloranta (2000) and Whiteman et al. (2001) (see Eq. 13 in Eloranta, 2000) and

<table><caption>Table A1. Fitting values of the EM parameters.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Ground-based lidar</th><th colspan="2">Space-borne lidar</th></tr><tr><th>Water cloud</th><th>JS cirrus</th><th>Water cloud</th><th>JS cirrus</th></tr></thead><tbody><tr><td>Θ<sub>s</sub> (rad)</td><td>0.01882</td><td>0.0064</td><td>0.01882</td><td>0.0064</td></tr><tr><td>Θ<sub>s</sub> (°)</td><td>1.078</td><td>0.3667</td><td>1.078</td><td>0.3667</td></tr><tr><td>P<sub>2</sub> (π) / P<sub>1</sub> (π)</td><td>0.502</td><td>0.7721</td><td>0.4993</td><td>0.8211</td></tr><tr><td>P<sub>3</sub> (π) / P<sub>1</sub> (π)</td><td>0.395</td><td>0.66</td><td>0.5</td><td>0.99</td></tr><tr><td>P<sub>4</sub> (π) / P<sub>1</sub> (π)</td><td>0.33</td><td>0.51</td><td>0.4</td><td>0.99</td></tr><tr><td>P<sub>5</sub> (π) / P<sub>1</sub> (π)</td><td>0.26</td><td>0.553</td><td>0.35</td><td>0.99</td></tr></tbody></table>

using our notations:

$$
\begin{equation}
\begin{split}
R_{\text{ntol}}(h) = \frac{P_n(\pi, h)}{P_1(\pi, h)} & \left[ 1 - \exp\left(-\frac{\rho_l^2}{\rho_l^2}\right) \right]^{-1} \\
& \cdot \left\{ \frac{\tau^{n-1}(h_b, h)}{(n-1)!} - \int_{h_b}^{h} \varepsilon(x_1) \int_{x_1}^{h} \varepsilon(x_2) \cdots \int_{x_{n-2}}^{h} \varepsilon(x_{n-1}) \exp \right. \\
& \left. \left( \frac{\rho_l^2 h^2}{(h-x_1)^2 \Theta_s^2(x_1) + (h-x_2)^2 \Theta_s^2(x_2) + \cdots + (h-x_{n-1})^2 \Theta_s^2(x_{n-1}) + \rho_l^2 h^2} \right) \right\} \\
& \times dx_{n-1} \cdots dx_2 dx_1 \},
\end{split}
\tag{A3}
\end{equation}
$$

where $\tau(h_b, h) = \int_{h_b}^{h} \varepsilon(x) dx$ and other notations have the same meaning as in Eq. (A1).

We used Eqs. (A1) and (A3) only in cases when $h \ge h_b$; we have verified that the first version of the EM code (i.e. Eq. A1), gives the same results as Eq. (A3) if the constant value $\gamma(x) = 1/2$ is assigned in Eq. (A1). We tested our codes of the EM against the data available in the literature. All tests were performed with $\gamma(x) = 1/2$. We obtained perfect agreement with Figs. 6–8 of the work by Eloranta (1998) using $P_n(\pi, h)/P_1(\pi, h) = 0.75, n = 2, 3, 4$, and with Fig. 15b of the work by Donovan and Van Lamberen (2001) using $P_n(\pi, h)/P_1(\pi, h) = 1.0, n = 2, 3, 4$.

## A2 Input parameters

The results reported in Fig. A1 were obtained for the same configuration that was described above (see Sect. 3.1). That is, a ground-based lidar is at an altitude of $H = 0$ km; the distance to the cloud base is 8 km. The full RFOV is 1.0 mrad; the full EFOV is 0.14 mrad. The emitter wavelength, $\lambda$, is 0.532 µm. The single-layer homogeneous cloud is within the altitude range of $H \in [8, 11]$ km and has an optical thickness of $\tau_p = 3.0$; i.e. cloud particles have an extinction coefficient of $\varepsilon_p(h) = 1.0 \text{ km}^{-1}$.
<figure><img src="image_14.png" /></figure>

**Figure A1.** Monte Carlo simulations of multiple-scattering, $R_{MSto1}$ (black points); double-scattering $R_{2to1}$ (red points); and the third-order, $R_{3to1}$ (d) (blue points), and fourth-order, $R_{4to1}$ (d) (green points), respectively, relative contributions to lidar signals. (a, c) Water cloud. (b, d) Jet-stream cirrus. Eloranta model simulations are shown by green (MS), blue (DS), red (third-order), and black (fourth-order) curves.

The results of the Monte Carlo modelling (points) and the EM's simulations (lines) are shown in Fig. A1 only for the in-cloud range of $d \in [0, 3]$ km – i.e. the altitude interval of $H \in [8, 11]$ km. Figure A1a and c correspond to the water cloud and Fig. A1b and d to the cirrus cloud. The microphysical and optical parameters of cloud particles are described in Sect. 2.2 above. The results of each case are divided in two; Fig. A1a and c show the relative contribution of multiple, $R_{MSto1}(d)$ (black points and green line), and double, $R_{2to1}(d)$ (red points and blue line), scattering; Fig. A1b and d show the relative contribution of the third, $R_{3to1}(d)$ (blue points and red line), and the fourth, $R_{4to1}(d)$ (green points and black line), orders. In other words, the lower panels are complementary to the corresponding upper panels. (Note that each panel in Fig. A1 has its own scale of the y axis.)

The good agreement between the MC and EM data is evident. That agreement was obtained by adjusting the EM parameters in the following way. The distance, $h_b$, to the cloud near edge; the extinction coefficient, $\epsilon_p(d)$; the half angle, $\rho_t$, of the receiver field of view; and the half angle, $\rho_l$, of the emitter divergence were assigned according to the configuration used for the MC simulations. First, we found the values of $\Theta_s$ and $P_2(\pi)/P_1(\pi)$ by fitting the MC data on $R_{2to1}(d)$ with the Eloranta model using the ordinary least-squares approach. Then we found the values of $P_3(\pi)/P_1(\pi)$ and $P_4(\pi)/P_1(\pi)$ by fitting the MC data to $R_{3to1}(d)$ and $R_{4to1}(d)$, respectively. We have limited our EM simulations by five orders of scattering. Consequently, the ratio $P_5(\pi)/P_1(\pi)$ was computed to fit the remaining part of the total multiple scattering.

The obtained values of the parameters are shown in Table A1. The value of $\Theta_s$ is in good agreement with Eq. (A2) for the water cloud; it is about 7% higher for the jet-stream cirrus. The important result is the fact that the values of the ratio $P_n(\pi)/P_1(\pi)$, $n = 2, \dots, 5$, are quite small (especially for the water cloud) and they decrease with the order of scattering increasing. If the recommendation of the work of Eloranta (1998) – i.e. “for typical phase functions, $P_n(\pi, h)/P_1(\pi, h)$ is between 0.5 and 1” – is applied, the effect of the multiple scattering will be largely overestimated in the case of the water cloud. In our opinion, the cause of the small values of the ratios $P_n(\pi)/P_1(\pi)$ and of the declining $R_{MSto1}(d)$ in the free atmosphere range is common – i.e. photons, scattered in forward and backward directions, escape the sampling volume.

We performed the same kind of work for the case of space-borne lidar and the configuration described in Sect. 4.1.1. In other words, we provide the values of parameters when the same cloud was probed by the ground-based and the space-borne lidars. The obtained values of parameters are shown in Table A1 as well. It is seen that the values of the ratio $P_n(\pi)/P_1(\pi)$ depend not only on the phase function of particles but also on the lidar configuration – especially on the distance, the RFOV, and the EFOV.

## Appendix B: An example of a simulation uncertainty analysis

It is well known that Monte Carlo modelling provides an estimate of the parameter of interest – i.e. the value that is affected by an error (an estimate minus a true value of the parameter). That error is a random variable whose mean is zero and whose width is characterized by the corresponding variance. The error is due to the statistical nature of Monte Carlo computations (see e.g. Kalos and Whitlock, 2008). In other words, a Monte Carlo estimate, like a measurement, has imperfections that give rise to random error. Such similarity provides grounds to perform a simulation uncertainty analysis using the same approach as in a measurement uncertainty analysis.

Table B1 shows the relative errors (REs) computed according to a document of the Joint Committee for Guides in Metrology (JCGM 100:2008) for our simulation data in the case of space-borne lidar (see Sect. 4.1.1 for details and the notations). The estimate of the variance was obtained using individual simulation replications (see e.g. Bicher et al., 2022). The altitudes of $H = 12.01$ and 6.01 km correspond to the cloud-free molecular atmosphere; the altitudes of $H = 10.99$ and 8.01 km correspond to the near edge and the far edge, respectively, of the cloud. The altitudes of 10.99, 8.01, and 6.01 km correspond to the distances of $d = 0.01, 2.99$, and 5.99 km, respectively, in Fig. 7. The MS contribution to lidar signals and the values of $R_{MSto1}$ are close to
<figure><img src="image_13.png" /></figure>

**Figure A1.** Monte Carlo simulations of multiple-scattering, $R_{MSto1}$ (black points); double-scattering $R_{2to1}$ (red points); and the third-order, $R_{3to1}$ (d) (blue points), and fourth-order, $R_{4to1}$ (d) (green points), respectively, relative contributions to lidar signals. (a, c) Water cloud. (b, d) Jet-stream cirrus. Eloranta model simulations are shown by green (MS), blue (DS), red (third-order), and black (fourth-order) curves.

The results of the Monte Carlo modelling (points) and the EM's simulations (lines) are shown in Fig. A1 only for the in-cloud range of $d \in [0, 3]$ km – i.e. the altitude interval of $H \in [8, 11]$ km. Figure A1a and c correspond to the water cloud and Fig. A1b and d to the cirrus cloud. The microphysical and optical parameters of cloud particles are described in Sect. 2.2 above. The results of each case are divided in two; Fig. A1a and c show the relative contribution of multiple, $R_{MSto1}(d)$ (black points and green line), and double, $R_{2to1}(d)$ (red points and blue line), scattering; Fig. A1b and d show the relative contribution of the third, $R_{3to1}(d)$ (blue points and red line), and the fourth, $R_{4to1}(d)$ (green points and black line), orders. In other words, the lower panels are complementary to the corresponding upper panels. (Note that each panel in Fig. A1 has its own scale of the y axis.)

The good agreement between the MC and EM data is evident. That agreement was obtained by adjusting the EM parameters in the following way. The distance, $h_b$, to the cloud near edge; the extinction coefficient, $\epsilon_p(d)$; the half angle, $\rho_t$, of the receiver field of view; and the half angle, $\rho_l$, of the emitter divergence were assigned according to the configuration used for the MC simulations. First, we found the values of $\Theta_s$ and $P_2(\pi)/P_1(\pi)$ by fitting the MC data on $R_{2to1}(d)$ with the Eloranta model using the ordinary least-squares approach. Then we found the values of $P_3(\pi)/P_1(\pi)$ and $P_4(\pi)/P_1(\pi)$ by fitting the MC data to $R_{3to1}(d)$ and $R_{4to1}(d)$, respectively. We have limited our EM simulations by five orders of scattering. Consequently, the ratio $P_5(\pi)/P_1(\pi)$ was computed to fit the remaining part of the total multiple scattering.

The obtained values of the parameters are shown in Table A1. The value of $\Theta_s$ is in good agreement with Eq. (A2) for the water cloud; it is about 7% higher for the jet-stream cirrus. The important result is the fact that the values of the ratio $P_n(\pi)/P_1(\pi)$, $n = 2, \dots, 5$, are quite small (especially for the water cloud) and they decrease with the order of scattering increasing. If the recommendation of the work of Eloranta (1998) – i.e. “for typical phase functions, $P_n(\pi, h)/P_1(\pi, h)$ is between 0.5 and 1” – is applied, the effect of the multiple scattering will be largely overestimated in the case of the water cloud. In our opinion, the cause of the small values of the ratios $P_n(\pi)/P_1(\pi)$ and of the declining $R_{MSto1}(d)$ in the free atmosphere range is common – i.e. photons, scattered in forward and backward directions, escape the sampling volume.

We performed the same kind of work for the case of space-borne lidar and the configuration described in Sect. 4.1.1. In other words, we provide the values of parameters when the same cloud was probed by the ground-based and the space-borne lidars. The obtained values of parameters are shown in Table A1 as well. It is seen that the values of the ratio $P_n(\pi)/P_1(\pi)$ depend not only on the phase function of particles but also on the lidar configuration – especially on the distance, the RFOV, and the EFOV.

## Appendix B: An example of a simulation uncertainty analysis

It is well known that Monte Carlo modelling provides an estimate of the parameter of interest – i.e. the value that is affected by an error (an estimate minus a true value of the parameter). That error is a random variable whose mean is zero and whose width is characterized by the corresponding variance. The error is due to the statistical nature of Monte Carlo computations (see e.g. Kalos and Whitlock, 2008). In other words, a Monte Carlo estimate, like a measurement, has imperfections that give rise to random error. Such similarity provides grounds to perform a simulation uncertainty analysis using the same approach as in a measurement uncertainty analysis.

Table B1 shows the relative errors (REs) computed according to a document of the Joint Committee for Guides in Metrology (JCGM 100:2008) for our simulation data in the case of space-borne lidar (see Sect. 4.1.1 for details and the notations). The estimate of the variance was obtained using individual simulation replications (see e.g. Bicher et al., 2022). The altitudes of $H = 12.01$ and 6.01 km correspond to the cloud-free molecular atmosphere; the altitudes of $H = 10.99$ and 8.01 km correspond to the near edge and the far edge, respectively, of the cloud. The altitudes of 10.99, 8.01, and 6.01 km correspond to the distances of $d = 0.01, 2.99$, and 5.99 km, respectively, in Fig. 7. The MS contribution to lidar signals and the values of $R_{MSto1}$ are close to
<table><caption>**Table B1.** Relative error of MC simulations in percentages (%).</caption><thead><tr><th></th><th>12.01<br>km</th><th>10.99<br>km</th><th>8.01<br>km</th><th>6.01<br>km</th></tr></thead><tbody><tr><td>Water cloud S<sub>1</sub></td><td>0.37</td><td>0.03</td><td>1.33</td><td>0.77</td></tr><tr><td>Water cloud S<sub>MS</sub></td><td>0.37</td><td>0.03</td><td>6.78</td><td>8.78</td></tr><tr><td>Water cloud R<sub>MSto1</sub></td><td></td><td></td><td>8.68</td><td>15.84</td></tr><tr><td>Jet-stream cirrus S<sub>1</sub></td><td>0.37</td><td>0.05</td><td>2.19</td><td>0.80</td></tr><tr><td>Jet-stream cirrus S<sub>MS</sub></td><td>0.37</td><td>0.05</td><td>10.25</td><td>14.44</td></tr><tr><td>Jet-stream cirrus R<sub>MSto1</sub></td><td></td><td></td><td>11.23</td><td>31.24</td></tr></tbody></table>

zero at altitudes of 12.01 and 10.99 km. Therefore, the corresponding relative errors are not shown in Table B1.

The relative errors, RE, in Table B1 can be discussed on the basis of the following relationship (Rubino and Tuffin, 2009):

$$ RE \sim \frac{1}{\sqrt{N} \cdot \sqrt{\chi}}, \qquad (B1) $$

where N is the sample size – i.e. the number of photons emitted by the lidar ($2 \times 10^{11}$ and $4 \times 10^{11}$ for the water and the cirrus clouds, respectively). The variable $\chi$ is the probability of an event – i.e. the probability that an emitted photon will interact with cloud particles or the molecular atmosphere within the strobe at the altitude of interest and will be scattered towards the lidar within the RFOV. It is evident that MC modelling of lidar signals deals with rare events – i.e. the probability $\chi$ is extremely low. $\chi$ depends on particle characteristics; for example, it is inversely proportional to the lidar ratio. It is proportional to the extinction coefficient. Therefore, the REs are much better at the altitude of 10.99 km compared to 12.01 km. The probability $\chi$ is proportional to the strobe width. For example, it follows from Eq. (B1) that the relative error would be 2 times lower if the strobe were 4 times wider – i.e. the vertical resolution was 4 times worse.

To conclude this section, we underline that the quality of MC modelling can be judged by the dispersion of the data points in a figure and by values of relative errors. That is why we did not perform any smoothing of our MC data, leaving the reader to see the data quality in our figures.

**Data availability.** The data from the simulations are available from the corresponding author upon request.

**Supplement.** The supplement related to this article is available online at: https://doi.org/10.5194/amt-17-3011-2024-supplement.

**Author contributions.** VS, FS, GM, and CC contributed to developing the McRALI code, numerical simulations, and data treatment. VS wrote the paper with help from FS. FS and GM acquired the funding.

**Competing interests.** The contact author has declared that none of the authors has any competing interests.

**Disclaimer.** Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

**Acknowledgements.** This work is part of the French scientific community EECLAT project (Expecting EarthCARE, Learning from A-Train). The EECLAT community and research activities are supported by the National Centre for Space Studies (CNES) and the National Institute for Earth Sciences and Astronomy (INSU).

**Financial support.** This research has been partly supported by the National Institute for Earth Sciences and Astronomy (INSU grant).

**Review statement.** This paper was edited by Dmitry Efremenko and reviewed by two anonymous referees.

## References

Alkasem, A., Szczap, F., Cornet, C., Shcherbakov, V., Gour, Y., Jourdan, O., Labonnote, L. C., and Mioche, G.: Effects of cirrus heterogeneity on lidar CALIOP/CALIPSO data, J. Quant. Spectrosc. Ra., 202, 38–49, https://doi.org/10.1016/j.jqsrt.2017.07.005, 2017.

Bicher, M., Wastian, M., Brunmeir, D., and Popper N.: Review on Monte Carlo simulation stopping rules: how many samples are really enough?, SNE Simul. Notes Eur., 32, 1–8, https://doi.org/10.11128/sne.32.on.10591, 2022.

Bissonnette, L. R.: Lidar and Multiple Scattering, in: Lidar, edited by: Weitkamp, C., Springer Series in Optical Sciences, Springer, New York, NY, Vol. 102, 43–103, https://doi.org/10.1007/0-387-25101-4_3, 2005.

Bissonnette, L. R., Bruscaglioni, P., Ismaelli, A., Zaccanti, G., Cohen, A., Benayahu, Y., Kleiman, M., Egert, S., Flesia, C., Schwendimann, P., Starkov, A. V., Noormohammadian, M., Oppeel, U. G., Winker, D. M., Zege, E. P., Katsev, I. L., and Polonsky, I. N.: LIDAR multiple scattering from clouds, Appl. Phys. B., 60, 355–362, https://doi.org/10.1007/BF01082271, 1995.

Cornet, C., Labonnote, L. C., and Szczap, F.: Three-dimensional polarized Monte Carlo atmospheric radiative transfer model (3DMCPOL): 3D effects on polarized visible reflectances of a cirrus cloud, J. Quant. Spectrosc. Ra., 111, 174–186, https://doi.org/10.1016/j.jqsrt.2009.06.013, 2010.

Delanoë, J. and Hogan, R. J. A variational scheme for retrieving ice cloud properties from combined radar, lidar, and infrared radiometer, J. Geophys. Res., 113, D07204, https://doi.org/10.1029/2007JD009000, 2008.
<table><caption>**Table B1.** Relative error of MC simulations in percentages (%).</caption><thead><tr><th></th><th>12.01<br>km</th><th>10.99<br>km</th><th>8.01<br>km</th><th>6.01<br>km</th></tr></thead><tbody><tr><td>Water cloud S<sub>1</sub></td><td>0.37</td><td>0.03</td><td>1.33</td><td>0.77</td></tr><tr><td>Water cloud S<sub>MS</sub></td><td>0.37</td><td>0.03</td><td>6.78</td><td>8.78</td></tr><tr><td>Water cloud R<sub>MSto1</sub></td><td></td><td></td><td>8.68</td><td>15.84</td></tr><tr><td>Jet-stream cirrus S<sub>1</sub></td><td>0.37</td><td>0.05</td><td>2.19</td><td>0.80</td></tr><tr><td>Jet-stream cirrus S<sub>MS</sub></td><td>0.37</td><td>0.05</td><td>10.25</td><td>14.44</td></tr><tr><td>Jet-stream cirrus R<sub>MSto1</sub></td><td></td><td></td><td>11.23</td><td>31.24</td></tr></tbody></table>

zero at altitudes of 12.01 and 10.99 km. Therefore, the corresponding relative errors are not shown in Table B1.

The relative errors, RE, in Table B1 can be discussed on the basis of the following relationship (Rubino and Tuffin, 2009):

$$ RE \sim \frac{1}{\sqrt{N} \cdot \sqrt{\chi}}, \qquad (B1) $$

where N is the sample size – i.e. the number of photons emitted by the lidar ($2 \times 10^{11}$ and $4 \times 10^{11}$ for the water and the cirrus clouds, respectively). The variable $\chi$ is the probability of an event – i.e. the probability that an emitted photon will interact with cloud particles or the molecular atmosphere within the strobe at the altitude of interest and will be scattered towards the lidar within the RFOV. It is evident that MC modelling of lidar signals deals with rare events – i.e. the probability $\chi$ is extremely low. $\chi$ depends on particle characteristics; for example, it is inversely proportional to the lidar ratio. It is proportional to the extinction coefficient. Therefore, the REs are much better at the altitude of 10.99 km compared to 12.01 km. The probability $\chi$ is proportional to the strobe width. For example, it follows from Eq. (B1) that the relative error would be 2 times lower if the strobe were 4 times wider – i.e. the vertical resolution was 4 times worse.

To conclude this section, we underline that the quality of MC modelling can be judged by the dispersion of the data points in a figure and by values of relative errors. That is why we did not perform any smoothing of our MC data, leaving the reader to see the data quality in our figures.

**Data availability.** The data from the simulations are available from the corresponding author upon request.

**Supplement.** The supplement related to this article is available online at: https://doi.org/10.5194/amt-17-3011-2024-supplement.

**Author contributions.** VS, FS, GM, and CC contributed to developing the McRALI code, numerical simulations, and data treatment. VS wrote the paper with help from FS. FS and GM acquired the funding.

**Competing interests.** The contact author has declared that none of the authors has any competing interests.

**Disclaimer.** Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

**Acknowledgements.** This work is part of the French scientific community EECLAT project (Expecting EarthCARE, Learning from A-Train). The EECLAT community and research activities are supported by the National Centre for Space Studies (CNES) and the National Institute for Earth Sciences and Astronomy (INSU).

**Financial support.** This research has been partly supported by the National Institute for Earth Sciences and Astronomy (INSU grant).

**Review statement.** This paper was edited by Dmitry Efremenko and reviewed by two anonymous referees.

## References

Alkasem, A., Szczap, F., Cornet, C., Shcherbakov, V., Gour, Y., Jourdan, O., Labonnote, L. C., and Mioche, G.: Effects of cirrus heterogeneity on lidar CALIOP/CALIPSO data, J. Quant. Spectrosc. Ra., 202, 38–49, https://doi.org/10.1016/j.jqsrt.2017.07.005, 2017.

Bicher, M., Wastian, M., Brunmeir, D., and Popper N.: Review on Monte Carlo simulation stopping rules: how many samples are really enough?, SNE Simul. Notes Eur., 32, 1–8, https://doi.org/10.11128/sne.32.on.10591, 2022.

Bissonnette, L. R.: Lidar and Multiple Scattering, in: Lidar, edited by: Weitkamp, C., Springer Series in Optical Sciences, Springer, New York, NY, Vol. 102, 43–103, https://doi.org/10.1007/0-387-25101-4_3, 2005.

Bissonnette, L. R., Bruscaglioni, P., Ismaelli, A., Zaccanti, G., Cohen, A., Benayahu, Y., Kleiman, M., Egert, S., Flesia, C., Schwendimann, P., Starkov, A. V., Noormoammadian, M., Oppeel, U. G., Winker, D. M., Zege, E. P., Katsev, I. L., and Polonsky, I. N.: LIDAR multiple scattering from clouds, Appl. Phys. B., 60, 355–362, https://doi.org/10.1007/BF01082271, 1995.

Cornet, C., Labonnote, L. C., and Szczap, F.: Three-dimensional polarized Monte Carlo atmospheric radiative transfer model (3DMCPOL): 3D effects on polarized visible reflectances of a cirrus cloud, J. Quant. Spectrosc. Ra., 111, 174–186, https://doi.org/10.1016/j.jqsrt.2009.06.013, 2010.

Delanoë, J. and Hogan, R. J. A variational scheme for retrieving ice cloud properties from combined radar, lidar, and infrared radiometer, J. Geophys. Res., 113, D07204, https://doi.org/10.1029/2007JD009000, 2008.
Donovan, D. P.: The expected impact of multiple scattering on ATLID signals, EPJ Web of Conferences, 119, 01006, https://doi.org/10.1051/epjconf/201611901006, 2016.

Donovan, D. P. and Van Lammeren, A. C. A. P.: Cloud effective particle size and water content profile retrievals using combined lidar and radar observations: 1. Theory and examples, J. Geophys. Res., 106, 27425–27448, https://doi.org/10.1029/2001JD900243, 2001.

Eloranta, E.: Practical model for the calculation of multiply scattered lidar returns, Appl. Opt., 37, 2464–2472, https://doi.org/10.1364/AO.37.002464, 1998.

Eloranta, E. W.: A Practical Model for the Calculation of Multiply Scattered Lidar Returns-Errata and Extensions to Applied Optics paper, https://lidar.ssec.wisc.edu/multiple_scatter/AO_errata_16.pdf (last access: 22 April 2023), 2000.

Flesia, C. and Starkov, A. V.: Multiple scattering from clear atmosphere obscured by transparent crystal clouds in satellite-borne lidar sensing, Appl. Opt. 35, 2637–2641, https://doi.org/10.1364/AO.35.002637, 1996.

Freville, P., Montoux, N., Baray, J.-L., Chauvigné, A., Réveret, F., Hervo, M., Dionisi, D., Payen, G., and Sellegri, K.: LIDAR Developments at Clermont-Ferrand – France for Atmospheric Observation, Sensors, 15, 3041–3069, https://doi.org/10.3390/s150203041, 2015.

Garnier, A., Pelon, J., Vaughan, M. A., Winker, D. M., Trepte, C. R., and Dubuisson, P.: Lidar multiple scattering factors inferred from CALIPSO lidar and IIR retrievals of semi-transparent cirrus cloud optical depths over oceans, Atmos. Meas. Tech., 8, 2759–2774, https://doi.org/10.5194/amt-8-2759-2015, 2015.

Giannakaki, E., Balis, D. S., Amiridis, V., and Kazadzis, S.: Optical and geometrical characteristics of cirrus clouds over a Southern European lidar station, Atmos. Chem. Phys., 7, 5519–5530, https://doi.org/10.5194/acp-7-5519-2007, 2007.

Hogan, R. J.: Fast approximate calculation of multiply scattered lidar returns, Appl. Opt., 45, 5984–5992, https://doi.org/10.1364/AO.45.005984, 2006.

Hogan, R. J.: Fast lidar and radar multiple-scattering models: Part 1: Small-angle scattering using the photon variance-covariance method, J. Atmos. Sci., 65, 3621–3635, https://doi.org/10.1175/2008JAS2642.1, 2008.

JCGM 100:2008: Evaluation of measurement data – Guide to the expression of uncertainty in measurement, (GUM 1995 with minor corrections), Joint Committee for Guides in Metrology, https://www.bipm.org/documents/20126/2071204/JCGM_100_2008_E.pdf/cb0ef43f-baa5-11cf-3f85-4dcd86f77bd6 (last access: 13 May 2024), 2008.

Kalos, M. H. and Whitlock, P. A.: Monte Carlo Methods, Wiley-VCH Verlag GmbH & Co. KGaA, ISBN 9783527407606, 2008.

Marchuk, G. I., Mikhailov, G. A., Nazareliev, M., Darbinjan, R. A., Kargin, B. A., and Elepov, B. S.: The Monte Carlo methods in atmospheric optics, Springer-Verlag, Berlin, Heidelberg, Germany, https://doi.org/10.1007/978-3-540-35237-2, 2013.

Miller, S. D. and Stephens, G. L.: Multiple scattering effects in the lidar pulse stretching problem, J. Geophys. Res., 104, 22205–22219, https://doi.org/10.1029/1999JD900481, 1999.

Neto, F. D. M. and Neto, A. J. S.: An introduction to inverse problems with applications, Springer, Berlin Heidelberg, ISBN 978-3-642-32557-1, 2013.

Nakoudi, K., Stachlewska, I., and Ritter, C.: An extended lidar-based cirrus cloud retrieval scheme: first application over an Arctic site, Opt. Express, 29, 8553–8580, https://doi.org/10.1364/OE.414770, 2021.

NOAA: U.S. Standard Atmosphere: 1976, National Oceanographic and Atmospheric Administration, Washington, DC, USA, NOAA-S/T 76-1562, ASIN: B000MMZKDS, 1976.

Reverdy, M., Chepfer, H., Donovan, D., Noel, V., Cesana, G., Hoareau, C., Chiriaco, M., and Bastin, S.: An Earth-CARE/ATLID simulator to evaluate cloud description in climate models, J. Geophys. Res.-Atmos., 120, 11090–11113, https://doi.org/10.1002/2015JD023919, 2015.

Rodgers, C. D.: Inverse Methods for Atmospheric Sounding: Theory and Practice, World Scientific Pub. Co. Inc, ISBN 978-9810227401, 2000.

Rubino, G. and Tuffin, B. Introduction to Rare Event Simulation, in: Rare Event Simulation using Monte Carlo Methods, edited by: Rubino, G. and Tuffin, B., John Wiley & Sons, Ltd., ISBN 978-0-470-77269-0, 2009.

Seifert, P., Ansmann, A., Muller, D., Wandinger, U., Althausen, D., Heymsfield, A. J., Massie, S. T., and Schmitt, C.: Cirrus optical properties observed with lidar, radiosonde, and satellite over the tropical Indian Ocean during the aerosol-polluted northeast and clean maritime southwest monsoon, J. Geophys. Res., 112, D17205. https://doi.org/10.1029/2006JD008352, 2007.

Shcherbakov, V., Szczap, F., Alkasem, A., Mioche, G., and Cornet, C.: Empirical model of multiple-scattering effect on single-wavelength lidar data of aerosols and clouds, Atmos. Meas. Tech., 15, 1729–1754, https://doi.org/10.5194/amt-15-1729-2022, 2022.

Szczap, F., Alkasem, A., Mioche, G., Shcherbakov, V., Cornet, C., Delanoë, J., Gour, Y., Jourdan, O., Banson, S., and Bray, E.: McRALI: a Monte Carlo high-spectral-resolution lidar and Doppler radar simulator for three-dimensional cloudy atmosphere remote sensing, Atmos. Meas. Tech., 14, 199–221, https://doi.org/10.5194/amt-14-199-2021, 2021.

Van de Hulst, H. C.: Light scattering by small particles, Dover Publications Inc., New York, ISBN 9780486642284, 1981.

Vaughan, M. A., Powell, K. A., Winker, D. M., Hostetler, C. A., Kuehn, R. E., Hunt, W. H., Getzewich, B. J., Young, S. A., Liu, Z., and McGill, M. J.: Fully automated detection of cloud and aerosol layers in the CALIPSO lidar measurements, J. Atmos. Ocean. Tech., 26, 2034–2050, https://doi.org/10.1175/2009JTECHA1228.1, 2009.

Wang, Z., Zhang, J., and Gao, H.: Impacts of laser beam divergence on lidar multiple scattering polarization returns from water clouds, J. Quant. Spectrosc. Ra., 268, 107618, https://doi.org/10.1016/j.jqsrt.2021.107618, 2021.

Weinman, J. A.: Axially symmetric transfer of light through a cloud of anisotropically scattering particles, Icarus, 9, 67-73, https://doi.org/10.1016/0019-1035(68)90006-7, 1968.

Whiteman, D. N., Evans, K. D., Demoz, B., Starr, D. O’C., Eloranta, E. W., Tobin, D., Feltz, W., Jedlovec, G. J., Gutman, S. I., Schwemmer, G. K., Cadirola, M., Melfi, S. H., and Schmidlin, F. J.: Raman lidar measurements of water vapor and cirrus clouds during the passage of Hurricane Bonnie, J. Geophys. Res., 106, 5211–5225, https://doi.org/10.1029/2000JD900621, 2001.
Donovan, D. P.: The expected impact of multiple scattering on ATLID signals, EPJ Web of Conferences, 119, 01006, https://doi.org/10.1051/epjconf/201611901006, 2016.

Donovan, D. P. and Van Lammeren, A. C. A. P.: Cloud effective particle size and water content profile retrievals using combined lidar and radar observations: 1. Theory and examples, J. Geophys. Res., 106, 27425–27448, https://doi.org/10.1029/2001JD900243, 2001.

Eloranta, E.: Practical model for the calculation of multiply scattered lidar returns, Appl. Opt., 37, 2464–2472, https://doi.org/10.1364/AO.37.002464, 1998.

Eloranta, E. W.: A Practical Model for the Calculation of Multiply Scattered Lidar Returns-Errata and Extensions to Applied Optics paper, https://lidar.ssec.wisc.edu/multiple_scatter/AO_errata_16.pdf (last access: 22 April 2023), 2000.

Flesia, C. and Starkov, A. V.: Multiple scattering from clear atmosphere obscured by transparent crystal clouds in satellite-borne lidar sensing, Appl. Opt. 35, 2637–2641, https://doi.org/10.1364/AO.35.002637, 1996.

Freville, P., Montoux, N., Baray, J.-L., Chauvigné, A., Réveret, F., Hervo, M., Dionisi, D., Payen, G., and Sellegri, K.: LIDAR Developments at Clermont-Ferrand – France for Atmospheric Observation, Sensors, 15, 3041–3069, https://doi.org/10.3390/s150203041, 2015.

Garnier, A., Pelon, J., Vaughan, M. A., Winker, D. M., Trepte, C. R., and Dubuisson, P.: Lidar multiple scattering factors inferred from CALIPSO lidar and IIR retrievals of semi-transparent cirrus cloud optical depths over oceans, Atmos. Meas. Tech., 8, 2759–2774, https://doi.org/10.5194/amt-8-2759-2015, 2015.

Giannakaki, E., Balis, D. S., Amiridis, V., and Kazadzis, S.: Optical and geometrical characteristics of cirrus clouds over a Southern European lidar station, Atmos. Chem. Phys., 7, 5519–5530, https://doi.org/10.5194/acp-7-5519-2007, 2007.

Hogan, R. J.: Fast approximate calculation of multiply scattered lidar returns, Appl. Opt., 45, 5984–5992, https://doi.org/10.1364/AO.45.005984, 2006.

Hogan, R. J.: Fast lidar and radar multiple-scattering models: Part 1: Small-angle scattering using the photon variance-covariance method, J. Atmos. Sci., 65, 3621–3635, https://doi.org/10.1175/2008JAS2642.1, 2008.

JCGM 100:2008: Evaluation of measurement data – Guide to the expression of uncertainty in measurement, (GUM 1995 with minor corrections), Joint Committee for Guides in Metrology, https://www.bipm.org/documents/20126/2071204/JCGM_100_2008_E.pdf/cb0ef43f-baa5-11cf-3f85-4dcd86f77bd6 (last access: 13 May 2024), 2008.

Kalos, M. H. and Whitlock, P. A.: Monte Carlo Methods, Wiley-VCH Verlag GmbH & Co. KGaA, ISBN 9783527407606, 2008.

Marchuk, G. I., Mikhailov, G. A., Nazareliev, M., Darbinjan, R. A., Kargin, B. A., and Elepov, B. S.: The Monte Carlo methods in atmospheric optics, Springer-Verlag, Berlin, Heidelberg, Germany, https://doi.org/10.1007/978-3-540-35237-2, 2013.

Miller, S. D. and Stephens, G. L.: Multiple scattering effects in the lidar pulse stretching problem, J. Geophys. Res., 104, 22205–22219, https://doi.org/10.1029/1999JD900481, 1999.

Neto, F. D. M. and Neto, A. J. S.: An introduction to inverse problems with applications, Springer, Berlin Heidelberg, ISBN 978-3-642-32557-1, 2013.

Nakoudi, K., Stachlewska, I., and Ritter, C.: An extended lidar-based cirrus cloud retrieval scheme: first application over an Arctic site, Opt. Express, 29, 8553–8580, https://doi.org/10.1364/OE.414770, 2021.

NOAA: U.S. Standard Atmosphere: 1976, National Oceanographic and Atmospheric Administration, Washington, DC, USA, NOAA-S/T 76-1562, ASIN: B000MMZKDS, 1976.

Reverdy, M., Chepfer, H., Donovan, D., Noel, V., Cesana, G., Hoareau, C., Chiriaco, M., and Bastin, S.: An Earth-CARE/ATLID simulator to evaluate cloud description in climate models, J. Geophys. Res.-Atmos., 120, 11090–11113, https://doi.org/10.1002/2015JD023919, 2015.

Rodgers, C. D.: Inverse Methods for Atmospheric Sounding: Theory and Practice, World Scientific Pub. Co. Inc, ISBN 978-9810227401, 2000.

Rubino, G. and Tuffin, B. Introduction to Rare Event Simulation, in: Rare Event Simulation using Monte Carlo Methods, edited by: Rubino, G. and Tuffin, B., John Wiley & Sons, Ltd., ISBN 978-0-470-77269-0, 2009.

Seifert, P., Ansmann, A., Muller, D., Wandinger, U., Althausen, D., Heymsfield, A. J., Massie, S. T., and Schmitt, C.: Cirrus optical properties observed with lidar, radiosonde, and satellite over the tropical Indian Ocean during the aerosol-polluted northeast and clean maritime southwest monsoon, J. Geophys. Res., 112, D17205. https://doi.org/10.1029/2006JD008352, 2007.

Shcherbakov, V., Szczap, F., Alkasem, A., Mioche, G., and Cornet, C.: Empirical model of multiple-scattering effect on single-wavelength lidar data of aerosols and clouds, Atmos. Meas. Tech., 15, 1729–1754, https://doi.org/10.5194/amt-15-1729-2022, 2022.

Szczap, F., Alkasem, A., Mioche, G., Shcherbakov, V., Cornet, C., Delanoë, J., Gour, Y., Jourdan, O., Banson, S., and Bray, E.: McRALI: a Monte Carlo high-spectral-resolution lidar and Doppler radar simulator for three-dimensional cloudy atmosphere remote sensing, Atmos. Meas. Tech., 14, 199–221, https://doi.org/10.5194/amt-14-199-2021, 2021.

Van de Hulst, H. C.: Light scattering by small particles, Dover Publications Inc., New York, ISBN 9780486642284, 1981.

Vaughan, M. A., Powell, K. A., Winker, D. M., Hostetler, C. A., Kuehn, R. E., Hunt, W. H., Getzewich, B. J., Young, S. A., Liu, Z., and McGill, M. J.: Fully automated detection of cloud and aerosol layers in the CALIPSO lidar measurements, J. Atmos. Ocean. Tech., 26, 2034–2050, https://doi.org/10.1175/2009JTECHA1228.1, 2009.

Wang, Z., Zhang, J., and Gao, H.: Impacts of laser beam divergence on lidar multiple scattering polarization returns from water clouds, J. Quant. Spectrosc. Ra., 268, 107618, https://doi.org/10.1016/j.jqsrt.2021.107618, 2021.

Weinman, J. A.: Axially symmetric transfer of light through a cloud of anisotropically scattering particles, Icarus, 9, 67-73, https://doi.org/10.1016/0019-1035(68)90006-7, 1968.

Whiteman, D. N., Evans, K. D., Demoz, B., Starr, D. O’C., Eloranta, E. W., Tobin, D., Feltz, W., Jedlovec, G. J., Gutman, S. I., Schwemmer, G. K., Cadirola, M., Melfi, S. H., and Schmidlin, F. J.: Raman lidar measurements of water vapor and cirrus clouds during the passage of Hurricane Bonnie, J. Geophys. Res., 106, 5211–5225, https://doi.org/10.1029/2000JD900621, 2001.
Winker, D. M.: Accounting for multiple scattering in retrievals from space lidar, Proc. SPIE Int. Soc. Opt. Eng., 5059, 128–139, https://doi.org/10.1117/12.512352, 2003.

Winker, D. M., Vaughan, M. A., Omar, A., Hu, Y., Powell, K. A., Liu, Z., Hunt, W. H., and Young, S. A.: Overview of the CALIPSO mission and CALIOP data processing algorithms, J. Atmos. Ocean. Tech., 26, 2310–2323, https://doi.org/10.1175/2009JTECHA1281.1, 2009.

Yang, P. and Liou, K. N.: Geometric-optics-integral-equation method for light scattering by nonspherical ice crystals, Appl. Opt., 35, 6568–84, https://doi.org/10.1364/AO.35.006568, 1996.

Young, S. A. and Vaughan, M. A.: The retrieval of profiles of particulate extinction from Cloud-Aerosol Lidar Infrared Pathfinder Satellite Observations (CALIPSO) data: Algorithm description, J. Atmos. Ocean. Tech., 26, 1105–1119, https://doi.org/10.1175/2008JTECHA1221.1, 2009.

Young, S. A., Vaughan, M. A., Kuehn, R. E., and Winker, D. M.: The Retrieval of Profiles of Particulate Extinction from Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) Data: Uncertainty and Error Sensitivity Analyses, J. Atmos. Ocean. Tech., 30, 395–428, https://doi.org/10.1175/JTECH-D-12-00046.1, 2013.

Zhou, C. and Yang, P.: Backscattering peak of ice cloud particles, Opt. Express, 23, 11995–12003, https://doi.org/10.1364/OE.23.011995, 2015.
Winker, D. M.: Accounting for multiple scattering in retrievals from space lidar, Proc. SPIE Int. Soc. Opt. Eng., 5059, 128–139, https://doi.org/10.1117/12.512352, 2003.

Winker, D. M., Vaughan, M. A., Omar, A., Hu, Y., Powell, K. A., Liu, Z., Hunt, W. H., and Young, S. A.: Overview of the CALIPSO mission and CALIOP data processing algorithms, J. Atmos. Ocean. Tech., 26, 2310–2323, https://doi.org/10.1175/2009JTECHA1281.1, 2009.

Yang, P. and Liou, K. N.: Geometric-optics-integral-equation method for light scattering by nonspherical ice crystals, Appl. Opt., 35, 6568–84, https://doi.org/10.1364/AO.35.006568, 1996.

Young, S. A. and Vaughan, M. A.: The retrieval of profiles of particulate extinction from Cloud-Aerosol Lidar Infrared Pathfinder Satellite Observations (CALIPSO) data: Algorithm description, J. Atmos. Ocean. Tech., 26, 1105–1119, https://doi.org/10.1175/2008JTECHA1221.1, 2009.

Young, S. A., Vaughan, M. A., Kuehn, R. E., and Winker, D. M.: The Retrieval of Profiles of Particulate Extinction from Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) Data: Uncertainty and Error Sensitivity Analyses, J. Atmos. Ocean. Tech., 30, 395–428, https://doi.org/10.1175/JTECH-D-12-00046.1, 2013.

Zhou, C. and Yang, P.: Backscattering peak of ice cloud particles, Opt. Express, 23, 11995–12003, https://doi.org/10.1364/OE.23.011995, 2015.