Laser light propagation in a turbid medium: solution including multiple scattering effects

We have shown that solutions to the radiative transfer equation for a homogeneous slab yield a zenith radiance reflectance that for collimated beam incidence in the nadir direction can be expressed in terms of the lidar ratio, defined as the extinction coefficient divided by the 180∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^\circ $$\end{document} backscattering coefficient. The recently developed QlblC method, which allows one to quantify layer-by-layer contributions to radiances emerging from a slab illuminated with a collimated beam of radiation, was used to show explicitly that in the single-scattering approximation the attenuated backscatter coefficient estimated by the new QlblC method gives the same result as the lidar equation. Originally developed for the continuous wave (CW) lidar problem, we have extended the new QlblC method to apply to the pulsed lidar problem. A specific example is provided to illustrate the challenge encountered for ocean property retrievals from space observations due to the fact that a very significant fraction of the signal is due to aerosol scattering/absorption; typically only about 10% (or less) comes from the ocean.


Introduction
The United States National Academy of Sciences, Engineering, and Medicine 2017 Decadal Survey, 'Thriving on Our Changing Planet: A Decadal Strategy for Earth Observations from Space' calls for a lidar and polarimeter to accurately characterize vertically resolved absorbing and scattering properties of aerosols. The lidar is expected to be at least as capable as the Cloud Aerosol Lidar with Orthogonal Polarization (CALIOP) instrument, which has been shown to enable retrievals of aerosol, cloud, and ocean products from space. In addition, the next generation of airborne lidars includes high spectral resolution lidars (HSRL) with the capability for aerosol measurements at 1064, 532, and 355 nm [1] and ocean measurements at 532 and 355 nm [2]. These advanced, powerful lidar systems require new algorithms in order to accurately and efficiently account for multiple scattering in cirrus and water clouds as well as in open ocean waters and coastal waters with high levels of particulate scattering. For spaceborne lidar systems, multiple scattering leads to a significant contribution to lidar returns [3].
Although Monte Carlo (MC) algorithms exist for elastic backscatter lidar like CALIOP, it is highly desirable to have access to an efficient yet accurate forward radiative transfer model that will be significantly faster than existing MC algorithms in simulating lidar signals due to multiple scattering in cirrus and water clouds as well as coastal water.
Reliable interpretation of lidar/radar returns from the atmosphere (or lidar returns from the ocean) relies on accurate solutions of a forward/inverse problem. To solve the inverse problem, it is very useful to have access to a fast yet accurate forward model. For continuouswave lidar/radar beam illumination, one needs to solve the steady-state (time-independent) radiative transfer equation (RTE) for beam propagation including multiple scattering throughout the medium (atmosphere or ocean). However, most active remote sensing techniques rely on illumination using lidar/radar pulses, implying that one must solve the time-dependent RTE. For a space-based lidar, such as the CALIOP instrument, which has operated in space since 2006 on the Cloud Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) platform, one needs to solve the time-dependent RTE for the coupled atmosphereocean system in order to infer water parameters from space [4][5][6].
Many retrieval algorithms ignore or treat multiple scattering in an approximate manner, which may yield unreliable results. For example, it has been shown that if lidar multiple-scattering effects were omitted in combined radar-lidar retrievals of ice clouds from space, the retrieved optical depth might be underestimated by as much as 40% [7]. Therefore, to fully understand and exploit the information contained in received lidar/radar returns, one should have access to accurate methods to assess the error incurred by ignoring or not including multiple scattering effects in an accurate manner. Exploring to what extent multiple scattering can be used as an additional source of retrievable information about medium properties is also an important issue that deserves serious attention [8].
Radiative transfer involving lidar/radar (finite) beam illumination is a three-dimensional (3-D) problem. The solution of the 3-D RTE for a narrow finite laser beam (i.e., the so-called searchlight problem) is quite challenging and computationally demanding. Therefore, it has become customary to use a one-dimensional (1-D) approach instead, and most treatments of the lidar/radar problem rely on solving a 1-D RTE for both atmospheric [9,10] and oceanic [11] applications.
The advantage of solving the RTE in a slab geometry, which leads to the SSA result in Eq. (1), is that it yields valid results also when multiple-scattering effects prevail. However, a slab geometry, on which Eq. (1) is based, does not apply to a laser beam of finite lateral extent. Nevertheless, Eq. (1) is expected to be a reasonable approximation when the scattering mean free path is less than the lateral extent of the laser beam.
The Quantification of layer-by-layer Contribution (QlblC) method [12] can be used to estimate the backscatter attenuation coefficient for the lidar problem as outlined in Sect. 3

The QlblC method
Let's consider a medium consisting of two adjacent, horizontal, multilayered coupled slabs illuminated from the top of the upper slab by a collimated beam of irradiance F 0 in direction θ 0 with respect to the normal to the two adjacent slabs. As discussed elsewhere [12], the radiative transfer equation (RTE) for such a problem can be solved by integrating the source function S ± i (t, μ, φ) (i = n, or p (layer indices) in Eqs. (2) and (3)) layer by layer to yield the following expressions for the Stokes vector I ± (τ, μ, φ) of the diffuse total and polarized radiance (the + sign denotes the upward hemisphere while the − sign denotes the downward hemisphere): Here τ is the optical depth, μ = (cos θ) (θ is the polar angle), and φ is the azimuthal angle. Since the source function S ± i (t, μ, φ) in layer denoted by subscript i (= n, or p) can be evaluated analytically by the discrete ordinate method, the integrals in Eqs. (2) and (3) have analytic solutions.
Note that Eqs. (2) and (3) allow us to quantify the contributions from any specific layer in a single slab or any layer in either of two adjacent slabs with different refractive indices. For example, if we are interested in the contribution from layer M in either of the two adjacent slabs to the zenith radiance emerging from the top of the upper slab, we may proceed as follows: Alternatively, we may obtain the same result by simply setting the source function to zero in all layers except in layer M , which is the layer of interest.

A simple two-layer model
As a simple example, consider a single slab (as opposed to two adjacent, coupled slabs) consisting of two layers (L = 2): a "target" (layer 2) overlain by layer 1. The slab is assumed to be illuminated by a collimated beam at an angle θ 0 . For this two-layer system, the reflected signal would be given by:

Single-scattering approximation (SSA)
Assuming that the single-scattering approximation is applicable for each of layers 1 and 2, we have (see Eq. (A3)) where we have ignored polarization (S + i → S + i ) and μ 0 = cos θ 0 . If we consider both layers to be homogeneous and to scatter in accordance with the Henyey-Greenstein scattering phase function, p HG (−μ 0 , μ, g i ) with asymmetry factor g i , then or (setting τ 0 = 0 at the top of the upper slab) is the lidar ratio, defined as the extinction coefficient divided by the 180 • backscattering coefficient, and g i is the asymmetry factor of the scattering phase function, the only input parameter required for the HG phase function. Clearly, if the two layers have identical optical prop- as expected for a homogeneous slab of optical thickness

QlblC estimate
To mimic the lidar problem, let's assume that both layers consist of non-absorbing particles ( 1 = 2 = 1.0) and that layer 1 scatters isotropically (with asymmetry factor g 1 = 0) and has a small optical thickness (say τ 1 = 0.05), while layer 2 (the target layer) scatters according to a Henyey-Greenstein scattering phase function with asymmetry factor g 2 and optical thickness τ 2 . According to the QlblC method, the contribution to the radiance reflectance from layer 2, the target layer, i.e., the single-scattering estimation of the radiance reflectance R I,QlblC is given by the second term in Eq. (8): The attenuated backscatter from the target layer (layer 2) is obtained by dividing by its vertical thickness Δz:

Lidar equation estimate
The attenuated backscatter coefficient (assuming a range bin of vertical extent Δz) from lidar measurements is given by the lidar equation: Setting T 2 = e −2τ1 and αΔz ≈ τ 2 − τ 1 , we find that the attenuated backscatter results predicted by Eqs. (10) and (12) are the same, i.e., β QlblC = β lidar if 2 = 1.0 and S = S 2 .

Generalization to an N -layer model
Assuming that all N −1 layers above the target layer N have the same optical properties, then Eq. (10) should be replaced by and Eq. (12) by

Numerical examples
The QlblC method [12] allows one to quantify specific layer contributions of signals emerging from stratified, multilayered media based on solutions to the radiative transfer equation (RTE). As shown in the previous section, in the single-scattering approximation, the QlblC solutions agree with those obtained using the standard lidar approach, but the QlblC method also yields reliable results for optically thick, multiple-scattering aerosol and cloud layers. In the following, we use radiation reflected from (i) a multilayered atmosphere and (ii) a coupled atmosphere-ocean system to address the problem relevant for EM probing by a space-based lidar system. Hence, we will assume that the medium consists of either (i) one single, multilayered slab with a constant refractive index or (ii) two adjacent multilayered slabs separated by a smooth interface across which the real part of the refractive index m r changes from one constant value m r,1 in slab 1 to a different constant value m r,2 in slab 2 , such as for an atmosphere-water system. To deal with the vertical variation of the scattering and absorption properties of the atmospheric constituents, it is necessary to treat it as a multilayered medium. Therefore, based on previous experience, we divided the atmosphere into 15 layers from top to bottom as illustrated in Fig. 1. This vertical stratification of the atmosphere in horizontal layers is sufficient to produce accurate radiances for a clear (cloud-and aerosol-free) atmosphere.
The key to accomplishing this quantification of the emerging radiance originating from any specific layer (or location) in the medium lies in our ability to compute and integrate the source function. 1 Thus, we are able to quantify the contribution to the emerging radiance from any particular layer of a stratified medium. Below we consider three different cases: 1. A clear atmosphere with embedded aerosol layers. 2. Same as item 1 above but contrasting continuous wave (CW) and pulsed lidar results.  3. An atmosphere overlying a body of water.

A clear atmosphere with embedded aerosol layers: CW lidar
First, we consider the CW lidar problem in which the steady-state RTE is solved, implying that the source function for any particular layer includes multiply-scattered radiation from all layers below the layer of interest. Figure 2 shows the layer-by-layer (lbl) contributions to the reflected radiance at the top of the atmosphere (TOA) for wavelength λ = 532 nm for this CW lidar configuration. As expected for an atmosphere that obeys the barometric law of near-exponential decrease with altitude in the bulk molecular number concentration, the major contributions to the reflected radiance at the TOA come from the layers close to the surface. To mimic the lidar problem, we assumed nadir incidence, i.e., a beam zenith angle of 0 • , and a polar observation angle of 0 • . To demonstrate how aerosols impact the lbl contributions to the TOA radiance, we included an aerosol component in each of layers 14 and 15 (closest to the surface). This aerosol component was assumed to consist of weakly absorbing particles, homogeneously mixed with the molecules in layers 14 and 15. The lbl contributions to the TOA radiance for this case are shown in Fig. 3. Note that the largest lbl contribution to the TOA radiance comes from layer 15, and the next largest from layer 14 as expected since the aerosol optical thickness in layer 15 is twice that of layer 14.

Continuous wave (CW) versus Pulsed lidar
The results shown in Figs. 2 and 3 pertain to a CW lidar configuration. To address the pulsed lidar problem, we adopt the following approach: • In the CW lidar problem: the source function for each layer contains contributions also from all layers, including those below a given layer of interest. • To simulate the pulsed lidar problem, we must omit contributions from layers below the layer of interest. Hence, we may proceed as follows: 1. Starting at the top of the slab (atmosphere), we first solve a two-layer problem to get the lbl contribution from the second layer (counting from the top). The uppermost layer is optically very thin. Hence, we can ignore its possible contribution, and contributions from all layers below layer 2 are ignored. 2. Next, we solve the 3-layer problem to get the lbl contribution from layer 3, again ignoring contributions from all layers below layer 3. 3. Then, we continue in this manner until we have determined the lbl contributions from layer 4, layer 5, . . ., layer N for a N -layer slab.
Results for the CW and pulsed lidar configurations are shown in Fig. 4. Note that the pulsed lidar configuration (right panel) yields smaller lbl contributions than the CW lidar configuration (left panel) as expected.

An atmosphere overlying a body of water
Our final example pertains to a clear atmosphere overlying a body of oceanic water. The atmosphere is assumed to consist of 15 layers, and the ocean is assumed to consist of three layers with a total thickness of 100 m. Each of the two uppermost ocean layers has a thickness of 5 m and contain pure water mixed Below 10 m depth, the water is assumed to be pure. The left panel of Fig. 5 shows lbl contributions. Note that the largest contribution from the ocean comes for the upper layer (layer 16). Layer 17 (below) is effectively shielded by the overlying layer 16, and very little light penetrates to layer 18. The right panel of Fig. 5 shows the same results as in the left panel, but presented as reflectance profiles as is common in the atmospheric lidar community. Please note that in Figs. 5 and 6 the percentage contributions are per layer (not per km as in Figs. 2, 3, 4). Figure 6 shows the impact of adding aerosol particles in layers 14 and 15 at the bottom of the atmosphere. In the upper left panel of Fig. 6, no aerosol particles were added, but in the upper right panel a small amount (AOD = 0.003) was added. We note that this small amount has almost negligible impact. The contribution from the water to the TOA radiance is about 10.3%, while that from the two aerosol layers is about 17.6%. As the AOD is increased to 0.03 (lower left panel), those contributions change to about 10.2% (ocean) and 16.6% (aerosols). Finally, the lower right panel shows that for a total AOD abundance of 0.3 in layers 14 and 15, these contributions change to 4.6% (ocean) and 52.1% (aerosols). These numbers are in general agreement with observations, and they illustrate that in order to retrieve reliable results for ocean properties from optical instruments deployed on space platforms, the impact of aerosols must be dealt with very carefully, since their abundance change in space and time.

Conclusion
We have shown that solutions to the RTE for a homogeneous slab of optical thickness τ b yields a zenith radiance reflectance that for collimated beam incidence in the nadir direction can be expressed in terms of the lidar ratio S simply as R I,SSA = 1 2S (1 − e −2τ b ) in the single-scattering approximation. The lidar ratio, defined as the extinction coefficient divided by the 180 • backscattering coefficient, is given by S = S iso = 4π for isotropic scattering, S = S Ray = 8π 3 for Rayleigh scattering, and S = S HG = 4π(1+g) 2 (1−g) for Henyey-Greenstein (HG)scattering, where g = cos Θ is the scattering asymmetry factor and Θ is the scattering angle.
A simple two-layer model was used to show explicitly that in the single-scattering approximation an estimate of the attenuated backscatter coefficient based on the QlblC method gives the same result as the lidar equation. As described in a recent paper [12], the QlblC method applies to the continuous wave (CW) lidar problem. In Sect. 3.5.2, we extended the QlblC method to be applicable to the pulsed lidar problem and showed that as expected the pulsed lidar configuration yields smaller lbl contributions than the CW lidar configuration.
In Sect. 3.5.3, an example was provided (Figs. 5 and 6) to illustrate the challenge encountered for ocean property retrievals from space observations due to the fact that a very significant fraction of the signal is due to aerosol scattering/absorption; typically only about 10% (or less) comes from the ocean. The contribution from molecules is easier to deal with because it depends largely on the total column amount of molecules which depends on pressure through the ideal gas law. Pressure measurements are commonly available.