A novel approach to solve forward/inverse problems in remote sensing applications

1 Introduction

1.1 Background

The National Academy of Sciences, Engineering, and Medicine 2017 Decadal Survey Luvall et al. (2017) calls for a lidar and polarimeter system to accurately characterize vertically-resolved absorbing and scattering properties of aerosols. To meet the requirements called for, the lidar should be at least as capable as the Cloud Aerosol Lidar with Orthogonal Polarization (CALIOP) instrument (Winker et al., 2003; Powell, 2005), which has demonstrated the retrieval of cloud (Hu et al., 2001, 2007) and ocean (Hu, 2009; Behrenfeld et al., 2013, 2017) products from a spaceborne lidar system. The next generation of airborne lidars, include high spectral resolution lidars (HSRL) with the capability for aerosol measurements at 1,064, 532, and 355 nm (Burton et al., 2014) and ocean measurements at 532 and 355 nm (Hair et al., 2016). These advanced, powerful lidar systems require new algorithms in order to accurately and efficiently account for multiple scattering by cirrus and water clouds, as well as by high concentrations of particulate matter in coastal waters. For spaceborne lidar systems, multiple scattering has been noted to contribute significantly to lidar returns, implying deeper penetration into opaque layers (Thorsen et al., 2013). Although Monte Carlo algorithms exist for elastic backscatter lidars like CALIOP, it is desirable to develop an efficient yet accurate forward radiative transfer model that will also be capable of simulating multiple scattering-induced signals measured by HSRL for media such as water clouds, cirrus clouds, and the ocean. By design such a model should also include the full effects of polarization, and it is desirable to develop a fully deterministic model, that will be significantly faster than existing Monte Carlo algorithms. Such a model can then be used as a forward model for optimal estimation retrievals of atmospheric and ocean properties. Platnick et al. (2017) report a water cloud mean optical depth of about 5. Hence, lidar has the potential to retrieve significant information about cloud properties, particularly for low-lying marine clouds which can be difficult targets for radar. The capability of lidar to measure aerosol properties below thin clouds and to retrieve ocean properties will also be made possible with an accurate and efficient multiple-scattering treatment.

1.2 Main purpose

Active remote sensing of an absorbing and scattering medium by electromagnetic (EM) beam probing (using e.g. X-ray, radar, or lidar beam illumination) has been used extensively in a great variety of research applications (Nicolas et al., 1997; Mitra and Churnside, 1999; Hu et al., 2001; Hu, 2007; Hogan, 2008; Behrenfeld et al., 2013, 2017; Burton et al., 2014; Hair et al., 2016; Hostetler et al., 2018). Similarly, determining the radiation emitted from a medium due to internal sources (e.g. due to thermal emission) can be used in passive remote sensing to investigate the internal structure of the medium, like the temperature profile of a non-isothermal atmosphere (Rodgers, 2000).

The main purpose of this paper is to describe how to quantify specific layer contributions either to the reflected or transmitted EM signal emerging from a multilayered, stratified medium that is illuminated by a beam of EM radiation or to the emitted EM signal emerging from such a medium due to emission by internal sources. To simplify the discussion we will assume that the medium consists of either 1) one single, multilayered slab with a constant index of refraction or 2) two adjacent multilayered slabs (labeled slab₁ and slab₂) separated by a smooth interface across which the real part of the refractive index changes abruptly from one constant value (m_r1) in slab₁ to a different constant value (m_r2) in slab₂, such as for an atmosphere-water. Note that we use the word “water” generically to describe the solid phase (i.e., snow and ice) as well as the liquid phase. Hence, “water” could mean ocean, sea ice, snow, or a fresh liquid water body. 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 (see Section 2.1), because knowledge of the source function is equivalent to knowledge of the complete solution of the radiative transfer problem (Stamnes et al., 2017).

1.3 Brief history

The transport of radiation through an absorbing, emitting, and multiple-scattering stratified medium is described by the radiative transfer equation (RTE), which is an integro-differential equation with no known analytic solutions. A variety of techniques are available to solve the RTE for such a stratified medium including the discrete ordinate method (Stamnes and Stamnes, 2015; Stamnes et al., 2017). One advantage of the discrete ordinate radiative transfer (DISORT) method is that an analytic solution for the source function can be obtained in terms of the approximate solutions valid at the quadrature points—the so-called “discrete ordinates” or computational polar angles. From this analytic solution for the source function an analytic expression for the radiance at any arbitrary polar angle can be derived. This approach, known as the integration of the source function (ISF) method (Stamnes, 1982), is implemented numerically in the standard DISORT code (Stamnes et al., 1988) applicable to the scalar version of the RTE. The latest version of DISORT, which includes significant upgrades and improvements of both its speed and accuracy, is described by Lin et al. (2015).

A discrete ordinate solution to the vector version of the RTE that accounts for polarization effects was presented by Weng (1992). Significant improvements to this vector (VDISORT) version were reported by Schulz et al. (1999), and Schulz and Stamnes (2000). The ISF method has also been implemented in this VDISORT code (Lin et al., 2022) pertinent for a single slab medium with a constant index of refraction.

Yan and Stamnes (2003) extended the ISF method to apply to the scalar C-DISORT code pertinent for two coupled adjacent slabs with different refractive indices, like an atmosphere-water system (Jin and Stamnes, 1994). The C-DISORT ISF method has been implemented in the AccuRT platform designed for such coupled systems (Stamnes et al., 2018). Also, a Coupled Vector Discrete Ordinate (C-VDISORT) RTM was developed for pure Rayleigh scattering (Sommersten et al., 2010) as well as for scattering by particles in the Mie regime (Cohen et al., 2013).

2 Theoretical background

A stratified inhomogeneous medium (such as the Earth’s atmosphere, or a body of water) can be represented as a layered medium, i.e. as a multilayered slab of L layers in which each layer has constant inherent optical properties (IOPs), which may vary from layer to layer. The number of layers L should be large enough to adequately resolve the variation of the IOPs. If a multilayered medium is illuminated by a beam of infinite extent in the lateral direction normal to that of the stratification, we have a 1-dimensional (1-D) radiative transfer problem.

2.1 Solutions to the radiative transfer equation for a 1-D problem

2.1.1 The integration of the source function (ISF) method

As shown in Section 7, Supplementary Appendix S1, see Supplementary Appendix Eqs. 35, 36, the RTE can be solved by integrating the source function ${\mathbf{S}}_i^{\pm }(t,\mu ,\varphi )$ (i = n, or p in Eqs. 1, 2 layer by layer to yield the following expressions for the Stokes vector I^±(τ, μ, ϕ) of the diffuse total and polarized radiance (the + (plus) sign denotes the upward hemisphere while the − (minus) sign denotes the downward hemisphere):

In Eqs. 1, 2, which are obtained using the integration by the source function (ISF) method, t is a ‘dummy’ (integration) variable τ is the optical depth, μ = | cos θ| (θ is the polar angle), and ϕ is the azimuthal angle. Since the source function ${\mathbf{S}}_i^{\pm }(t,\mu ,\varphi )$ in layer denoted by subscript i (= n, or p) can be evaluated analytically by the discrete ordinate method, the integrals in Eqs. 1, 2 also have analytic solutions. Such solutions are valid for a single slab as well as for two adjacent slabs separated by an interface across which the refractive index is allowed to change (like in an atmosphere-water system). Evaluating Eq. 1 at τ = 0 and μ = 1 (ϕ irrelevant), we obtain the upward total and polarized radiance in the zenith direction (relevant for spaceborne lidar deployments), while evaluating Eq. 2 at τ = τ_L and μ = 1 (ϕ irrelevant), we obtain the downward total and polarized radiance in the nadir direction (relevant for ground-based lidar deployments).

2.1.2 The essence of the QlblC method

Note that the ISF method Eqs. 1, 2 allows 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. This quantification is the essence of the QlblC method. 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:

• Compute I⁺(τ = 0, μ = 1) by integrating all layers from the top of the upper slab including layer M;

• Repeat the computation by integrating all layers from the top of the upper slab including layer M − 1;

• Subtract the latter result from the former to obtain the contribution from layer M.

Alternatively, (and more efficiently) we can apply the ISF method to obtain the same result by simply setting the source function to zero in all layers except the one we are interested in. We should emphasize that the QlblC method consists simply of using the ISF to quantify the contribution by individual layers to the signal reflected from or transmitted through the medium. By adding the contributions from all layers we obtain the total signal as required.

Note that Eqs. 1, 2 are not new (see Stamnes and Stamnes (2015); Stamnes et al. (2017) and references therein), but using them to quantify contributions to the radiance emerging at the top or bottom of a scattering/absorbing medium from any specific layer in the medium is a novel idea, behind the QlblC method. Although this idea may seem obvious in retrospect, we are not aware of any previous use of this approach to quantify the contribution from any specific layer in a medium to the emerging radiances. As we demonstrate in Section 4.3, this approach allows one to compute attenuated backscatter profiles for the lidar problem appropriate for space-based lidar systems such as the CALIOP lidar deployed on the CALIPSO satellite (Winker et al., 2003). In Section 4.5 we show that the results agree well with attenuated backscatter profiles derived by independent Monte Carlo (MC) simulations. The discrete ordinate solutions employed in AccuRT are much faster than MC simulations and they provide an efficient and accurate QlblC solution to the lidar problem for space-borne (e.g. CALIPSO and IceSat) lidar systems as well as for ground-based lidar systems.

3 Application to the lidar problem

3.1 The lidar equation

Consider an air-water system that is illuminated by a laser beam of finite lateral extent. In the case of a pulsed light beam, time-gating can be used to keep track of the depth (or layer) in the medium from which the observed backscattered radiance originates. A simple solution to the observed backscattered radiance is obtained by assuming that it is due to a single backscattering event in the layer of interest (target layer) multiplied by one exponential attenuation of the signal going from the source to the target (before the single backscattering event) and another exponential attenuation going back from the target to the observer (after the single backscattering event). Here backscattering means that the signal reverses direction by 180°. Hence, in the conventional lidar approach the backscatter contribution ${\widehat{\beta }}_{\text{b,\hspace{0.17em}lid}}(r)$ coming from a target layer located at range r is given simply by

Here τ(r) is the optical path of the medium between the lidar and a target (scattering layer) at range r, T²(r) is the two-way transmittance, and β_{b, lid}(r) is the backscattering coefficient [m⁻¹sr⁻¹] of the target layer located at range r. For a vertically stratified medium like the atmosphere-water system illuminated by a lidar beam normal to the stratification, the range is simply the vertical distance $\mathcal{L}+z$ between the lidar and the scattering (target) layer of thickness z located at depth $\mathcal{L}$, and the optical path τ(r(z)) is replaced by the vertical optical depth $\tau (z)={\int }_0^{\mathcal{L}+z}\gamma (z^{\prime })d{z}^{\prime }$, where γ(z) = α(z) + β(z) [m⁻¹], with γ(z), α(z), and β(z) being depth dependent extinction, absorption, and scattering coefficients of the medium between the lidar source and the target.

The optical thickness of the scattering (target) layer is $\mathrm{\Delta }\tau =\gamma (\mathcal{L})z$.

3.2 Radiative transfer

The Radiative Transfer Equation (RTE) involving lidar (finite) beam illumination describes a 3-dimensional (3-D) problem, and 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, most treatments of the lidar problem rely on solving a 1-dimensional (1-D) RTE [i.e. Supplementary Appendix Eq. 22 in Section 7] for both atmospheric (Hogan, 2008; Hogan and Battaglia, 2008) and aquatic (Mitra and Churnside, 1999) applications.

3.2.1 Single-scattering approximation

The single-scattering approximation is valid when the single-scattering albedo, ϖ, and the optical thickness of the medium, Δτ, satisfy the condition ϖΔτ ≪ 1. In the single-scattering approximation we ignore the multiple-scattering term in Supplementary Appendix Eq. 22, which then permits analytic solutions of the RTE. If the beam is incident perpendicularly on a homogeneous slab of optical thickness τ_b, then [see Supplementary Appendix S2, Section 8, Supplementary Appendix Eq. 48] for isotropic (conservative) scattering (ϖ = 1.0) the radiance reflectance (in the upward direction perpendicular to the slab) is given by (τ_b ≪ 1)

where g is the scattering asymmetry factor. For conservative scattering according to a Henyey-Greenstein phase function one gets [see Supplementary Appendix Eq. 51]

3.2.2 Multiple scattering

We now assume that the propagation of one particular pulsed laser beam is separated in time from the next pulsed laser beam so that there is no time overlap between them. Then we can focus attention on the fate of the energy deposited by this single laser beam as it propagates through the medium. This energy deposition can be determined by solving the steady-state 1-D RTE [i.e. Supplementary Appendix Eq. 22], but how can we quantify the contribution to the observed backscattered radiance that originates at a specific target layer in the medium? If we know the distance from the source (laser transmitter) to the target layer of interest, we know (through time-gating) when the observed signal will arrive at the receiver, and if the two-way transmittance represented by T²(r) in Eq. 3 applies, we have an analytic solution. But how do we determine the effective backscattering coefficient

i.e. the factor C₁ multiplying β_{b, lid}(r) in Eq. 6 if the medium surrounding the target layer giving rise to the observed backscattered radiance does not satisfy the assumption of a single backscattering event in the target multiplied by the two-way transmittance factor T²(r) in Eq. 3 based on Beer’s law attenuation of the incident beam from the source to the target and back from the target to the detector?

The answer to this question is that the solutions I^±(τ, μ, ϕ) given by Eqs. (1), (2) allow us to quantify the contribution from individual layers to the observed signal as explained in Section 2.1. The evaluation of the source function in Eqs. 1, 2 does not rely on the assumptions invoked to arrive at the lidar equation [Eq. 3], implying that the contribution from any layer of interest includes all orders of multiple scattering.

In principle, one could divide the medium into a large number of vertically stratified horizontal layers so that for each optically thin layer the single-scattering criterion is satisfied. Then an analytic single-scattering solution could be applied for each of these layers, but since the total optical depth for the entire system will be large the single-scattering approximation would not be applicable.

3.3 Assessment of the impact of multiple scattering

As mentioned above, if the single-scattering albedo ϖ and optical thickness Δτ of the medium satisfy ϖΔτ ≪ 1, then the single-scattering approximation is valid. Otherwise, multiple scattering cannot be ignored. When multiple scattering occurs the backscattered light appears to be delayed compared to the singly scattered light.

An approximate estimate of the importance of multiple scattering is obtained by considering the ratio $C_1={\beta }_{\text{b,\hspace{0.17em}lid}}^{\prime }(r)/{\beta }_{\text{b,\hspace{0.17em}lid}}(r)$ [Eq. 6] between the multiple scattering result (obtained for the actual slab optical thickness) and the single-scattering result (obtained for an optically thin medium, i.e. when the condition ϖΔτ ≪ 1 is satisfied). The tail (delayed return) observed experimentally between multiple-scattering results and single-scattering results should be approximately proportional to this ratio C₁. The delayed return should be proportional to the increased pathlength due to multiple scattering. In the case of single scattering event, the mean free path of a light beam (i.e. the mean distance traveled by a photon between collisions) can be estimated as ${\ell }_{\text{mfp}}\propto \frac{1}{\beta }$ where β [m⁻¹] is the scattering coefficient. In the presence of multiple scattering, ℓ_mfp should be replaced by the transport mean free path, which can be estimated as ${\ell }_{\text{trans}}\propto \frac{1}{(1-g)\beta }$, where g is the scattering asymmetry factor (i.e. the first moment of the scattering phase function). Hence, a time delay between multiple-scattering and single-scattering results is expected because 0 < g < 1 and therefore ℓ_trans > ℓ_mfp.

The travel time t_ss for singly scattered photons (no multiple scattering) can be estimated as ℓ_mfp ∝ 1/β = t_ssv_c, where v_c is the speed of light, while the travel time t_ms for multiply scattered photons can be estimated from ℓ_trans ∝ 1/(1 − g)β = t_msv_c. These estimates imply that t_ms = (ℓ_trans/ℓ_mfp)t_ss = t_ss/(1 − g). Hence, the time for multiply scattered light to reach the sensor would be delayed by approximately a factor 1/(1-g) compared to the travel time for singly scattered light.

4 Simulations relevant for the lidar problem

We start by providing an overview of the lidar problem from the point of view of a space-borne system designed for atmospheric applications (Section 4.1). Then in Section 4.2 we give a brief introduction to AccuRT, our forward radiative transfer model. In Section 4.3 we provide three examples of how to quantify layer-by-layer contributions to the reflected radiance for 1) a Rayleigh scattering atmosphere (Section 4.3.1), 2) a Rayleigh scattering atmosphere with an embedded aerosol layer (Section 4.3.2), and 3) a Rayleigh scattering atmosphere with an embedded cloud layer (Section 4.3.3). In Section 4.4 we discuss the impact/importance of multiple scattering, and finally, in Section 4.5 we present comparisons with Monte Carlo simulations.

4.1 Basic description of the lidar problem for a single slab atmospheric system

We will now consider a spaceborne lidar system like the CALIOP lidar deployed on the CALIPSO satellite (Winker et al., 2003). To be specific we will assume that the lidar emits a pulsed beam in the nadir direction, and receives backscattered radiation in the zenith direction. The receiver is assumed to be placed in close proximity to the transmitter on the same platform.

Also, the transmitter is assumed to emit a lidar pulse of EM radiation in the nadir direction and the receiver is assumed to detect the radiation backscattered from the target in the zenith direction.

The radiance P(r(z), λ) that is backscattered to a receiver on a spacecraft situated at an altitude z above the mean sea level from a sample volume of atmosphere at range r is given by (Powell, 2005)

where P_Bkg(r(z), λ) is the background radiance and P_Bsc(r(z), λ) is the lidar radiance backscattered to the observer (receiver). The background radiance is due to molecular (Rayleigh) scattering caused by sunlight illumination, while the lidar radiance is described by the lidar equation (strictly valid only in the single-scattering limit):

Here C₂(λ) is a calibration coefficient (counts⋅m³) which contains the characteristics of the lidar transmitter. The attenuated backscatter profile ${\beta }_{\text{b}}^{\prime }(z,\lambda )$ (m⁻¹ sr⁻¹) describes the scattering properties of the atmospheric (molecular, aerosol, cloud) components. It is given by

The total backscattering coefficient due to altitude dependent molecular and particulate (aerosol and cloud) contributions, given by

specifies the 180° lidar radiance backscattered to the receiver. The two-way transmittance between the lidar and the sample volume is given by

where the optical depth τ(z, λ) is defined as

Here γ(z, λ) (km⁻¹) is the total extinction coefficient formed by summing the altitude dependent molecular, particulate, and ozone extinction coefficients

Note that potential absorption by trace gases other than ozone has been ignored in Eq. 13. For wavelengths other than 532 and 1,064 nm used in the CALIOP lidar system, absorption by other trace gases may be important and should be taken into account. For CALIOP the ozone absorption coefficient γ_{O_{3}}(z, λ) is negligible at 1,064 nm and is only applied to the 532 nm signals. The lidar ratio S(z, λ) (sr⁻¹) is defined as the ratio between the total extinction and total backscattering coefficients:

The CALIOP lidar measures two orthogonal polarization components of the light backscattered at 532 nm. The beam emitted by the lidar at the wavelength of 532 nm is linearly polarized in the parallel (‖) direction, but the backscattered light generally contains a perpendicular (⊥) as well as a ‖ component, although light backscattered from spherical particles remains polarized in the ‖ direction. The CALIOP polarization cube separates the backscattered light at 532 nm into its ‖ and ⊥ components, given by

where δ_x is the depolarization ratio and the subscript x refers to the molecular (x = m) or particulate (x = p) component of the backscattered light. The depolarization ratio is the ratio of the backscattered radiance in the ⊥ channel to that in the ‖ channel. At 532 nm, the backscattering profiles β_b,‖(z, λ) and β_b,⊥(z, λ) for respectively the parallel and perpendicular channels are defined as

The lidar power P(h) received by the lidar detector in space from a target layer at range $r=\mathcal{L}+h$, where $\mathcal{L}$ is the target range, i.e. the distance from the lidar source to the top of the target, and h is the penetration depth into the target (cloud or aerosol layer), can be described as follows

In Eq. 18 P₀ is the power received at the top of the target layer; K is an instrument constant; v_c is the speed of light in vacuum; τ is the pulse width; A is the collection area of the telescope; β_π is the volume backscattering coefficient at π rad (180°); k_lidar is the effective attenuation coefficient, and T_a is the transmittance through the atmosphere given by

where γ_a(h′) is the extinction coefficient of the atmosphere. According to Eq. 18 the power P(h) due to a target layer of thickness h (cloud or aerosol layer) depends on lidar hardware parameters (K, A), target range $(\mathcal{L})$, and target optical properties β_π and k_lidar. We may define the normalized power as

Hence,

For a sufficiently large field-of-view (FOV), the lidar signal returned from depth h in the target can be approximated by setting k_lidar ≈ γ, where γ is the target extinction coefficient assumed to be independent of depth h. Then we may rewrite Eq. 19 as

4.2 AccuRT—A tool for radiative transfer simulations

AccuRT is a tool developed to simulate the radiative transfer (RT) process in a coupled atmosphere-water system. It simulates irradiances as well as radiances in any user-specified direction and at any desired vertical location in the coupled system consisting of two adjacent horizontal slabs in which the real part of the refractive index m_r is assumed to have a constant value in the upper slab (m_{r, air} = 1.0) and a different constant value in the lower slab (m_{r, water} ≈ 1.34, but wavelength dependent).

AccuRT has two distinct features (Stamnes et al., 2018):

1. It deals with a coupled system that includes both the atmosphere and the underlying water body (modeled as two adjacent slabs) and it allows for a change of the refractive index across the atmosphere-water interface. Most RT models only consider either atmosphere or water and are not coupled.

2. It is very flexible, allowing the user to specify 1) location and direction of input and output, 2) inherent optical properties (IOPs) of the materials (dissolved or particulate matter) inside each layer in each slab, and 3) parameters affecting the IOPs (including particle size, density, impurities, etc.) of the materials.

The QlblC method described in Section 2.1.2 has been implemented in AccuRT, which can be used to carry out the QlblC computations.

AccuRT is a very efficient computational tool for a coupled atmosphere-water system, and its use of the integration of the source function method to quantity specific layer contributions to the radiance emerging from a scattering/absorbing medium is unique. Readers interested in the numerical code may contact the first author (KS) in order to obtain a demo version of AccuRT. The DISORT and VDISORT codes are publicly available at www.rtatmocn.com. Other methods of solutions to the radiative transfer equation could in principle be used to compute the source function which in turn could be integrated (numerically as need be) to produce results similar to those presented in this paper. The advantage of the discrete ordinate method is that such results are already accessible (in the DISORT and VDISORT codes publicly available at www.rtatmocn.com) and can be used to produce the desired layer-specific contributions to the emerging radiances by the method described in this paper.

4.3 Contributions to the backscattered radiance from a particular layer

By solving the lidar problem, one is able to quantify the contribution to the backscattered light from any particular layer of a stratified medium. Below we consider three different cases:

1. A Rayleigh scattering atmosphere (Section 4.3.1);

2. A Rayleigh scattering atmosphere with moderate aerosol loading, AOD = 0.3 at 500 nm (Section 4.3.2).

3. A Rayleigh scattering atmosphere with an opaque liquid water cloud layer (Section 4.3.3).

4.3.1 A Rayleigh scattering atmosphere

First we consider the radiance reflected from an aerosol- and cloud-free molecular US Standard atmosphere. We divided the atmosphere into 15 layers as follows from top to bottom (in km): 100–70 (30 km, layer 1), 70–50 (20 km, layer 2), 50–40 (10 km, layer 3), 40–30 (10 km, layer 4), 30–25 (5 km, layer 5), 25–20 (5 km, layer 6), 20–15 (5 km, layer 7), 15–13 (2 km, layer 8), 13–10 (3 km, layer 9), 10–8 (2 km, layer 10), 8–6 (2 km, layer 11), 6–4 (2 km, layer 12), 4–2 (2 km, layer 13), 2–1 (1 km, layer 14), and 1–0 (1 km, layer 15). Experience has shown that this vertical stratification of the atmosphere in horizontal layers is sufficient to resolve vertical variations in the IOPs required to produce accurate radiances. In fact, even in the absence of aerosol and cloud droplets the atmosphere is vertically inhomogeneous due gaseous absorption by molecular oxygen (O₂) as well as ozone (O₃) and other trace gases. Hence, layering is necessary to account for the vertical variation of the IOPs. The aerosol and cloud IOPs used in layers 14 and 15 were assumed to be vertically homogeneous in the simulations carried out in this paper. However, if partially absorbing aerosols were to be mixed with the exponentially varying density distribution of molecules, some vertical variation in the effective (combined) IOPs would result that would require a subdivision of layers 14 and 15 in order to obtain accurate radiances. For the results presented in this paper involving non-absorbing cloud droplets and weakly-absorbing aerosol particles, the 1 km layer thickness of layers 14 and 15 is expected to be adequate.

To mimic the lidar problem, we assumed nadir incidence, i.e. a beam zenith angle of 0°, and a polar observation angle of 0°. For simplicity, the results shown in Figures 1–5 pertain to a molecular atmosphere assumed to be non-absorbing. Figure 1 shows the layer-by-layer (lbl) contributions to the reflected radiance at the top of the atmosphere (TOA) for wavelength λ = 532 nm. Clearly, the major contributions to the reflected radiance at the TOA comes from the layers close to the surface as expected for an atmosphere that obeys the barometric law of near-exponential decrease in molecular number concentration with altitude.

FIGURE 1. Layer-by-layer (lbl) contributions to the reflected radiance for a non-absorbing Rayleigh scattering atmosphere. Left panel: Percentage lbl contributions. Right panel: lbl contributions to the scattering coefficient. The total Rayleigh scattering optical depth at 532 nm was 0.134.

FIGURE 2. Same as Figure 1 except that aerosol particles were mixed with the molecules in layers 14 and 15. The total AOD of layers 14 and 15 was taken to be 0.3.

FIGURE 3. Same as Figure 1 except that cloud droplets were mixed with the molecules in layers 14 and 15. The optical depth of the cloud droplets was assumed to be 2.25 in each of layers 14 and 15, giving the total cloud optical depth was 4.5.

FIGURE 4. Left panel: Same as left panel of Figure 2 except that the total AOD for layers 14 and 15 was decreased from 0.3 to 0.015. Middle panel: Same as left panel of Figure 2, but for a total AOD of 0.15 for layers 14 and 15. Right panel: Same as left panel of Figure 2 i.e. for a total AOD of 0.3 for layers 14 and 15.

FIGURE 5. Left panel: Same as the left panel of Figure 3 except that the total cloud optical depth (COD) for layers 14 and 15 was decreased by a factor 100 to 0.045. Middle panel: Same as the left panel of Figure 3 except that total COD for layers 14 and 15 was decreased by a factor 10 to 0.45. Right panel: Same as left panel of Figure 3 for a total COD of 4.5 in layers 14 and 15.

4.3.2 A Rayleigh scattering atmosphere with moderate aerosol loading—AOD = 0.3 at 500 nm

To illustrate the impact of aerosols on the lbl contributions to the TOA radiance, we included an equal 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. We adopted a bimodal aerosol distribution (Ahmad et al., 2010) with a fine-mode fraction of 0.5. The relative humidity was assumed to be 80%. The lbl contributions to the TOA radiance for this case are shown in Figure 2. The optical depth of the aerosol particles at 532 nm was assumed to be 0.15 in each of layers 14 and 15. Hence the total aerosol optical depth (AOD) was 0.3. Note that the largest lbl contribution to the TOA radiance comes from layer 14, and the next largest from layer 15. The reason is that the aerosol particles in layer 14 effectively shield layer 15 from exposure to radiation due to attenuated contributions from the incident direct beam and multiply scattered diffuse sources.

4.3.3 A Rayleigh scattering atmosphere with an opaque liquid water cloud layer

Instead of adding aerosol particles in layers 14 and 15, we now investigate the impact of adding cloud droplets. Thus, we assumed a liquid water cloud to be present in layers 14 and 15. The optical depth at 532 nm of the cloud droplets was assumed to be 2.25 in each of the two layers. Hence, the total cloud optical depth of layers 14 and 15 was 4.5. The lbl contributions in Figure 3 show that the reflected TOA radiance is dominated by the contribution from the cloud droplets in layer l4. Very little of the incident direct beam radiation and the multiply scattered diffuse radiation penetrates layer 14, explaining why the contribution from layer 15 is negligible.

The results presented here apply equally well as those provided by the traditional lidar approach for cases in which the single-scattering limit prevails, but they also give correct results when multiple scattering must be taken into account. Hence, we may conclude that the QlblC approach is robust, versatile, and computationally efficient (especially compared with Monte Carlo simulations).

4.4 Assessing the importance of multiple scattering

To assess the importance of multiple scattering, we compare in Figure 4 layer-by-layer (lbl) contributions to the reflected radiance at the top of the atmosphere (TOA) obtained for different values of the total AOD for layers 14 and 15. The results shown in the left panel of Figure 4 for a total AOD of 0.015 for layers 14 and 15 may be considered to be representative for those that would be produced by the standard lidar approach based on the single-scattering assumption (ϖΔτ ≪ 1). The middle panel in Figure 4 shows results obtained when the total AOD for layers 14 and 15 is increased by a factor of 10 to 0.15, and the right panel is the same as the left panel in Figure 2, corresponding to a total AOD of 0.3 for layers 14 and 15.

The left panel of Figure 5 shows results obtained for a total cloud optical depth (COD) of 0.045 for layers 14 and 15 (a 100-fold COD reduction compared to the COD in Figure 3). The middle panel of Figure 5 shows results similar to those in the left panel of Figure 3 except that the total COD of layers 14 and 15 was reduced by a factor of 10 to 0.45. For reference the right panel of Figure 5 is the same as Figure 3 for a total COD equal to 4.5 in layers 14 and 15.

These results demonstrate that our new QlblC approach makes it possible to account for multiple-scattering effects in the analysis of lidar attenuated backscatter profiles. In fact, the factor C₁ in Eq. 6 can be used to quantify the increased pathlength and thus the observed tail (delayed return) due to multiple scattering. Figures 4, 5 demonstrate that the QlblbC method can be used to quantify the impact of multiple scattering. For the case of an aerosol layer (Figure 4) as well as a cloud layer (Figure 5), the middle and right panels show results for multiple scattering scenarios, while the left panels show results for situations where the single scattering approximation may be valid.

4.5 Comparisons with Monte Carlo simulations

The Monte Carlo (MC) code described by Hu et al. (2001) was modified/extended as follows to be compatible with the AccuRT simulations reported in Section 4.3:

1. While the MC code described by Hu et al. (2001) applies to polarized radiative transfer simulations, modifications were made to make it suitable for applications to the unpolarized cases reported in Section 4.3.

2. To be consistent with the QlblC results the location (altitude) instead of the time of the last scattering event was recorded.

3. The Monte Carlo method was modified in order to simulate the direct solution of the RTE.

In the backward Monte Carlo method, we can consider the direct solution as the convolution of the particulate single scattering of the laser beam and the solution of the adjoint RTE, which is the radiative transfer for a plane-parallel medium with an incident beam in the reversed direction of the viewing angle. In the semi-analytical Monte Carlo simulation (Hu et al., 2001), we can simulate the similar physical process but with reversed light propagation. We can consider the statistical distributions of the locations of the photons as the adjoint solution of the reversed RTE problem, and estimates of the probability that a photon will be scattered directly to the receiver are obtained from the reversed process of the single scattering of the laser beam. Due to reciprocity, the solutions of the two approaches are equivalent.

In Figure 6 we compare AccuRT multiple scattering results with MC simulations. This simple case is for a homogeneous slab of non-absorbing particles scattering according to the Henyey-Greenstein phase function. Note that the AccuRT and MC results are essentially identical for isotropic scattering (g = 0.0). The results are close for g = 0.8, but the agreement deteriorates with slab optical thickness τ_b and the difference between AccuRT and MC results gets as large as about 2.7% at τ_b = 0.3. Also, both MC and AccuRT results agree with the prediction of Eq. 5 valid in the single-scattering limit when the slab optical thickness is smaller than about τ_b = 0.05.

FIGURE 6. Radiance reflectance for a slab of non-absorbing particles scattering according to the Henyey-Greenstein phase function for two values of the asymmetry factor, g. The AccuRT and MC results are compared with the single-scattering predictions of Eq. 5. The results are shown as a function of slab optical thickness, τ_b, increasing from a small value to τ_b = 0.3.

The backscatter reflectance profile is obtained by plotting results such as those displayed in Figures 1–5 as a function of position (location) in the slab. In the lidar community, these profiles are commonly referred to as “attenuated backscatter” in units per steradian per kilometer. To compare attenuated backscatter profiles generated by the QlblC and MC simulations, we assumed that an “aerosol layer”, consisting of an ensemble of particles that scattered in accordance with the HG phase function, was homogeneously distributed in a slab of geometrical thickness 1,000 m. The total optical thickness of this “aerosol layer” was assumed to be 0.3 (the largest optical thickness considered in Figure 6). Figure 7 shows a comparison of attenuated backscatter profiles obtained by the QlblC and MC simulations for this “aerosol layer”. Note that there is almost perfect agreement for isotropic scattering (g = 0.0), and for g = 0.5, while for g = 0.8 there is an average difference of about 3% which is similar to the difference in the radiance reflectance shown in Figure 6 at slab optical thickness τ_b = 0.3 ($\sim 2.7$%). The average difference of 3% for the attenuated backscatter profile for g = 0.8 in Figure 7 was calculated as the mean value over 10 depths (50 m, 150 m, …, 950 m), with differences varying from 5.76% at 50 m to 0.8% at 950 m.

FIGURE 7. Attenuated backscatter profiles computed by the MC code and the AccuRT QlblC code for the case displayed in Figure 6 for a 1,000 m thick slab of total optical thickness τ_b = 0.3.

5 Discussion

As alluded to in Section 3.1, for a pulsed light beam, time-gating can be used to keep track of the depth (or layer) in the medium from which the observed backscattered radiance originates. Hence, it may be sufficient to solve the steady-state radiative transfer equation (RTE). Also, to address the time dependence, the steady-state RTE can be solved sequentially by first ignoring all layers except the uppermost layer, then solve for the two uppermost layers, and so on.

As briefly discussed in Section 3.2, the RTE involving lidar (finite) beam illumination describes a 3-D problem. The solution of the 3-D RTE for a narrow finite laser beam is quite challenging and computationally demanding. Therefore, most treatments of the lidar problem rely on solving a 1-D RTE for both atmospheric (Hogan, 2008; Hogan and Battaglia, 2008) and aquatic (Mitra and Churnside, 1999) applications. We have adopted the 1-D point of view, but we recognize that a study of the validity of this approach needs further attention, which is beyond the scope of our paper, although we expect the 1-D approach to be a reasonable approximation when the scattering mean free path is less than the lateral extent of the lidar beam.

The factor 1/(1—g) introduced in Section 3.3 can be used to approximate the enhanced pathlength of multiply scattered light. The asymmetry factor g determines the degree of forward scattering. Hence, a larger g-value will increase the mean free path and thus lead to a time delay for multiply scattered light compared to singly scattered light.

The multiple scattering between aerosols and clouds may also influence the standard lidar approach based on the single scattering approximation to derive the contribution of each aerosol or cloud layer as shown in Figure 6. The capability of lidar to measure aerosol properties below thin clouds is indeed an important topic that merits further study to quantify its significance. Such an investigation would be a good target for a follow-up study. Curves similar to those shown in Figure 6 but for the absolute and/or relative differences between single-scattering, QlblC, and MC methods would help better understand the multiple scattering relationship with the asymmetry g and optical depth. One possible reason for the difference between MC and QlblC results for the larger g-value (Figure 7) could be the different treatment of the scattering phase function in the MC and QlblC forward radiative transfer simulations. In the discrete-ordinate method, the forward scattering peak is truncated, while such truncation is not used in the MC method.

6 Summary and conclusion

When a scattering/absorbing medium is illuminated by a beam of electromagnetic (EM) radiation a fraction of the incident EM beam will be reflected by the medium. For a stratified medium, like a multilayered slab of material, that is illuminated by a beam of infinite extent in the direction normal to that of the stratification, another fraction of the incident EM beam will be transmitted. We have described how to quantify specific layer contributions to 1) the reflected or transmitted EM signal emerging from a multilayered, stratified medium that is illuminated by a beam of EM radiation and/or 2) the emitted EM signal emerging from such a medium due to emission by internal sources. To simplify the discussion we assumed that the medium consisted of either 1) one single, multilayered slab with a constant index of refraction or 2) two adjacent multilayered slabs (labeled slab₁ and slab₂) separated by a smooth interface across which the real part of the refractive index changes abruptly from one constant value in slab₁ to a different constant value in slab₂, such as for an atmosphere-water system. This quantification of the emerging radiance originating from any specific layer (or location) in the medium relies on our ability to compute and integrate the source function as explained in Section 2.1 [see Eqs. 1, 2].

In Section 2 we provided the theoretical background, in Section 3 we described how the methodology can be applied to the lidar problem, while in Section 4 we provided numerical examples by applying the method to space-borne lidar probing of a Rayleigh scattering atmosphere that was 1) cloud- and aerosol-free, 2) had an embedded aerosol layer, and 3) had an embedded cloud layer. Finally, we investigated how multiple scattering affects the results for the case in which the air molecules in the two atmospheric layers closest to the surface (layers 14 and 15) were assumed to be mixed with either aerosol particles or cloud droplets. For small values of the total optical depth of aerosols (AOD, see Figure 4) or cloud droplets (COD, see Figure 5) in layers 14 and 15, we obtained results for which the standard lidar equation applies [see Eqs. 19, 21 valid in the single-scattering limit]. For large values of total AOD or COD in layers 14 and 15 (ϖΔτ > 1), we found the reflected TOA radiance to be significantly influenced (and for values of ϖΔτ ≫ 1 dominated) by multiple-scattering effects.

A comparison of QlblC results and Monte Carlo simulations in Section 4.5 shows that the two approaches give almost identical attenuated backscatter profiles for isotropic scattering and small differences (less than 3%) for anisotropic scattering.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

KS provided the theoretical concept and wrote the first draft of the paper. WL and YF did the AccuRT/QlblC computations under KS guidance. YH carried out the Monte Carlo simulations. All co-authors contributed with critical feedback and ideas during paper writing process.

Funding

The authors wish to thank the NASA Remote Sensing Theory program and the NASA ESTO’s IIP program for supporting this research.

Acknowledgments

Constructive suggestions made by two reviewers that helped improve this paper are gratefully acknowledged.

Conflict of interest

Authors CW and JL were employed by Ball Aerospace Corp.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/frsen.2022.1025447/full#supplementary-material

References