Physical model for multiple scattered space-borne lidar returns from clouds

Kaori Sato; Hajime Okamoto; Hiroshi Ishimoto

doi:10.1364/OE.26.00A301

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² + y²), azimuthal coordinate (tan⁻¹(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 2t_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.

Fig. 1 Schematic of the cloud area, the satellite position and the foot print sizes of the beam and field-of-view at cloud position.

Download Full Size | PDF

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 2t 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).

${\bar{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,

{\bar{I}}_{i, j} = \frac{1}{t_{j + 1 / 2} - t_{j - 1 / 2}} \frac{1}{z_{i + 1 / 2} - z_{i - 1 / 2}} \int_{t_{j - 1 / 2}}^{t_{j + 1 / 2}} \int_{z_{i - 1 / 2}}^{z_{i + 1 / 2}} I (z, t) d z d t

In the physical model,

{\bar{I}}_{i, j}

is approximated by,

{\bar{I}}_{i, j} \approx I (z_{i}, t_{j})

The total backscattered irradiance from all of the cloud layers at 2t_j is obtained as:

I_{tot} (t_{j}) = \sum_{i = 1}^{j} I (z_{i}, t_{j})

In the PM, 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]:

I_{tot} (t_{j}) = I_{ss} (t_{j}) + I_{ms,tot} (t_{j})

At 2t_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₁ 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} (t_{j}) = I_{ms, tot, on} (t_{j}) + I_{ms, tot, off} (t_{j})

I_ms,tot,on(t_j) is the backscattered irradiance from z_j at t_j, and I_ss scattered at z₁ 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₁~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_{ms, tot, off} (t_{j}) = \sum_{i = 1}^{j} I_{ms, off} (z_{i}, t_{j})

Here, 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.

Fig. 2 Relation between I_ss (black), I_ms,tot,on (blue) and I_ms,off (green). The direction of the arrows indicates the sources and sinks. Since all of the backscattered irradiances are defined by the maximum z-axis value of the photon trajectory, the arrows are all looking downstream. The z-t dependences of I are indicated by the numbers in parentheses.

Download Full Size | PDF

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

I_{tot} (t_{j}) = I_{ss} (t_{j}) + I_{ms, tot, on} (t_{j}) + \sum_{i = 1}^{j} I_{ms, off} (z_{i}, t_{j})

From Eqs. (3) and (7),

I (z_{i}, t_{j}) = {\begin{matrix} I_{ss} (t_{j}) + I_{ms,tot,on} (t_{j}) + I_{ms,off} (z_{j}, t_{j}), (i = j) \\ I_{ms,off} (z_{i}, t_{j}), (i \neq j) \end{matrix}

Example of I_tot in Eq. (7) at t₃ is shown in Fig. 3. I_tot(t₃) consists from five components, I_ss(t₃) + I_ms,tot,on(t₃) + I_ms,off(z₁,t₃) + I_ms,off(z₂,t₃) + I_ms,off(z₃,t₃), i.e., the single scattering component, the on-beam ms backscattered irradiance form z₃ and those caused by pulse stretching within layers z₁, z₂ and z₃. 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).

Fig. 3 (a.) Schematic of the components of I_tot(t₃):① I_ss, ② I_ms,tot,on and ③,④,⑤ I_ms,off from different cloud layers, but observed at the same time 2t₃. (b.) The components ①-⑤ shown in (a) are modeled by introducing the geometrical scattering angles Θ_on and Θ_off that represent the photon trajectories. The returning paths of I are abbreviated by simple arrows.

Download Full Size | PDF

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}) = (\cos^{- 1} [(t_{j} - t_{1}) c / {(t_{j + 1 / 2} - t_{1}) c}])

The cosine of the minimum and maximum values of Θ_on(t_j) provide the ratio between the total photon path length from the z-axis position of the source z₁ to z_j and the distance between z₁ and z_j at t_j and t_{j + 1/2}, respectively. The fraction of I_ss scattered at z₁ within the angle 0<Θ<Θ_on becomes the source for I_ms,tot,on(t_j) as in Fig. 4_.

Fig. 4 Example of Θ_on and Θ_off in Eqs. (9,10 and 11). Θ_on and Θ_off are determined by the geometrical lengths,①, ②and,③,④,⑤, respectively. The fraction of I_ss(t₁) scattered at z₁ in the direction 0<Θ<Θ_on(t₃) becomes the source for I_ms,tot,on(t₃). I_ms,off (z₄, t₆) and I_ms,off (z₄, t₇) are the time delayed returns from z₄ at different times t₆ and t₇, respectively. The fraction of I_ms,tot,on(t₃) scattered at z₃ within the angle Θ_{off,3 →4}(t₅)<Θ<Θ_{off,3 →4}(t₆) and Θ_{off,3 →4}(t₆)<Θ<Θ_{off,3 →4}(t₇) becomes the source for I_ms,off (z₄, t₆) and I_ms,off (z₄, t₇), respectively. The returning paths of I are abbreviated by simple arrows.

Download Full Size | PDF

Θ_off is defined as follows:

$Θ_{off, k \to i} (t_{j - 1}) = {\begin{matrix} 0, (j = i) \\ \cos^{- 1} [(t_{i} - t_{k}) c / {(t_{j - 1 / 2} - t_{k}) c}], (j > i) \end{matrix}$ (10)

Θ_{off, k \to i} (t_{j}) = \cos^{- 1} [(t_{i} - t_{k}) c / {(t_{j + 1 / 2} - t_{k}) c}]

The fraction of 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,

β_{w} = I_{w} / (I_{0} C_{s})

where

C_{s} = A c τ_{s} / (2 R_{s}^{2})

τ_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,

β_{tot} (t_{j}) = β_{ss} (t_{j}) + β_{ms, tot} (t_{j})

In this study, we omitted absorption along the path for simplicity and therefore σ_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:

p_{n} (Θ) = \sum_{l = 0}^{\infty} \frac{2 l + 1}{4 π} {(A_{l})}^{n} P_{l} (\cos Θ)

where (2l+1)A_l/(4π) and P_l are the 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.

Fig. 5 Normalized n-th order scattering phase function at 532nm for a liquid particle with an effective radius of 10µm. The colors of the lines indicate the scattering order n = 1,3,5,7,10,20,30; a transition in the scattering phase function from a forward-peaked shape to a more isotropic shape was observed. Enlarged figure of the scattering phase functions for n = 20 and 30 is shown in the upper right corner.

Download Full Size | PDF

2.3 Theoretical procedure for I_ss and I_ms,tot with path integral approach

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:

I_{ss} (t_{j}) = C_{s} I_{o} \exp [- 2 τ (z_{j - 1 / 2})] \cdot [\frac{\exp {- 2 σ_{ext} (Z_{j}) l_{j}} - 1}{- 2 σ_{ext} (Z_{j}) l_{j}}] \cdot σ_{sca} (z_{j}) p_{ss} (z_{j}, π) / (4 π)

2.3.2 Step 2: I_ms,tot

The amount of the irradiance scattered out of the initial incident direction at layer z₁ during step 1 is

Δ I_{s s} = I_{o} [1 - {\exp (- 2 σ_{ext} (z_{1}) l_{1}) - 1} / {[- 2 σ_{ext} (z_{1}) l_{1}]}]

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:

I_{ms,tot,on} (t_{j}) = B_{on} (t_{j}) \cdot [\frac{\exp {- 2 σ_{ext,on} (Z_{j}) l_{j}} - 1}{- 2 σ_{ext,on} (Z_{j}) l_{j}}] σ_{sca} (z_{j}) p_{ss} (z_{j}, π) / (4 π)

B_{on} (t_{j}) = C_{s} Δ I_{ss} F_{on} (Θ_{on} (t_{j})) \exp [- 2 \sum_{i = 1}^{j - 1} σ_{ext,on} (z_{i}) l_{i}]

B_on represents the source term of I_ms,tot,on. In Eq. (19), ΔI_ssF_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 term
The 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₁ and Θ_on, F_on is estimated by:
$F_{on} (Θ_{on} (t_{j})) = \frac{1}{4 π} \int_{0}^{2 π} \int_{0}^{Θ_{on} (t_{j})} p_{1} (z_{1}, Θ) \sin Θ d Θ d ϕ$
(2) σ_ext,on: effective extinction
The 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:
(21) $σ_{ext,on} (z_{j}) = σ_{ext} (z_{j}) (1 - \int_{0}^{2 π} \int_{0}^{θ_{d, j}} p_{j} (θ) \sin θ d θ d ϕ)$

σ_ext(z₁) and p₁ 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:
(22) $θ_{d, j} = θ_{ds, j} = t a n^{-} (r_{fov} (z_{j}) / l_{j}), (j \leq 2)$

(23) $\begin{array}{l} θ_{d, j} = \sqrt{θ_{ds, j}^{2} + (θ_{\max}^{2} - θ_{ds, j}^{2}) \tanh ((j - 1) \cdot θ_{ds, 3} / {〈 θ 〉}_{3})}, (3 \leq j) \\ θ_{\max} = \max (θ_{ds, j}, {〈 θ 〉}_{j}) \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⁻¹ [∫_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).

Fig. 6 (a) Geometry for defining θ_ds,j in Eq. (22). (b) Schematic of the dependence of θ_d,j on scattering order j when θ_max = 〈θ〉_j in Eq. (23). For satellites, θ_ds,j is approximately constant for homogenous profiles, but varies for inhomogeneous profiles.

Download Full Size | PDF

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₁≤z_k≤z_i. It is also approximated in analogy to Eq. (16) as follows:

\begin{array}{l} I_{ms,off} (z_{i}, t_{j}) = \sum_{_{k = 1}}^{i} [B_{on} (t_{k}) F_{off} (Θ_{off, k \to i} (t_{j})) R (z_{i}) \\ \times \exp [- 2 \sum_{q = k}^{i} σ_{ext,off} (z_{q}, z_{k}, z_{i}, t_{j}) \cdot (l_{q} + Δ d (z_{q}, z_{k}, z_{i}, t_{j}))] σ_{sca} (z_{i}) p_{i} (π) / (4 π) \end{array}

B_on = C_sI_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_{off} (Θ_{off, k \to i} (t_{j})) = C_{f} \frac{1}{4 π} \int_{0}^{2 π} \int_{Θ_{off, k \to i} (t_{j - 1})}^{Θ_{off, k \to i} (t_{j})} p_{k} (Θ) \sin Θ d Θ d ϕ$
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_→_i) in Eqs. (10) and (11) but inversely proportional to the mean number of scattering events 〈n_q〉 on the basis of the path integral approach [22], which is expressed by the 1/N_q term in the following equation:
$σ_{ext,off} (z_{q}, z_{k}, z_{i}, t_{j}) = σ_{ext,on} (z_{q}) / N_{q}, (k \leq q \leq i)$
N_q is defined as follows:
$N_{q} = 〈 n_{q} 〉 / [C_{0, i} \cos^{2} Θ_{off, k \to i} (t_{j - 1 / 2}) + C_{1, i} \cos Θ_{off, k \to i} (t_{j - 1 / 2}) + C_{2, i}]$
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,
$〈 n_{q} 〉 = \sum_{h = 1}^{q} g_{h - 1}$
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:
$C_{0, i} = C_{1, i} = 1 - {1 / τ}_{i}$ $C_{2, i} = 1 + 1 / τ_{i} - 1 / \sum_{h = 1}^{i} g_{h - 1}$
For the same phase function shape, a photon path with a larger scattering angle Θ_off, _k_→_i has a smaller cosΘ_off, _k_→_i value and the loss of backscattered irradiance is reduced compared with the more forward scattering paths. This effect is absent at τ_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,
$Δ d (z_{q}, z_{k}, z_{i}, t_{j}) = (t_{j} - t_{i}) 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 ratio
Some 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 (z_{i}) = \frac{1}{4 π} \int_{0}^{2 π} \int_{π - θ_{d, i}}^{π} p_{i} (θ) \sin θ d θ d ϕ$
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,

\begin{array}{l} I_{tot} (t_{j}) = I_{ss} (t_{j}) + I_{ms,tot} (t_{j}) \\ = I_{ss} (t_{j}) + I_{ms, tot,on} (t_{j}) + \sum_{i = 1}^{j} I_{ms, off} (z_{i}, t_{j}) \end{array}

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

β_{tot} (t_{j}) = β_{ss} (t_{j}) + β_{ms, tot,on} (t_{j}) + \sum_{i = 1}^{j} β_{ms, off} (z_{i}, t_{j})

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.

Fig. 7 Vertical profile of (a.) cloud r_eff and (b.) σ_ext for Cases 1, 2, and 3.

Download Full Size | PDF

・ Case 1: a homogeneous case with constant LWC = 0.1g/m³ and r_eff = 10µm (σ_ext~15.7km⁻¹) 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³ (σ_ext = 3km⁻¹) and 0.255g/m³ (σ_ext = 40km⁻¹).
・ 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³.
・ 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⁸ 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₁) + β_ms,tot,on(t₁) + β_ms,off(z₁,t₁)(=0), while the other purple circles correspond to β_ms,off(z₁,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).

Fig. 8 (a, b) Comparison between the physical model (PM) and the Monte Carlo (MC) simulation for Case 1 for β(z,t). In (a) and (b), the colors indicate β(z,t) from cloud layers 1~7 and 8~17, respectively. (c.) The contributions of ss, ss + ms,on and ms,off to β_tot. (d) Same as (a,b) but with reduced number of lines.

Download Full Size | PDF

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.

Fig. 9 Comparison between β_tot obtained by the physical model (PM: red solid and dotted lines,) and the Monte Carlo simulation (MC: solid black and gray lines) for (a) Case 1 as well as simulation with a 10-times larger field of view angle (10θ_fov) and for (b) water cloud cases with the same r_eff = 10 µm as (a) but with extinction coefficient of 3km⁻¹ and 40km⁻¹ at the 532nm CALIPSO lidar specification. In (a) and (b), β_tot simulated by the OFGA method [19] (blue solid and dotted lines) are also shown.

Download Full Size | PDF

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₁₀β_tot,PM and log₁₀β_tot,MC in Fig. 10, which was forced to pass through the origin, had a slope s = 1.008 and correlation coefficient r² = 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.

Fig. 10 Comparison between log₁₀β_tot,MC and log₁₀β_tot,PM in Fig. 9(a) for Case1. Black and red lines correspond to a 1:1 ratio and the slope (s) of the linear regression line forced through the origin with correlation coefficient r², respectively.

Download Full Size | PDF

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,

{ERR}_{PM or OFGA} = (β_{PM or OFGA} - β_{MC}) / β_{MC}

The mean relative error (100 × ERR_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.

Fig. 11 Relative error of β_tot,PM and β_tot,OFGA corresponding to (a,b) Fig. 9(a) and (c,d) Fig. 9(b) as functions of the optical thickness τ. In (c), τ range larger than that corresponding to Fig. 9(b) as indicated by the black dotted line is shown. The mean relative errors <ERR> and standard deviations (sd) in percentage for each case at 4 different τ ranges are shown in each figure.

Download Full Size | PDF

Figure 9(b) shows comparisons among the three methods for r_eff = 10µm as in Case 1 but with σ_ext = 3km⁻¹ and 40km⁻¹. β_tot,PM more closely traced β_tot,MC than β_tot,OFGA did. β_tot,OFGA constantly overestimated β_tot for the 3km⁻¹ case. For the σ_ext = 40km⁻¹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⁻¹ and 40km⁻¹ 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:

p_{n} (Θ) = \sum_{l = 0}^{\infty} \frac{2 l + 1}{4 π} {\prod_{j = 1}^{n^{'}} {({A^{'}}_{l} (z_{j}))}^{{n^{'}}_{j}}} P_{l} (\cos Θ)

A′_l is the same as A_l but corresponds to the expansion coefficient of 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_j,t) from layer j is calculated at the time resolution corresponding to layer j, and then β(z_j,t) is interpolated to the finest time step within the profile to estimate β_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).

Fig. 12 (a,b,d) Same as Fig. 8(a,b,d) but for Case 2. (c) Same as Fig. 9 but for Case2. (e) Same as Fig. 11 but the relative errors of β_tot corresponding to (c) are shown.

Download Full Size | PDF

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.

Fig. 13 (a,b,c,d,e) Same as Fig. 12(a,b,c,d,e) but for Case 3.

Download Full Size | PDF

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.

Fig. 14 Same as Fig. 13(a), but for a thick molecular layer.

Download Full Size | PDF

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⁻¹, 15.7km⁻¹, 40km⁻¹, 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.

Table 1. Experiment setting

View Table

Fig. 15 Relative error of β_tot againtst Monte Calro simulation for (a-d) the PM and for (e-h) the OFGA method for 14 experiments summarized in Table 1 and categorized according to the optical thickness, i.e., (a,e), 0≤τ (b,f) τ≤2, (c,g) 2<τ<8, (d,h) 8≤τ. Relative errors out of the vertical axis range were indicated in parentheses. Different shaded areas denote simulations performed with different σ_ext profiles. The maximum τ considered are about 48 for the σ_ext = 40km⁻¹ cases and about 18 for the others. The overall mean relative errors 〈ERR〉 and standard deviations [%] listed in the upper right corner of each figure are estimated from all 12 homogeneous profiles.

Download Full Size | PDF

〈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⁻¹ 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.”

Physical model for multiple scattered space-borne lidar returns from clouds

Abstract

1. Introduction