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
There is a large diversity among the representations of clouds by the major general circulation models (GCMs) , with uncertainties in water clouds responsible for the largest variability in future climate predictions . 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) . 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 .
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 reff and liquid water content LWC) . 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  and , respectively. A solution for polarized radiative transfer was developed by . 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  and by diffusion models . 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 reff 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  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 pn 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 126.96.36.199). 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 . The path integral approach was first developed by Feynman  to re-formulate Quantum Mechanics and has been widely applied to many fields. Ref.  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 188.8.131.52). In section 3, we conducted extensive error analyses of the estimated β by the PM for homogeneous and inhomogeneous clouds against Backward Monte Carlo  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, reff 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)
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 (x2 + y2), azimuthal coordinate (tan−1(y/x)) and height, respectively, the satellite is located at (0,0,-zsat) and it emits a laser beam in the positive z direction. rbd(z) and rfov(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 z1/2 ( = 0) and the lateral extent of the cloud layer is assumed to be larger than rfov(z). zj (j = 1,n) corresponds to the mid-point of the j-th cloud layer, and zj-1/2 and zj + 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., lj = 1/σsca(zj). σsca(zj) is the scattering coefficient at zj, 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 Io. The in-cloud time of the backscattered irradiance is counted from the time a photon of the laser enters the cloud top t1/2 ( = 0) to the time it exits the cloud. The instant value of the backscattered irradiance I is estimated at a discrete time step 2tj (j = 1, n). tj is the half way in-cloud time of the photon trajectory and it is related to zj by zj = tj ⋅c, where c is the speed of light in vacuum. We defined Δtj = tj -tj-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 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 184.108.40.206). Therefore, r and ϕ are hereafter omitted from the expression of I for simplicity, i.e., I(z,t).
represents the backscattered irradiances from layer zi at tj, which is the sum over I(z,t) in the range zi-1/2<z< zi + 1/2 and tj-1/2<t< tj + 1/2 as,12]:Fig. 2(a). On the other hand, Ims,tot,off(tj) is the backscattered irradiance due to pulse stretching from z1~zj at tj, and basically Ims,tot,on scattered at previous layers at earlier times become the source for Ims,tot,off as depicted in Fig. 2(b). Ims,tot,off(tj) is the sum of its components Ims,off (zi,tj) from all of the cloud layers as:Eq. (6) to categorize the backscattered irradiances that emerged from Ims,tot,on into the “off-beam” component. Ims,off do not become a source for another Ims,off.
By combining Eqs. (4), (5) and (6),Eqs. (3) and (7),
Example of Itot in Eq. (7) at t3 is shown in Fig. 3. Itot(t3) consists from five components, Iss(t3) + Ims,tot,on(t3) + Ims,off(z1,t3) + Ims,off(z2,t3) + Ims,off(z3,t3), i.e., the single scattering component, the on-beam ms backscattered irradiance form z3 and those caused by pulse stretching within layers z1, z2 and z3. 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 Ims,tot,on(tj) and Ims,off(zi,tj), respectively.
Θon is defined as follows:Fig. 4.
Θoff is defined as follows:
The attenuated backscattering coefficient β due to multiple scattering is introduced by an analogous expression to the single scattering lidar equation . In this study, it is assumed that Io and I are independent of the incident locations within the beam projected area. Therefore, β is related to I as,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 pn rather than tracking the complete scattering processes. pn(Θ) 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  as follows:
Figure 5 shows an example of pn estimated at different scattering order n in the visible wavelength for effective radius  reff 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, pn gradually evolves from the initial forward peaked shape into an isotropic shape as the number of scattering events increases.
2.3 Theoretical procedure for Iss and Ims,tot with path integral approach
The pn(Θ) 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 Ims,tot,on and Ims,off, respectively. This treatment of the effective extinction corresponds to the constant “η factor” approximation to reduce the extinction coefficient introduced in 
2.3.1 Step 1: Iss
Iss 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  as follows:
2.3.2 Step 2: Ims,tot
The amount of the irradiance scattered out of the initial incident direction at layer z1 during step 1 is
Part of ΔIss is re-scattered and becomes the source for Ims,tot . Because pn(Θ) determines only the normalized angular distribution of the irradiance, the irradiance must be determined. Here, the procedure to estimate Ims,tot,on is first described (see subsection 220.127.116.11.) Then Ims,tot,off is estimated from Ims,tot,on (see subsection 18.104.22.168).
22.214.171.124 Estimation of Ims,tot,on
Analogously to Iss, Ims,tot,on is estimated in a similar form with the single scattering approximation in Eq. (16) as follows:Eq. (19), ΔIssFon and σext,on are used instead of Io and σext in Eq. (16), respectively. Here, we define Fon and σext,on.
- (1) Fon: fraction of the source term
The fraction of ΔIss scattered in the range 0<Θ<Θon(tj) as in Fig. 4 assures that the scattered irradiance enters layer zj at tj. Using p1 and Θon, Fon is estimated by:
- (2) σext,on: effective extinction
The trajectory of the photons contributing to Ims,tot,on does not significantly deviate from those contributing to Iss. To estimate σext,on(zj), it is assumed that the extinction coefficient at zj is reduced from the single scattering value σext(zj) by the amount of photons scattered at zj-1 into the effective scattering volume at zj  as follows:
σext(z1) and p1 in Eq. (21) are treated as the scaled extinction coefficient and geometrically truncated phase function, respectively . The number of scattering events at tj equals j from the definition of lj, 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:
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 rfov(zj) and lj 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 zj 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π pj(Θ)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).
126.96.36.199 Estimation of Ims,tot,off
Ims,off (zi, tj) (i≤j) is the sum of Bon(tk) that is scattered at z1≤zk≤zi. It is also approximated in analogy to Eq. (16) as follows:
- (1) Foff: fraction of the source term
Foff(Θoff,k→i (tj)) is the fraction of Ims,tot,on that are scattered at (zk, tk) in the direction Θoff,k→i (tj-1) < Θ < Θoff,k→i (tj) defined in Eqs. (10) and (11) and incident at (zi, tj). It is estimated as in Eq. (20) as follows:
Cf is an additional term to Eq. (20). Cf = 1 for tj≠ti while Cf becomes proportional to the asymmetry factor gi for tj = ti when the FOV is small (θd,k ≤ 〈θ〉k) and the source term for Ims,off decreases.
- (2) σext,off (zq, zk, zi, tj): the decay rate from zk to zi.
Here the procedure to estimate the loss due to scattering from layer k to layer i is described. Bon in Eq. (19) already includes the loss along the path from layer 1 to k-1. σext,off at zq (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 〈nq〉 on the basis of the path integral approach , which is expressed by the 1/Nq term in the following equation:
Nq is defined as follows:
The mean z position of the scattered photons at tq+1/2 is as in  and 〈nq〉 is defined by,
Equation (28) is rounded up to an integer. Because  considered semi-infinite geometry, in this study the dependence of σext,off on τ were incorporated into the coefficients Co,i, C1,i, and C2,i in Eq. (27) as:
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 (Co,i = C1,i = 0) and becomes effective at large τi. With regard to C2,i, C2,i = 1 at τi = 1, and C2,i~0 and C2, i~gss 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(zq, zk, zi, tj) in Eq. (24) is the excessive effective photon path length from the geometrical layer thickness lq at layer q and Δd(zq, zk, zi, tj) is obtained by,
Here, the total excessive path length of the photons from layer k to layer i, (tj-ti)c, is divided equally among the number of layers lying between them, and then multiplied by the asymmetry factor of the q-th order gq. When the source term and Ims,off is from the same layer (k = i), the Nq and gq terms in Eqs. (27) and (31) are set to be 1, respectively, and σext,off = σext,on.
- (3) R(zi): backscattered ratio
Some fraction R(zi) of Ims,tot,on incident at layer i are backscattered within the angle π-θd,i ≤θ<π at tj. R(zi) are considered to remain within the FOV, i.e.,
The photons within the FOV that are collected by the receiver are determined by σsca(zi)pi(π)/(4π) in Eq. (24), where pi(π) is the fraction of photons moving in the direction π from the z-axis after the i-th scattering event.
In summary, Itot(tj) can be obtained by substituting Iss, Ims,tot,on and Ims.off estimated by Eqs. (16), (18) and (24) into the following equation,Eq. (33) is obtained through Eq. (12) as,
In the PM, the single scattering properties are calculated for the given water cloud reff and LWC profiles by Mie theory, and βss(tj), βms,tot,on(tj) and βms,off(zi,tj) are then estimated to obtain βtot(tj).
3. Evaluation by a Monte Carlo simulation
The Backward Monte Carlo method  (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 , and has been used to analyze the CALIPSO lidar signal . 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 -zsat = 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/m3 and reff = 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/m3 (σext = 3km−1) and 0.255g/m3 (σext = 40km−1).
- ・ Case 2: inhomogeneous cloud particle size profiles with reff = 20µm, 5µm, and 10 µm from the cloud top bounded at 240m and 480m with LWC = 0.1 g/m3.
- ・ Case 3: an inhomogeneous case with two cloud layers with reff = 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 × 108 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(t1) + βms,tot,on(t1) + βms,off(z1,t1)(=0), while the other purple circles correspond to βms,off(z1,tj>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 reff 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 pss. It is noted that instead of the exact pss, the OFGA method parameterized pss by double Gaussians, and here the option for the representation of the anisotropic phase function in the near-180 degree direction  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 log10βtot,PM and log10βtot,MC in Fig. 10, which was forced to pass through the origin, had a slope s = 1.008 and correlation coefficient r2 = 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.
Figure 9(b) shows comparisons among the three methods for reff = 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−1case, βtot,OFGA underestimated βtot at the 2<τ<8 regime, which was also found in Case 1. (100 × )ERRPM 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 pss change vertically. The n-th order phase function for inhomogeneous case is estimated analogously to Eq. (15) as follows:
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 pss 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., pss is in the Rayleigh regime, to investigate the applicability of PM when pss was nearly isotropic despite the number of scattering events, but sufficient multiple scattering effects exists  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, reff and FOVs, i.e., σext = 3km−1, 15.7km−1, 40km−1, reff = 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 〈ERRPM 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.
〈ERRPM〉 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 reff. On the other hand, 〈ERROFGA〉 was 65.5% on average and had the largest error in the σext = 3km−1 case. 〈ERROFGA〉 also had larger standard deviation than 〈ERRPM〉. For the same experiment setting, both 〈ERRPM〉 and 〈ERROFGA〉 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 〈ERRPM 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.
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 pn 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, reff 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  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 . 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 pn 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 .
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).
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]
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]
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]
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.”