## Abstract

Abstract: A practical model for determining the time-dependent lidar attenuated backscattering coefficient *β* was developed for application to global lidar data. An analytical expression for the high-order phase function was introduced to reduce computational cost for simulating the angular distribution of the multiple scattering irradiance. The decay rate of the multiple scattering backscattered irradiance was expressed by incorporating the dependence on the scattering angle and the scattering order based on the path integral approach. The estimated *β* over time and the actual range showed good agreement with Monte Carlo simulations for vertically homogeneous and inhomogeneous cloud profiles, resulting in about 15% mean relative error corresponding to 4 times improved accuracy against the Ornstein–Fürth Gaussian approximation method.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

There is a large diversity among the representations of clouds by the major general circulation models (GCMs) [1], with uncertainties in water clouds responsible for the largest variability in future climate predictions [2]. Observations of the vertical profiles of liquid clouds are currently obtained on a global basis by the CloudSat 94-GHz cloud radar and CALIOP lidar onboard the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO), and will be further facilitated by Doppler Cloud Radar and 355-nm high resolution spectral lidar (ATLID) onboard the forthcoming JAXA/ESA EarthCARE mission (2019) [3]. An analysis of water cloud detection by CloudSat and CALIOP has indicated that many water clouds detected by CALIOP are missed by CloudSat; thus, CALIOP-based water cloud retrieval is of great interest. In the Arctic, where super-cooled water and mixed phased clouds have a high frequency of occurrence, quantitative estimates of the microphysical and radiative properties of these clouds by lidar are important to understand their roles in the Arctic climate [4].

For CALIOP, it is known that the attenuated backscattering coefficient *β* and depolarization ratio *δ* from water clouds are strongly affected by multiple scattering events due to the wide footprint size. *δ*–*β* diagrams have proven to be powerful tools for cloud phase discrimination [5,6]. A preliminary case study using Monte Carlo-based multiple scattering look-up-tables has indicated that the combination of CALIOP *β* and *δ* is also effective for retrieving water cloud microphysics (i.e., effective radius *r*_{eff} and liquid water content LWC) [7]. However, in contrast to ice cloud retrievals [8,9], the combination of CALIOP *β* and *δ* values obtained from water clouds have not yet been fully analyzed on a grid-by-grid basis for such purposes at the global scale, due in part to the treatment of multiple scattering.

There have been many attempts to model lidar multiple scattering. For example, a practical model that can simulate the total lidar backscattered intensity in inhomogeneous media was developed by several research groups [10–12]. These were tested against Monte Carlo simulations and have been applied to ground-based lidar inversions. A multiple-scattering solution for a multiple field of view (FOV) lidar system to derive microphysical properties from the total backscattered intensity and secondary-polarized lidar signal was developed by [13] and [14], respectively. A solution for polarized radiative transfer was developed by [15]. The inversion of the vertical profile of the extinction coefficient and cloud macro-properties of optically thick clouds from multiple-scattering lidar has been conducted by a look-up table approach [16] and by diffusion models [17]. For space-borne lidar, a fast forward model for the total attenuated *β* was developed by [18,19] based on the time dependent two stream approximation with the Ornstein–Fürth Gaussian approximation to express the small angle and wide angle scattering (hereafter denoted as the OFGA method).

In this study, the main objective was to construct a practical physical model (PM) for the time-dependent *β* and *δ* from water clouds to capture their transition from the single scattering regime to the multiple scattering regime and to the diffusion regime. It is intended to be applied to global CALIOP and ATLID data to derive vertical profiles of *r*_{eff} and LWC. For this purpose, the PM was evaluated by Monte Carlo simulation. Further, intensive error analyses were conducted to investigate the accuracy of the PM compared to the OFGA method [19] for space-borne lidar configuration. This paper reports on the PM for the *β* part and is organized as follows. In section 2.1 and 2.2, we describe the configuration of the model. In section 2.3, we separate the backscattering coefficient into a contribution from that derivable from the single scattering approximation and a contribution from multiple scattering. The multiple scattering backscattered irradiance was modeled analogously to the single scattering approximation by further introducing the following features. First, instead of the single scattering phase function, a high-order phase function *p*_{n} was introduced to express the normalized angular distribution of the scattered irradiance with respect to the incident direction of the laser beam, according to the number of scattering events. This was based on multiple scattering theory without the small-angle approximation applicable to uniform media [20,21], and was extended to inhomogeneous cloud profiles (section 2.3.2.1). Second, the decay rate of the multiple scattering backscattered irradiance was modeled in an analogous manner to that derived using the path integral approach for a homogenous semi-infinite geometry and with limited boundary conditions for medical applications [22]. The path integral approach was first developed by Feynman [23] to re-formulate Quantum Mechanics and has been widely applied to many fields. Ref. [22] showed that the exponential reduction rate of the diffuse backscattered irradiance was proportional to the first and second power of the normalized photon velocity projected in the direction parallel to that of the incident beam, but was inversely proportional to the number of scattering events; from this, the relationship was modeled (section 2.3.2.2). In section 3, we conducted extensive error analyses of the estimated *β* by the PM for homogeneous and inhomogeneous clouds against Backward Monte Carlo [24] calculations. Results of the error analysis for the PM and the OFGA method over total of 14 different cloud profiles in terms of extinction coefficient *σ*_{ext}, *r*_{eff} and FOV footprint size were characterized for 3 different optical thickness *τ* ranges bounded at *τ* = 2 and *τ* = 8. In addition, the PM was applied to an optically thick molecule layer as an extreme case. Finally, the work is summarized in section 4.

## 2. Physical model (PM)

#### 2.1 Configuration

Before the details of the PM procedure are given, the configuration should be considered. For convenience, a description of the PM is provided here, in which a nadir-pointing space-borne lidar with a laser beam divergence half angle θ_{bd} and FOV half angle θ_{fov} is described as Fig. 1. In a cylindrical polar coordinate system (*r*,*ϕ*,*z*), where *r*,*ϕ* and *z* are the radial distance (*x*^{2} + *y*^{2}), azimuthal coordinate (tan^{−1}(*y*/*x*)) and height, respectively, the satellite is located at (0,0,-*z*_{sat}) and it emits a laser beam in the positive *z* direction. *r*_{bd}(*z*) and *r*_{fov}(*z*) are the radius of the laser beam and FOV footprint size at *z*, respectively. *x* and *y* are orthogonal to the *z*-axis in Cartesian coordinate, and the incident laser beam is linearly polarized in the *x*-direction. The cloud top is located at *z*_{1/2} ( = 0) and the lateral extent of the cloud layer is assumed to be larger than *r*_{fov}(*z*). *z*_{j} (*j* = 1,n) corresponds to the mid-point of the *j*-th cloud layer, and *z*_{j-1/2} and *z*_{j + 1/2} are the bottom and top boundaries of layer *j*, respectively. The geometrical thickness of each layer is determined by the scattering mean free path, i.e., *l*_{j} = 1/*σ*_{sca}(*z*_{j}). *σ*_{sca}(*z*_{j}) is the scattering coefficient at *z*_{j}, and each layer has unit optical thickness from definition. The beam illuminates an area at the cloud top enclosed by the laser beam footprint size with irradiance *I*_{o}. The in-cloud time of the backscattered irradiance is counted from the time a photon of the laser enters the cloud top *t*_{1/2} ( = 0) to the time it exits the cloud. The instant value of the backscattered irradiance *I* is estimated at a discrete time step 2*t*_{j} (*j* = 1, n). *t*_{j} is the half way in-cloud time of the photon trajectory and it is related to *z*_{j} by *z*_{j} = *t*_{j} ⋅c, where c is the speed of light in vacuum. We defined Δ*t*_{j} = *t*_{j} *-t*_{j-1}.

*I* depends on (*r*,*ϕ*,*z,t*) at the maximum z-coordinate position of the photon trajectories. A photon is approximated to reach its maximum *z*-axis position at the half way time of the trajectory, and *I* at 2*t* is denoted by the half way time *t* as *I*(*r*,*ϕ*,*z,t*). In this study, we assumed that *I* is azimuthally independent. In addition, the dependence of *I* on *r* is considered when the dependence of *I* on the FOV is estimated (see section 2.3.2.1). Therefore, *r* and *ϕ* are hereafter omitted from the expression of *I* for simplicity, i.e., *I(z*,*t*).

${\overline{I}}_{i,j}$ represents the backscattered irradiances from layer *z*_{i} at *t*_{j}, which is the sum over *I*(*z,t*) in the range *z*_{i-1/2}<*z*< *z*_{i + 1/2} and *t*_{j-1/2}<*t*< *t*_{j + 1/2} as,

*t*

_{j}is obtained as:

*I*

_{tot}(

*t*

_{j}) is decomposed into the single scattering (ss) component

*I*

_{ss}(

*t*

_{j}) and the multiple scattering (ms) component

*I*

_{ms,tot}(

*t*

_{j}) as follows [12]:

*t*

_{j},

*I*

_{ss}(

*t*

_{j}) only has a contribution from

*z*

_{j}. In contrast

*I*

_{ms,tot}(

*t*

_{j}) has contribution from not only

*z*

_{j}but also from

*z*

_{1}to

*z*

_{j-1}according to the trajectory of the photons contributing to

*I*

_{ms,tot}(

*t*

_{j}). In the PM,

*I*

_{ms,tot}(

*t*

_{j}) is further divided into the on-beam multiple scattering (ms,on) backscattered irradiance

*I*

_{ms,tot,on}(

*t*

_{j}) and the off-beam multiple scattering (ms,off) backscattered irradiance

*I*

_{ms,tot,off}(

*t*

_{j}) according to the contributions as follows:

*I*

_{ms,tot,on}(

*t*

_{j}) is the backscattered irradiance from

*z*

_{j}at

*t*

_{j}, and

*I*

_{ss}scattered at

*z*

_{1}becomes the source for

*I*

_{ms,tot,on}as shown in Fig. 2(a). On the other hand,

*I*

_{ms,tot,off}(

*t*

_{j}) is the backscattered irradiance due to pulse stretching from

*z*

_{1}

*~z*

_{j}at

*t*

_{j}, and basically

*I*

_{ms,tot,on}scattered at previous layers at earlier times become the source for

*I*

_{ms,tot,off}as depicted in Fig. 2(b).

*I*

_{ms,tot,off}(

*t*

_{j}) is the sum of its components

*I*

_{ms,off}(

*z*

_{i}

*,t*

_{j}) from all of the cloud layers as:

*i*=

*j*is included in Eq. (6) to categorize the backscattered irradiances that emerged from

*I*

_{ms,tot,on}into the “off-beam” component.

*I*

_{ms,off}do not become a source for another

*I*

_{ms,off}.

By combining Eqs. (4), (5) and (6),

Example of *I*_{tot} in Eq. (7) at *t*_{3} is shown in Fig. 3. *I*_{tot}(*t*_{3}) consists from five components, *I*_{ss}(*t*_{3}) + *I*_{ms,tot,on}(*t*_{3}) + *I*_{ms,off}(*z*_{1}*,t*_{3}) + *I*_{ms,off}(*z*_{2}*,t*_{3}) + *I*_{ms,off}(*z*_{3}*,t*_{3}), i.e., the single scattering component, the on-beam ms backscattered irradiance form *z*_{3} and those caused by pulse stretching within layers *z*_{1}, *z*_{2} and *z*_{3}. In the PM, the scattering angles of the photon trajectories *Θ*_{on} and *Θ*_{off} were introduced to determine the *z*-*t* dependence of *I*, such that the photon trajectories with larger scattering angles contribute to more off-beam returns as depicted in Fig. 3(b).

The geometrical scattering angles *Θ*_{on} and *Θ*_{off} were used to estimate the source terms for *I*_{ms,tot,on}(*t*_{j}) and *I*_{ms,off}(*z*_{i}*,t*_{j}), respectively.

*Θ*_{on} is defined as follows:

*Θ*

_{on}(

*t*

_{j}) provide the ratio between the total photon path length from the

*z*-axis position of the source

*z*

_{1}to

*z*

_{j}and the distance between

*z*

_{1}and

*z*

_{j}at

*t*

_{j}and

*t*

_{j + 1/2}, respectively. The fraction of

*I*

_{ss}scattered at

*z*

_{1}within the angle 0<

*Θ*<

*Θ*

_{on}becomes the source for

*I*

_{ms,tot,on}(

*t*

_{j}) as in Fig. 4

_{.}

*Θ*_{off} is defined as follows:

${\Theta}_{\text{off,}k\to i}\left({t}_{j-1}\right)=\{\begin{array}{c}0,\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\left(j=i\right)\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\\ {\mathrm{cos}}^{-1}\left[\left({t}_{i}-{t}_{k}\right)c/\left\{\left({t}_{j-1/2}-{t}_{k}\right)c\right\}\right],\text{\hspace{1em}}\left(j>i\right)\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{1em}}\end{array}$(10)

*I*

_{ms,tot,on}(

*t*

_{k}) scattered at

*z*

_{k}within the angle

*Θ*

_{off,k→i}(

*t*

_{j-1})<

*Θ*<

*Θ*

_{off,k→i}(

*t*

_{j}) becomes the source for

*I*

_{ms,off}(

*z*

_{i},

*t*

_{j}). The subscript

*k*→

*i*indicates the scattering from

*z*

_{k}to

*z*

_{i}.

The attenuated backscattering coefficient *β* due to multiple scattering is introduced by an analogous expression to the single scattering lidar equation [25]. In this study, it is assumed that *I*_{o} and *I* are independent of the incident locations within the beam projected area. Therefore, *β* is related to *I* as,

_{s},

*R*

_{s}and A are the temporal pulse length, the range from the source to the target, and the cross-sectional area at

*R*

_{s}, respectively.

*β*

_{w}corresponds to

*I*

_{w}with the same subscript and denotes one of the attenuated backscattering coefficients,

*β*

_{ss}(

*t*

_{j}),

*β*

_{ms,tot}(

*t*

_{j}),

*β*

_{ms,tot,on}(

*t*

_{j}),

*β*

_{ms,tot,off}(

*t*

_{j}),

*β*

_{ms,off}(

*z*

_{i},

*t*

_{j}),

*β*

_{tot}(

*t*

_{j}).

*β*

_{tot}(

*t*

_{j}) is given by,

*σ*

_{ext}(

*z*

_{j}) =

*σ*

_{sca}(

*z*

_{j}). Note that it is straightforward to include absorption [22].

#### 2.2 Introduction of a high-order phase function

The basic concept behind the derivation of *β* lies in the use of the n-th order scattering phase function *p*_{n} rather than tracking the complete scattering processes. *p*_{n}(*Θ*) has been derived by solving the time-dependent scalar transport equation in homogeneous media and provides a normalized angular distribution of the scattered irradiance at angle *Θ* from the initial incident direction *z* after the n-th scattering event. It has been expressed in the Legendre polynomial expansion form by [20,21], and [26] as follows:

*l*+1)A

*/(4π) and*

_{l}*P*are the

_{l}*l*-th order expansion coefficient of the single scattering phase function

*p*

_{ss}(

*z*,

*Θ*) and the Legendre polynomial, respectively.

Figure 5 shows an example of *p*_{n} estimated at different scattering order n in the visible wavelength for effective radius [9] *r*_{eff} of 10µm assuming lognormal size-distribution with the logarithmic standard deviation of 0.33, a typical water cloud particle size. In the visible wavelength, *p*_{n} gradually evolves from the initial forward peaked shape into an isotropic shape as the number of scattering events increases.

#### 2.3 Theoretical procedure for I_{ss} and I_{ms,tot} with path integral approach

_{ss}

_{ms,tot}

The *p*_{n}(*Θ*) is defined relative to the *z*-axis and not to the incident direction from the n-1 th scattering event. Therefore we considered that the backscattered irradiance in the multiple scattering condition could be treated analogously to the single-scattering approximation framework. However, different treatment of the attenuation is needed. That is, here the loss of irradiance with respect to time was estimated by introducing the effective extinction *σ*_{ext,on}(*z*) and *σ*_{ext,off}(*z*) due to multiple scattering for *I*_{ms,tot,on} and *I*_{ms,off}, respectively. This treatment of the effective extinction corresponds to the constant “η factor” approximation to reduce the extinction coefficient introduced in [10]

### 2.3.1 Step 1: *I*_{ss}

*I*_{ss} is obtained by the two-way attenuation along the incident direction in the single scattering approximation for finite layer thickness as in Eq. (10) of [27] as follows:

### 2.3.2 Step 2: *I*_{ms,tot}

The amount of the irradiance scattered out of the initial incident direction at layer *z*_{1} during step 1 is

Part of Δ*I*_{ss} is re-scattered and becomes the source for *I*_{ms,tot} [14]. Because *p*_{n}(*Θ*) determines only the normalized angular distribution of the irradiance, the irradiance must be determined. Here, the procedure to estimate *I*_{ms,tot,on} is first described (see subsection 2.3.2.1.) Then *I*_{ms,tot,off} is estimated from *I*_{ms,tot,on} (see subsection 2.3.2.2).

### 2.3.2.1 Estimation of *I*_{ms,tot,on}

Analogously to *I*_{ss}, *I*_{ms,tot,on} is estimated in a similar form with the single scattering approximation in Eq. (16) as follows:

*B*

_{on}represents the source term of

*I*

_{ms,tot,on}. In Eq. (19), Δ

*I*

_{ss}

*F*

_{on}and

*σ*

_{ext,on}are used instead of

*I*

_{o}and

*σ*

_{ext}in Eq. (16), respectively. Here, we define

*F*

_{on}and

*σ*

_{ext,on}.

- (1)
*F*_{on}: fraction of the source termThe fraction of Δ

*I*_{ss}scattered in the range 0<*Θ*<*Θ*_{on}(*t*_{j}) as in Fig. 4 assures that the scattered irradiance enters layer*z*_{j}at*t*_{j}. Using*p*_{1}and*Θ*_{on},*F*_{on}is estimated by: - (2) σ
_{ext,on}: effective extinctionThe trajectory of the photons contributing to

*I*_{ms,tot,on}does not significantly deviate from those contributing to*I*_{ss}. To estimate*σ*_{ext,on}(*z*_{j}), it is assumed that the extinction coefficient at*z*_{j}is reduced from the single scattering value*σ*_{ext}(*z*_{j}) by the amount of photons scattered at*z*_{j-1}into the effective scattering volume at*z*_{j}[28] as follows:$${\sigma}_{\text{ext,on}}({z}_{j})={\sigma}_{\text{ext}}({z}_{j})\left(1-{\displaystyle {\int}_{0}^{2\pi}{\displaystyle {\int}_{0}^{{\theta}_{d,j}}{p}_{j}\left(\theta \right)}}\mathrm{sin}\theta d\theta d\varphi \right)$$*σ*_{ext}(*z*_{1}) and*p*_{1}in Eq. (21) are treated as the scaled extinction coefficient and geometrically truncated phase function, respectively [29]. The number of scattering events at*t*_{j}equals*j*from the definition of*l*_{j}, and the second term on the right-hand side of Eq. (21) represents the collection of photons scattered within the effective scattering volume defined by the angle*θ*_{d,j}.*θ*_{d,j}changes according to the transition between the geometrical volume determined by geometrical angle*θ*_{ds,j}and the scattering volume determined by the angle of the mean cosine of the*j*-th order phase function 〈*θ*〉_{j}by the tanh function, and is defined as follows:$${\theta}_{\text{d,}\text{\hspace{0.17em}}j}={\theta}_{\text{ds,}\text{\hspace{0.17em}}j}=ta{n}^{-}\left({r}_{\text{fov}}\left({z}_{j}\right)/{l}_{j}\right)\text{\hspace{0.17em}},\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{0.17em}}\text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{0.17em}}\text{\hspace{1em}}\text{\hspace{1em}}\left(j\le 2\right)$$$$\begin{array}{l}{\theta}_{\text{d,}\text{\hspace{0.17em}}j}=\sqrt{{\theta}_{\text{ds,}\text{\hspace{0.17em}}j}{}^{2}+\left({\theta}_{\mathrm{max}}{}^{2}-{\theta}_{\text{ds,}\text{\hspace{0.17em}}j}{}^{2}\right)\mathrm{tanh}\left((j-1)\cdot {\theta}_{\text{ds,}\text{\hspace{0.17em}}3}/{\u3008\theta \u3009}_{3}\right)},\text{\hspace{1em}}\left(3\le j\right)\\ \text{\hspace{1em}}\text{\hspace{1em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\theta}_{\mathrm{max}}=\mathrm{max}\left({\theta}_{\text{ds,}\text{\hspace{0.17em}}j},{\u3008\theta \u3009}_{j}\right)\end{array}$$The two regimes for

*θ*_{d,j}is explained as follows:- ・ When the optical thickness τ is small (≤2),
*θ*_{d,j}is the geometrical angle*θ*_{ds,j}determined by*r*_{fov}(*z*_{j}) and*l*_{j}as in Fig. 6(a) and Eq. (22). In this case,*θ*_{d,j}determines the angle at which the photons are not scattered out of the footprint volume and then lost, but remain inside the footprint volume when passing through*z*_{j}under the single scattering condition. - ・ As τ becomes larger, gradually the amount of photons scattered into the effective scattering volume increases. In this case, the geometrical angle
*θ*_{ds,j}is insufficient to obtain*σ*_{ext,on}, and the scattering volume should be taken into account, i.e., 〈*θ*〉_{j}= cos^{−1}[∫_{4π}*p*_{j}(Θ)cosΘ*d*Ω].*θ*_{d,j}is obtained by Eq. (23) according to the scattering order*j*and the magnitude of 〈*θ*〉_{j}to*θ*_{ds,j}, i.e.,*θ*_{d,j}=*θ*_{ds,j}when*j*is small and*θ*_{d,j}approaches the maximum of either*θ*_{ds,j}or 〈*θ*〉_{j}as*j*increases as in Fig. 6(b).

### 2.3.2.2 Estimation of *I*_{ms,tot,off}

*I*_{ms,off} (*z*_{i}, *t*_{j}) (i*≤*j) is the sum of *B*_{on}(*t*_{k}) that is scattered at *z*_{1}*≤z*_{k}*≤z*_{i}. It is also approximated in analogy to Eq. (16) as follows:

*B*

_{on}= C

_{s}

*I*

_{o}when

*i*= 1. Definition and the physical meaning of

*F*

_{off},

*σ*

_{ext,off}, and

*R*are provided as follows.

- (1)
*F*_{off}: fraction of the source term*F*_{off}(*Θ*_{off,}_{k}_{→}(_{i}*t*_{j})) is the fraction of*I*_{ms,tot,on}that are scattered at (*z*_{k},*t*_{k}) in the direction*Θ*_{off,}_{k}_{→}(_{i}*t*_{j-1}) <*Θ*<*Θ*_{off,}_{k}_{→}(_{i}*t*_{j}) defined in Eqs. (10) and (11) and incident at (*z*_{i},*t*_{j}). It is estimated as in Eq. (20) as follows:$${F}_{\text{off}}\left({\Theta}_{\text{off},k\to i}\left({t}_{j}\right)\right)={C}_{\text{f}}\frac{1}{4\pi}{\displaystyle {\int}_{0}^{2\pi}{\displaystyle {\int}_{{\Theta}_{\text{off,}\text{\hspace{0.17em}}k\to i}\left({t}_{j-1}\right)}^{{\Theta}_{\text{off,}\text{\hspace{0.17em}}k\to i}\left({t}_{j}\right)}{p}_{k}\left(\Theta \right)\mathrm{sin}\Theta d\Theta d\varphi}}$$C

_{f}is an additional term to Eq. (20). C_{f}= 1 for*t*_{j}≠*t*_{i}while C_{f}becomes proportional to the asymmetry factor*g*_{i}for*t*_{j}*= t*_{i}when the FOV is small (*θ*_{d,k}≤ 〈*θ*〉_{k}) and the source term for*I*_{ms,off}decreases. - (2)
*σ*_{ext,off}(*z*_{q},*z*_{k},*z*_{i},*t*_{j}): the decay rate from*z*_{k}to*z*_{i}.Here the procedure to estimate the loss due to scattering from layer

*k*to layer*i*is described.*B*_{on}in Eq. (19) already includes the loss along the path from layer 1 to*k-*1.*σ*_{ext,off}at*z*_{q}(*k*≤*q*≤*i*) is modeled to be proportional to the first and second power of the cosine of the scattering angle (cos*Θ*_{off,}_{k}_{→}) in Eqs. (10) and (11) but inversely proportional to the mean number of scattering events 〈_{i}*n*_{q}〉 on the basis of the path integral approach [22], which is expressed by the 1/*N*_{q}term in the following equation:$${\sigma}_{\text{ext,off}}({z}_{q},{z}_{k},{z}_{i},{t}_{j})={\sigma}_{\text{ext,on}}({z}_{q})/{N}_{q},\text{\hspace{1em}}\left(k\le q\le i\right)$$*N*_{q}is defined as follows:$${N}_{q}=\u3008{n}_{q}\u3009/\left[{C}_{0,\text{\hspace{0.17em}}i}{\mathrm{cos}}^{2}{\Theta}_{\text{off,}\text{\hspace{0.05em}}\text{\hspace{0.17em}}k\to i}\left({t}_{j-1/2}\right)+{C}_{1,\text{\hspace{0.17em}}i}\mathrm{cos}{\Theta}_{\text{off,}\text{\hspace{0.05em}}\text{\hspace{0.17em}}k\to i}\left({t}_{j-1/2}\right)+{C}_{2,\text{\hspace{0.17em}}i}\right]$$The mean

*z*position of the scattered photons at*t*_{q+}_{1/2}is $\sum}_{h=1}^{q}{l}_{h}{g}_{h-1$ as in [30] and 〈*n*_{q}〉 is defined by,Equation (28) is rounded up to an integer. Because [22] considered semi-infinite geometry, in this study the dependence of

*σ*_{ext,off}on*τ*were incorporated into the coefficients C_{o,i}, C_{1,i}, and C_{2,i}in Eq. (27) as:For the same phase function shape, a photon path with a larger scattering angle

*Θ*_{off,}_{k}_{→}has a smaller cos_{i}*Θ*_{off,}_{k}_{→}value and the loss of backscattered irradiance is reduced compared with the more forward scattering paths. This effect is absent at_{i}*τ*_{i}= 1 (C_{o,i}= C_{1,i}= 0) and becomes effective at large τ_{i}. With regard to C_{2,i}, C_{2,i}= 1 at*τ*_{i}= 1, and C_{2,i}~0 and C_{2, i}~g_{ss}at larger τ for the isotropic and forward-peaked phase functions, respectively, and caused*σ*_{ext,off}to be smaller for the isotropic case than the anisotropic scattering case. Δ*d*(*z*_{q},*z*_{k},*z*_{i},*t*_{j}) in Eq. (24) is the excessive effective photon path length from the geometrical layer thickness*l*_{q}at layer*q*and Δ*d*(*z*_{q},*z*_{k},*z*_{i},*t*_{j}) is obtained by,$$\Delta d\left({z}_{q},{z}_{k},{z}_{i},{t}_{j}\right)=\left({t}_{j}-{t}_{i}\right)c\cdot \frac{{g}_{q}}{i-k+1}$$Here, the total excessive path length of the photons from layer

*k*to layer*i*, (*t*_{j}-*t*_{i})c, is divided equally among the number of layers lying between them, and then multiplied by the asymmetry factor of the*q*-th order*g*_{q}. When the source term and*I*_{ms,off}is from the same layer (*k*=*i*), the N_{q}and*g*_{q}terms in Eqs. (27) and (31) are set to be 1, respectively, and*σ*_{ext,off}=*σ*_{ext,on}. - (3)
*R*(*z*_{i}): backscattered ratioSome fraction

*R*(*z*_{i}) of*I*_{ms,tot,on}incident at layer*i*are backscattered within the angle π-*θ*_{d,i}≤*θ*<π at*t*_{j}.*R*(*z*_{i}) are considered to remain within the FOV, i.e.,$$R\left({z}_{i}\right)={\scriptscriptstyle \frac{1}{4\pi}}{\displaystyle {\int}_{0}^{2\pi}{\displaystyle {\int}_{\pi -{\theta}_{d,i}}^{\pi}{p}_{i}\left(\theta \right)}}\mathrm{sin}\theta d\theta d\varphi $$The photons within the FOV that are collected by the receiver are determined by

*σ*_{sca}(*z*_{i})*p*_{i}(π)/(4π) in Eq. (24), where*p*_{i}(π) is the fraction of photons moving in the direction π from the*z*-axis after the*i*-th scattering event.

In summary, *I*_{tot}(*t*_{j}) can be obtained by substituting *I*_{ss}, *I*_{ms,tot,on} and *I*_{ms.off} estimated by Eqs. (16), (18) and (24) into the following equation,

*β*corresponding to Eq. (33) is obtained through Eq. (12) as,

In the PM, the single scattering properties are calculated for the given water cloud *r*_{eff} and LWC profiles by Mie theory, and *β*_{ss}(*t*_{j}), *β*_{ms,tot,on}(*t*_{j}) and *β*_{ms,off}(*z*_{i},*t*_{j}) are then estimated to obtain *β*_{tot}(*t*_{j}).

## 3. Evaluation by a Monte Carlo simulation

The Backward Monte Carlo method [24] (hereafter referred to as MC) used in this study compared well with analytical solutions as well as other Monte Carlo codes tested against the same condition [24], and has been used to analyze the CALIPSO lidar signal [6]. The *β* for the MC was obtained with the same definition as those for the PM. The results obtained by the PM and MC were denoted by the subscripts PM and MC, respectively. A simulation was performed for the CALIPSO lidar specification at a wavelength of 532 nm, i.e., θ_{fov} = 130µrad, θ_{bd} = 100µrad, and the distance from cloud top -*z*_{sat} = 702km. The incident beam was set to nadir for simplicity.

The following three liquid cloud cases were considered. For all three cases, the cloud top layer was located at a height of 3km from the earth’s surface, and cloud vertical extent of 1.2km was considered as in Fig. 7.

- ・ Case 1: a homogeneous case with constant LWC = 0.1g/m
^{3}and*r*_{eff}= 10µm (*σ*_{ext}~15.7km^{−1}) was considered. In addition, simulations were also performed for the same cloud but with a 10-times larger θ_{fov}(10θ_{fov}case) as well as with different LWC profiles, i.e., LWC~0.019g/m^{3}(*σ*_{ext}= 3km^{−1}) and 0.255g/m^{3}(*σ*_{ext}= 40km^{−1}). - ・ Case 2: inhomogeneous cloud particle size profiles with
*r*_{eff}= 20µm, 5µm, and 10 µm from the cloud top bounded at 240m and 480m with LWC = 0.1 g/m^{3}. - ・ Case 3: an inhomogeneous case with two cloud layers with
*r*_{eff}= 10µm at 0–240 m and 480–1200m from the cloud top, and a molecular layer between 240 and 480 m. The optical properties of the molecular layer were calculated for a standard mid-latitude atmospheric profile at 3km.

The validity of the PM were tested for *β*_{tot}(*t*) and *β*(*z*,*t*) against MC simulations with about 3.4 × 10^{8} photons and further compared to the OFGA method.

#### 3.1 Homogeneous profile: Case 1 and sensitivity studies

Here we report on the performance of the PM for Case 1. The *z*-*t* dependence of *β* is shown in Figs. 8(a) and 8(b) with circle symbols and solid lines for the PM and MC, respectively. Colors of the symbols and lines indicate cloud layer numbers, and results for layer 1~7 and 8~17 are shown in Figs. 8(a) and 8(b), respectively. For example, purple color in Fig. 8(a) is *β* from the first layer, and the horizontal axis indicates the apparent *z* position of the return (i.e., time). The first purple circle at *τ* = 1 in Fig. 8(a) corresponds to the sum *β*_{ss}(*t*_{1}) + *β*_{ms,tot,on}(*t*_{1}) + *β*_{ms,off}(*z*_{1},*t*_{1})(=0), while the other purple circles correspond to *β*_{ms,off}(*z*_{1},*t*_{j>1}) in Eq. (34). Figures 8(a) and 8(b) show that *β*_{PM}(*z*,*t*) could capture the features seen in *β*_{MC}(*z*,*t*).

Using the PM, the contributions of ss, ms,on and ms,off to the total *β* were investigated in Fig. 8(c). It was found that the results cloud be categorized into the following three regimes bounded at about τ = 2 and τ = 8: (1) ss dominant, (2) ms,on dominant and (3) ms,off dominant regime, where the contribution from pulse stretching becomes dominant. The slope of *β*(*z*,*t*) simulated by the MC becomes remarkably parabolic for layer number larger than 8 (τ≥8) as seen in Fig. 8(b). This indicated that the contribution of the time-delayed return becomes larger than the on-site returns from about τ≥8, which corresponded to the region where *β*_{ms,tot,off} become to dominate over the other mechanisms in Fig. 8(c).

In Fig. 9(a), *β*_{tot} obtained by the PM for Case 1 as well as for the 10θ_{fov} case were compared to those obtained by the MC and by the OFGA method (*β*_{tot,OFGA}). *β*_{tot,OFGA} was obtained at 1m resolution as the MC. In contrast to the PM where *r*_{eff} and LWC profiles are the input parameters for cloud microphysics, the OFGA method requires the equivalent area radius, *σ*_{ext}, asymmetry factor and the lidar ratio to be specified. The OFGA model simulations were performed for the single scattering properties that were consistent with those used for the PM and MC simulations, expect for the single scattering phase function *p*_{ss}. It is noted that instead of the exact *p*_{ss}, the OFGA method parameterized *p*_{ss} by double Gaussians, and here the option for the representation of the anisotropic phase function in the near-180 degree direction [19] was further selected.

On overall, *β*_{tot,PM} had good agreement with *β*_{tot,MC} for both FOVs as seen in Fig. 9(a). The slope of *β*_{tot,PM} changed according to the three regimes in Fig. 8(c). For Case 1, the regression line for the 1:1 ratio between log_{10}*β*_{tot,PM} and log_{10}*β*_{tot,MC} in Fig. 10, which was forced to pass through the origin, had a slope s = 1.008 and correlation coefficient *r*^{2} = 0.96. In the case of *β*_{tot,OFGA}, deviation from *β*_{tot,MC} occurred from the ms, on dominant regime (2<τ<8) for both FOV cases as seen in Fig. 9(a). In this regime, *β*_{tot,OFGA} underestimated *β*_{tot}.

The relative errors (ERR) of the PM and OFGA corresponding to Fig. 9(a) in respect to MC are compared in Figs. 11(a) and 11(b). ERR is defined as,

_{PM}) and standard deviation (sd) in percentage were 6.39 ± 4.76(%), 6.00 ± 3.00(%), 21.8 ± 14.6(%) at τ≤2, 2<τ<8, 8≤τ, respectively, and 14.4 ± 13.4(%) on average. ERR

_{PM}and sd were further reduced when 10θ

_{fov}was considered, i.e., 8.08 ± 7.07(%). Errors in PM were generally smaller than those of the OFGA method. The (100 × )ERR

_{OFGA}and sd for Case1 were 8.93 ± 6.34%, 31.7 ± 5.12%, 14.7 ± 8.04% and 19.2 ± 11.2% at τ≤2, 2<τ<8, 8≤τ and on average, respectively. Contrary to PM, (100 × )ERR

_{OFGA}and sd slightly increased for the 10θ

_{fov}case, i.e., 23.1 ± 16.8(%). Consequently, the PM especially had about 5 and 7 times better accuracy than the OFGA method for Case1 and 10θ

_{fov}case at the ms,on dominant regime (2<τ<8), respectively, where the transition occurred from single scattering to multiple scattering and to diffusion.

Figure 9(b) shows comparisons among the three methods for *r*_{eff} = 10µm as in Case 1 but with *σ*_{ext} = 3km^{−1} and 40km^{−1}. *β*_{tot,PM} more closely traced *β*_{tot,MC} than *β*_{tot,OFGA} did. *β*_{tot,OFGA} constantly overestimated *β*_{tot} for the 3km^{−1} case. For the σ_{ext} = 40km^{−1}case, *β*_{tot,OFGA} underestimated *β*_{tot} at the 2<τ<8 regime, which was also found in Case 1. (100 × )ERR_{PM} was about 15% on average for both σ_{ext} profiles as shown in Figs. 11(c) and (d). The PM had about 8 times and 4 times better accuracy on average than the OFGA method for the σ_{ext} = 3km^{−1} and 40km^{−1} cases considering τ up to about 18 and 48, respectively.

#### 3.2 Inhomogeneous profile: Case 2

We considered vertically inhomogeneous profiles of cloud microphysics as well as cloud and molecular layers. When the profile is inhomogeneous, the single scattering optical parameters such as *σ*_{ext} and *p*_{ss} change vertically. The *n*-th order phase function for inhomogeneous case is estimated analogously to Eq. (15) as follows:

*A′*is the same as

_{l}*A*but corresponds to the expansion coefficient of

_{l}*p*

_{ss}(

*z*

_{j},

*Θ*), and

*n′*

_{j}is the number of scattering events occurring within layer

*j*(1

*≤j≤n′*). Because the time sampling resolution Δ

*t*

_{j}is determined by the scattering mean free path of the medium, Δ

*t*

_{j}is also not constant. In the PM,

*β*(

*z*) from layer

_{j},t*j*is calculated at the time resolution corresponding to layer

*j*, and then

*β*(

*z*) is interpolated to the finest time step within the profile to estimate

_{j},t*β*

_{tot}(

*t*), contrary to section 3.1.

Figure 12 shows the simulated *β*(*z*,*t*), *β*_{tot} and the relative errors of *β*_{tot} for the PM and the OFGA for inhomogeneous Case 2. In Case 2, *β*_{tot} had a discontinuity in their structures. The PM could reproduce *β*_{tot,MC} and also *β*_{MC}(*z*,*t*) on overall. The second inhomogeneous layer from the cloud top had a larger *σ*_{sca} than the first inhomogeneous layer as shown in Fig. 7(b), and both *β*_{tot,MC} and *β*_{tot,PM} increased at the boundary. *β*_{tot,OFGA} also reproduced the structures seen in *β*_{tot,MC}, but similar tendencies with Case1 were noticed for *β*_{tot,OFGA} to underestimate *β*_{tot,MC} in the ms,on dominant regime, and to overestimate them in other regimes as seen in Figs. 12(c) and 12(e). The mean relative error of PM (12.9%) was smaller than that of the OFGA method (22.3%) on average, and the largest difference in accuracy was found at τ≤2 (4.33% for PM vs. 22.2% for OFGA).

#### 3.3 Inhomogeneous profile: Case 3

*β*_{tot} from the cloud and molecular layers in Case 3 shown in Fig. 13 had more dramatic change in their profiles than Case 2. *β*_{tot,MC} and *β*_{tot,PM} decreased similarly for the molecular layer at 240<z<480m. In this range, ms,off from cloud layer 4 (dark green circle and line in Fig. 13(a)) had major contribution to *β*_{tot}. The mean relative errors in the PM and the OFGA method were 19.0% vs. 69.6%, respectively, and the largest errors in *β*_{tot,OFGA} occurred at the molecule layer between cloud layers as seen in Figs. 13(c) and 13(e). The smaller errors in *β*_{tot,OFGA} in the 600<z<900m region could be explained as follows, i.e., for a homogeneous profile with the same optical properties, underestimation of *β*_{tot} for the OFGA method was found at the optical thickness range corresponding to 600<z<900m as shown in Fig. 9(a) and the effect of the overestimation at 240<z<480m in *β*_{tot,OFGA} could have compensated for the overestimation.

#### 3.4 Rayleigh case with large σ_{ext}

In Cases 1, 2 and 3, the performance of the PM was investigated for *p*_{ss} in the geometrical optics regime, which gradually become isotropic as scattering event increases. In Fig. 14, an optically thick molecular layer was considered as another extreme case, i.e., *p*_{ss} is in the Rayleigh regime, to investigate the applicability of PM when *p*_{ss} was nearly isotropic despite the number of scattering events, but sufficient multiple scattering effects exists [31] on the lidar return. The *σ*_{ext} of the molecule layer was forced to be the same as the cloud in Case 1 in Fig. 8. In this case, *β*(*z*,*t*) for the same color increased with time, which differed from the cases discussed previously. *β*_{tot} in Fig. 14 was larger than that in Case 1, and *β*_{tot,MC} was well reproduced by the PM. *β*_{MC}(*z*,*t*) was also well captured by the PM on overall for layer number larger than 3.

#### 3.5 Characterization of the PM and comparison with the OFGA

The accuracy of the PM was investigated by conducting 14 experiments, and the parameter settings for each experiment were tabulated in Table.1. Experiments 1-12 corresponded to homogeneous cloud profiles with different *σ*_{ext}, *r*_{eff} and FOVs, i.e., *σ*_{ext} = 3km^{−1}, 15.7km^{−1}, 40km^{−1}, *r*_{eff} = 5µm, 10µm, 20µm, and *θ*_{fov} for CALIPSO and 10*θ*_{fov} case. Experiments 13 and 14 corresponded to the inhomogeneous Case2 and Case3, respectively. The overall relative errors 〈ERR_{PM or OFGA}〉 [%] were calculated for different *τ* ranges by averaging |ERR _{PM or OFGA}| estimated at all data points in consideration for the 12 homogeneous cases, and listed in Fig. 15.

〈ERR_{PM}〉 were 15.2% on average and generally around or much smaller than 15% for all optical thickness ranges as shown in Figs. 15(a)-15(d). The accuracy of the PM did not depend on *σ*_{ext} or on *r*_{eff}. On the other hand, 〈ERR_{OFGA}〉 was 65.5% on average and had the largest error in the *σ*_{ext} = 3km^{−1} case. 〈ERR_{OFGA}〉 also had larger standard deviation than 〈ERR_{PM}〉. For the same experiment setting, both 〈ERR_{PM}〉 and 〈ERR_{OFGA}〉 tended to increase with penetration depths, but more noticeable increase was found for the OFGA method. In summary, *β*_{tot,PM} had 4.3 times better accuracy and 11 times smaller standard deviation compared to *β*_{tot,OFGA} on overall, and about 3 times, 4 times and 4.5 times increase in accuracy at τ≤2, τ<8 and 8≤τ, respectively. It is noted that 〈ERR_{PM or OFGA}〉 did not change significantly by including the two inhomogeneous cases, i.e., 15.3% for PM and 63.2% for the OFGA method.

## 4. Summary

A physical model was developed to provide accurate and fast estimates of the time-dependent total multiple-scattered backscattered irradiance and its *z*-*t* dependent components. The following three mechanisms were introduced to capture the major physical processes of lidar multiple scattering in the model.

- (1) The backscatter irradiance was decomposed into three components: ss, ms,on and ms,off.
- (2) The n-th order phase function
*p*_{n}was used. - (3) The path integral formulation was used to estimate the effective extinction.

Evaluation of the PM by the MC simulation was performed for 14 cases: 12 homogenous cases with different *σ*_{ext} profiles, *r*_{eff} and FOVs, case with inhomogeneous cloud particle size profiles, and two-layer cloud case. The comparisons showed overall mean relative error in *β*_{tot} of about 15% for the PM in reference to the MC. Good agreement with the PM and MC in *β*(*z*,*t*) for homogeneous and inhomogeneous cloud profiles was also found. The accuracy of the PM did not strongly depend on the cloud microphysical properties or on the FOV. Further comparison with the Ornstein–Fürth Gaussian approximation (OFGA) method for *β*_{tot} [19] resulted in about 4 times larger error and 11 times larger standard deviation for the OFGA method than the PM. The PM was much more effective than the MC, and could be applied for the analysis of satellites.

The contribution of ss, ms,on and ms,off to *β*_{tot} was investigated by the PM. It was found that *β*_{tot} could be categorized into three regimes according to *τ*. The PM was also tested against an extreme case with particles in the Rayleigh regime (a thick molecule case), and it revealed that the PM had wide applicability.

In future studies, further evaluation of the model outputs and assumptions by MC and speedup of the PM by implementation of more effective integration schemes will be conducted. It is considered that a simulation with higher *z-t* resolution may further improve the accuracy of the proposed method, and for that, a more precise treatment of the probability distribution of the number of scattering event with respect to time [20,21,26] is necessary. Recently, a new type of multiple scattering lidar called the Multiple-FOV Multipe-Scattering Polarization Lidar (MFMSPL) has been developed [32]. This instrument consists of a collection of telescopes that can observe *β* in a comparable manner to those obtained from large satellite footprints with much higher resolution. An evaluation of the estimated *β* by the PM against the MFMSPL observation is also planned for the application of the PM to CALIPSO and ATLID global data.

Further extension of the concept of *p*_{n} to the other scattering phase matrix elements is possible to obtain the parallel and perpendicular components of *β*, and therefore *δ*, for the three regimes, ss, ms,on and ms,off, in a consistent manner to *β* from the PM, and these will be reported in a forthcoming paper [33].

## Funding

JSPS Kakenhi (JP15K17762, JP17H06139); The Japan Aerospace Exploration Agency for EarthCARE Research Announcement; The Arctic Challenge for Sustainability (ArCS); Collaborated Research Program of Research Institute for Applied Mechanics, Kyushu University (Fukuoka, Japan).

## Acknowledgments

The simulations for the OFGA method in Figs. 9,11,12, 13 and 15 were performed by the multiscatter-1.2.11 code developed and distributed by Dr. R. Hogan at University of Reading.

## References and links

**1. **D. E. Waliser, J.-L. F. Li, C. P. Woods, R. T. Austin, J. Bacmeister, J. Chern, A. D. Genio, J. H. Jiang, Z. Kuang, H. Meng, P. Minnis, S. Platnick, W. B. Rossow, G. L. Stephens, S. Sun-Mack, T. Szedung, W. K. Tao, A. M. Tompkins, D. G. Vane, C. Walker, and D. Wu, “Cloud ice: A climate model challenge with signs and expectations of progress,” J. Geophys. Res. **114**, D00A21 (2009).

**2. **S. Bony and J.-L. Dufresne, “Marine boundary layer clouds at the heart of tropical cloud feedback uncertainties in climate models,” Geophys. Res. Lett. **32**(20), L20806 (2005). [CrossRef]

**3. **A. J. Illingworth, H. W. Barker, A. Beljaars, M. Ceccaldi, H. Chepfer, N. Clerbaux, J. Cole, J. Delanoë, C. Domenech, D. P. Donovan, S. Fukuda, M. Hirakata, R. J. Hogan, A. Huenerbein, P. Kollias, T. Kubota, T. Nakajima, T. Y. Nakajima, T. Nishizawa, Y. Ohno, H. Okamoto, R. Oki, K. Sato, M. Satoh, M. W. Shephard, A. Velázquez-Blázquez, U. Wandinger, T. Wehr, and G. van Zadelhoff, “The EarthCARE Satellite: The Next Step Forward in Global Measurements of Clouds, Aerosols, Precipitation, and Radiation,” Bull. Am. Meteorol. Soc. **96**, 1311–1332 (2015).

**4. **J. E. Kay, T. L’Ecuyer, H. Chepfer, N. Loeb, A. Morrison, and G. Cesana, “Recent Advances in Arctic Cloud and Climate Research,” Curr. Clim. Change Rep. **2**(4), 159–169 (2016). [CrossRef]

**5. **Y. Hu, “Depolarization ratio-effective lidar ratio relation: Theoretical basis for space lidar cloud phase discrimination,” Geophys. Res. Lett. **34**(11), L11812 (2007). [CrossRef]

**6. **R. Yoshida, H. Okamoto, Y. Hagihara, and H. Ishimoto, “Global analysis of cloud phase and ice crystal orientation from Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) data using attenuated backscattering and depolarization ratio,” J. Geophys. Res. **115**, D00H32 (2010). [CrossRef]

**7. **K. Sato, H. Okamoto, and H. Ishimoto, “Modeling Lidar Multiple Scattering,” EPJ Web of Conferences **119**, 21005 (2016). [CrossRef]

**8. **H. Okamoto, K. Sato, and Y. Hagihara, “Global analysis of ice microphysics from CloudSat and CALIPSO: Incorporation of specular reflection in lidar signals,” J. Geophys. Res. **115**(D22), D22209 (2010). [CrossRef]

**9. **K. Sato and H. Okamoto, “Refinement of global ice microphysics using spaceborne active sensors,” J. Geophys. Res. **116**(D20), D20202 (2011). [CrossRef]

**10. **C. M. Platt, “Remote sounding of High Clouds. III: Monte Carlo Calculations of Multiple-Scattered Lidar Returns,” J. Atmos. Sci. **38**(1), 156–167 (1981). [CrossRef]

**11. **E. W. Eloranta, “Practical model for the calculation of multiply scattered lidar returns,” Appl. Opt. **37**(12), 2464–2472 (1998). [CrossRef] [PubMed]

**12. **L. R. Bissonnette, “Multiple-scattering lidar equation,” Appl. Opt. **35**(33), 6449–6465 (1996). [CrossRef] [PubMed]

**13. **L. R. Bissonnette, G. Roy, L. Poutier, S. G. Cober, and G. A. Isaac, “Multiple-scattering lidar retrieval method: tests on Monte Carlo simulations and comparisons with in situ measurements,” Appl. Opt. **41**(30), 6307–6324 (2002). [CrossRef] [PubMed]

**14. **G. Roy, L. Bissonnette, C. Bastille, and G. Vallée, “Retrieval of droplet-size density distribution from multiple-field-of-view cross-polarized lidar signals: theory and experimental validation,” Appl. Opt. **38**(24), 5202–5211 (1999). [CrossRef] [PubMed]

**15. **E. P. Zege and L. I. Chaikovskaya, “New approach to the polarized radiative transfer problem,” J. Quant. Spectrosc. Radiat. Transf. **55**(1), 19–31 (1996).

**16. **R. F. Cahalan, M. McGill, and J. Kolasinski, “THOR-cloud thickness from offbeam lidar returns,” J. Atmos. Ocean. Technol. **22**(6), 605–627 (2005). [CrossRef]

**17. **A. B. Davis, “Multiple-scattering lidar from both sides of the clouds: Addressing internal structure,” J. Geophys. Res. **113**(D14), D14S10 (2008). [CrossRef]

**18. **R. J. Hogan, “Fast Lidar and Radar Multiple-Scattering Models. Part I: Small-Angle Scattering Using the Photon Variance–Covariance Method,” J. Atmos. Sci. **65**(12), 3621–3635 (2008). [CrossRef]

**19. **R. J. Hogan and A. Battaglia, “Fast Lidar and Radar Multiple-Scattering Models. Part II: Wide-Angle Scattering Using the Time-Dependent Two-Stream Approximation,” J. Atmos. Sci. **65**(12), 3636–3651 (2008). [CrossRef]

**20. **S. Goudsmit and J. L. Saunderson, “Multiple Scattering of Electrons,” Phys. Rev. **57**(1), 24–29 (1940). [CrossRef]

**21. **H. W. Lewis, “Multiple scattering in an infinite medium,” Phys. Rev. **78**(5), 526–529 (1950). [CrossRef]

**22. **L. T. Perelman, J. Wu, I. Itzkan, and M. S. Feld, “Photon Migration in Turbid Media Using Path Integrals,” Phys. Rev. Lett. **72**(9), 1341–1344 (1994). [CrossRef] [PubMed]

**23. **R. P. Feynman, “Space-Time Approach to Quantum Electrodynamics,” Phys. Rev. **76**(6), 769 (1949). [CrossRef]

**24. **H. Ishimoto and K. Masuda, “A Monte Carlo approach for the calculation of polarized light: application to an incident narrow beam,” J. Quant. Spectrosc. Radiat. Transf. **72**(4), 467–483 (2002). [CrossRef]

**25. **K. N. Liou, “*An Introduction to Atmospheric Radiation*,” (Academic Press, 2002).

**26. **J. M. Fernández-Varea, R. Mayol, J. Baró, and F. Salvat, “On the theory and simulation of multiple elastic scattering of electrons,” Nucl. Instrum. Methods Phys. Res. B **73**(4), 447–473 (1993). [CrossRef]

**27. **H. Okamoto, S. Iwasaki, M. Yasui, H. Horie, H. Kuroiwa, and H. Kumagai, “An algorithm for retrieval of cloud microphysics using 95-GHz cloud radar and lidar,” J. Geophys. Res. **108**(4226), D7 (2003).

**28. **C. F. Bohren and D. R. Huffman, *Absorption and Scattering of Light by Small Particles* (Wiley, 1983).

**29. **H. Iwabuchi and T. Suzuki, “Fast and accurate radiance calculations using truncation approximation for anisotropic scattering phase functions,” J. Quant. Spectrosc. Radiat. Transf. **110**(17), 1926–1939 (2009). [CrossRef]

**30. **F. O. Nicolas, L. R. Bissonnette, and P. H. Flamant, “Lidar effective multiple-scattering coefficients in cirrus clouds,” Appl. Opt. **36**(15), 3458–3468 (1997). [CrossRef] [PubMed]

**31. **J. E. Hansen, “Absorption-line formation in a scattering planetary atmosphere: A test of Van de Hulst’s similarity relations,” Astrophys. J. **158**, 337–349 (1969). [CrossRef]

**32. **H. Okamoto, K. Sato, T. Nishizawa, N. Sugimoto, T. Makino, Y. Jin, A. Shimizu, T. Takano, and M. Fujikawa, “Development of a multiple-field-of-view multiple-scattering polarization lidar: comparison with cloud radar,” Opt. Express **24**(26), 30053–30067 (2016). [CrossRef] [PubMed]

**33. **K. Sato, Research Institute for Applied Mechanics, Kyushu University, Kasuga, Fukuoka 816–8580, H. Okamoto, and H. Ishimoto are preparing a manuscript to be called “Modeling the depolarization of space-borne lidar signals.”