Tomography is a retrieval technique inverting a 2D spatial density from a dataset of measured “slices” (τoμoς, tomos means “slice” in Greek), which are the integrals of this density over transect lines (chords). The collection of slices parameterized by angle and offset of the chords (relative a selected central point) constitutes a 2D dataset called tomogram (or “ray sum” in medicine). The original 2D density can be restored back from its tomogram using the inverse Radon transform, introduced in 1917 by Austrian mathematician Johann Radon (Radon, 1917; Radon and Parks, 1986). This transform serves as the mathematical basis for the X-ray tomography (commonly known as the CT scan in medicine, with “CT” referring to “computer tomography”).

While topography is an analysis of essentially active measurements (when a ray sent through the object is captured by a detector on its other side), this term sometimes is used in cloud remote sensing in a wider sense: as any technique of inversion of cloud interior structure from remote optical measurements. Previously developed retrieval algorithms of this kind (e.g., Martin and Hasekamp, 2018; Levis et al., 2020) rely on least-square-fit (LSF) inversions in highly-multi-dimensional space of possible cloud configurations and use extensive 3D RT computations on each iterative step. Levis et al. (2020) call such algorithms “passive tomography” using the word “tomography” in that wider sense (while “slice” measurements or Radon transform may not be necessarily used). Below we will use this definition when addressing LSF techniques as well as our own.

We present a new algorithm for inversion of the internal cloud structure. It is ideologically closer than LSF approaches to “tomography” in the traditional sense, while also relying on passive optical measurements. Such measurements obviously cannot directly provide a proper tomogram (directional optical thickness as a function of angle and offset of the viewing ray), but allow for estimation of it using a nested family of “cloud shapes” corresponding to an array of thresholds in the measured total reflectance. Once such tomogram is obtained, we proceed following the standard tomographic procedure and apply inverse Radon transform to it deriving a 2D field of the extinction coefficients (which requires calibration using independent measurements). So this technique may be called “semi-tomographic”. Note that unlike the LSF techniques, our method does not require iterative 3D RT computations, but relies on a number of empirical assumptions which we validate using simulated data. In fact, we only use simulated data (LES + 3D RT) for testing of our algorithm and for selection of an appropriate “proxy” formula relating the measured reflectance to the corresponding directional COT.

The presented algorithm has been developed for the Research Scanning Polarimeter (RSP) measurements. The RSP (Cairns et al., 1999) is an airborne sensor which high-resolution along-track scanning provides sufficient number of view rays to constrain cloud shapes (Alexandrov et al., 2016a) and also to derive droplet size distribution (DSD) profiles along the cloud side (Alexandrov et al., 2020). The RSP served as a prototype for the satellite Aerosol Polarimetery Sensor (APS) built for the NASA Glory Project (Mishchenko et al., 2006; Mishchenko et al., 2007). Currently the cloud retrieval algorithms developed for the RSP are being adopted for use with the measurements of the spaceborne polarimeters on the upcoming Plankton, Aerosol, Cloud Ocean Ecosystem (PACE) satellite mission (Werdell et al., 2019). During the past decades the RSP has been deployed during numerous NASA field campaigns (see e.g., Alexandrov et al., 2015; Alexandrov et al., 2016a; Alexandrov et al., 2016b; Sinclair et al., 2017; Alexandrov et al., 2018; Sinclair et al., 2019).

This study was inspired by the RSP measurements made during NASA’s Cloud, Aerosol and Monsoon Processes Philippines Experiment (CAMP²Ex) conducted in 2019. This campaign took place off Philippines coast during the monsoon season, which is characterized by abundance of moisture and strong convections producing wide variety of cloud types and shapes which present interesting cases for application of our tomographic analysis. We will report analyses of the actual RSP data from CAMP²Ex in subsequent publications following this paper.

The RSP has high angular resolution of 14 mrad field of view (FOV) with measurements made at 0.8° intervals within ±60° from nadir. The translation of these parmeters into spatial resolutions of the RSP datasets depends on the speed and the altitude of the aircraft (Alexandrov et al., 2016a) and can range from about 200 m for high-altitude NASA ER-2 down to less than 100 m for more conventional NASA airplanes (P3-B, B-200, UC-12, C-130) which fly slower and at lower altitude.

The RSP measures both total and polarized reflectances in nine spectral channels. Its measurement geometry and the view-line aggregation types used for different kinds of retrievals are schematically presented in Figure 1. The polarized reflectances in the rainbow (cloud bow) scattering range (135°–165°) have been routinely used for retrieval of droplet size distributions in liquid-water clouds (Alexandrov et al., 2012a; Alexandrov et al., 2012b).

The COT retrievals from RSP-measured total reflectances at nadir view are made using a modification of the legacy bi-spectral technique (Nakajima and King, 1990). In the modified algorithm no absorbing spectral channels are used. Instead, the droplet effective radius is retrieved from the polarized reflectance (Alexandrov et al., 2012a) and used to derive the COT value from the look-up table (LUT) computed for non-absorbing 863 nm channel. This LUT was built using plain-parallel radiative transfer computations, thus, this method can produce biases in COT values in the presence of 3D radiative effects (such as light escape from broken cloud’s sides or shadowing at low Sun angles).

The standard operational RSP cloud top height (CTH) product is based on a stereo block-correlation algorithm developed by A. Wasilewski. This technique has been generalized to multi-layer cloud scenes by Sinclair et al. (2017).

For illustration and validation of the proposed tomographic technique we use simulated RSP measurements generated by the 3D radiative transfer (RT) model called “Monte Carlo code for the phYSically correct Tracing of photons In Cloudy atmospheres” (MYSTIC Mayer, 2009; Emde et al., 2010) applied to a LES dataset with 100 m horizontal and 40 m vertical resolution. This dataset is an idealized representation of the shallow, maritime convection cloud fields (Ackerman et al., 2004) observed during the Rain in Cumulus over the Ocean project (RICO van Zanten et al., 2011). In the 3D RT computations the actual bin-by-bin LES DSDs were replaced by gamma distributions with the same effective radius as in the microphysical model, while the effective variance was uniformly set to 0.1 throughout the cloud (see Hansen and Travis, 1974; Alexandrov et al., 2018 for definitions). In this simulation the virtual RSP makes measurements at 555 nm wavelength and is assumed to fly in the solar principal plane (corresponding to the y-axis of the LES grid) at 2.4 km altitude. The solar zenith angle is set to 40° and the surface albedo—to 5%. These simulations have been already used in our previous work (Alexandrov et al., 2010; Alexandrov et al., 2012a; Alexandrov et al., 2016a).

We selected the last (No. 3) cloud from the segment used in Alexandrov et al. (2016a). Figure 2A presents the distribution of the droplet number concentration N _c (measured in cm⁻³) within the 2D cross section of this cloud. Figure 2B shows the corresponding field of the effective radius r _eff (measured in μm) of the droplet size distribution, this plot is restricted to the area with N _c ≥ 10 cm⁻³ for consistency with Figure 2A. The effective variance v _eff in the RT model is constant and set to 0.1. The effective radius and variance of cloud droplet DSD are routinely retrieved from the RSP measurements of the polarized reflectance using the parametric algorithm (Alexandrov et al., 2012a).

The extinction coefficient field for the 2D cross section of this cloud is presented in Figure 4A. It is derived from the cloud microphysical parameters using the following formula k e x t = Q e x t N c π ⟨ r 2 ⟩ × 1 0 − 6 , (1) where N _c is the droplet number concentration (in cm⁻³), r is the droplet radius (in μm), and Q _ext is the extinction efficiency (Hansen and Travis, 1974). The factor 10^–6 makes the units of k _ext to be m⁻¹. However, for our purposes we can use a simplified version of Eq. 1 k e x t ≈ 2 π N c r e ff 2 1 − v e ff 1 − 2 v e ff × 1 0 − 6 , (2) where we have assumed that Q _ext ≈ 2 for large cloud droplets and expressed the second moment ⟨r ²⟩ of the droplet size distribution in terms of r _eff and v _eff [assuming that the DSD has gamma-distribution shape (Hansen and Travis, 1974)].

In this study we use the implementation of the Radon transform from the standard Interactive Data Language (IDL v.8.4) library (RADON function) with ramp filter added. The 2D field to be determined is the spatial distribution of the cloud extinction coefficient k _ext (x, y) (measured in m⁻¹), whose integral over a linear transect (chord) through the cloud yields the directional COT τ d ψ , ρ = ∫ − ∞ ∞ k e x t ρ ⁡ cos ⁡ ψ − s ⁡ sin ⁡ ψ , ρ ⁡ sin ⁡ ψ + s ⁡ cos ⁡ ψ d s , (3) where s is the distance along the chord. The geometric setup of the method is schematically shown in Figure 3A. The resulting value of the directional COT is attributed to a point in the tomogram (the Radon space) with the coordinates (ψ, ρ) as depicted in Figure 3B.

The inverse Radon transform deriving the extinction coefficient distribution from the tomogram is called “filtered backprojection.” It consists of two parts. First, the ramp filter is applied to the tomogram τ _d (ψ, ρ) in ρ-direction transforming it into its filtered version τ ̃ d ψ , ρ = F F T − 1 | f | ⋅ τ ̂ d ψ , f , (4) where τ ̂ d ψ , f = F F T ρ τ d ψ , ρ . (5)

Here f is frquency and FFT is the Fast Fourier Transform operator (implemented in IDL by the function FFT). The ramp filter is a high pass filter with V-shaped kernel |f| in frequency domain. It prevents low-frequency blurring of the image. While filtering is an important part of the algorithm, it, unfortunately, has been omitted in the IDL v.8.4 distribution, so we designed our own version.

After filtering the backprojection itself is applied yielding the inverted extinction field: k e x t ′ x , y = ∫ 0 π τ ̃ d ψ , x ⁡ cos ⁡ ψ + y ⁡ sin ⁡ ψ d ψ . (6)

Note that inverse Radon transform is not exact, so the inverted k e x t ′ ( x , y ) may be somehow spatially different in structure from the original k _ext (x, y). The inversion is also defined up to an unknown constant factor, so that k e x t x , y = C k e x t ′ x , y . (7)

The factor C should be determined using a calibration procedure based on an independent data source (we will discuss this in detail in Section 8).

Before considering the conversion of simulated remote sensing measurements into a 2D extinction distribution, we first want to check that our implementation of Radon transform works well on the LES extinction field itself. This means that the consecutive application of direct and inverse Radon transforms should yield the field (almost) identical to the initial one. This test is illustrated in Figure 4. Application of direct transform Eq. 3 to the extinction distribution k _ext (x, y) shown in Figure 4A provides the directional COT tomogram τ _d (ψ, ρ) shown in Figure 4B. It is a function of angle and offset of chords passing through the cloud domain as defined in Figure 3A. Then, (after filtering) inverse transform Eq. 6 is applied to this tomogram resulting in the backprojection field k e x t ′ ( x , y ) . This field is defined up to an arbitrary constant factor and needs to be calibrated. In order to do this we need an independent piece of information that cannot be derived from the tomogram itself. For this purpose we use the (vertical) COT τ(x) derived from the initial LES extinction k _ext (x, y), assuming that it is known from a hypothetical independent observation. Then, we compute COT τ′(x) from the backprojection field k e x t ′ ( x , y ) and scale it by appropriate constant factor (the ratio C between the maxima of τ(x) and τ′(x)) to match the initial LES COT values. Figure 4D shows that this adjustment leads to a perfect match in COT. The scaled backprojection field C ⋅ k e x t ′ ( x , y ) shown in Figure 4C appears to be almost identical to the initial extinction distribution k _ext (x, y) in Figure 4A. Thus, we conclude that the IDL Radon transform code (with filtering added) works very well and can be used for further remote sensing applications.

The RSP’s uniquely high angular resolution coupled with the high frequency of measurements can provide geometric constraints on the cumulus cloud’s 2D cross section and yield cloud’s geometric shape estimates (Alexandrov et al., 2016a). Use of a clear/cloud-separation threshold in the measured total reflectance allows for creation of a 1D cloud mask in each RSP scan distinguishing between the bright cloudy part of it and the darker background. The viewing angles at the edges of such cloud mask correspond to the view lines “tangent” to the cloud shape (given that the position of the aircraft relative to the cloud is accurately known). Figure 5 shows how the collection of all such tangent lines obtained during the overflight effectively “cuts out” a polygon (containing the cloud shape) in the 2D plane below the flight path (red and blue lines in Figure 5 correspond to different edges of 1D cloud masks). Then a disc-inscription technique is used to create a realistically-looking cloud shape within this polygon. In this method for each vortex of the polygon we create the largest disk with the center on the bisector of the corresponding corner, requiring that this disk is entirely contained inside the polygon’s interior (so it is tangent to some of the polygon’s edges). The boundary of the union of all such (usually overlapping) disks is then declared to be the shape of the cloud. In some cases (especially of elongated polygons) this technique does not work well, so the polygon itself provides a better representation of the cloud shape (Alexandrov et al., 2020).

A similar “space carving” methodology has been developed in a 3D case for the Multi-angle Imaging SpectroRadiometer (MISR) measurements (Lee et al., 2018). This technique is also used for the first-step definition of the cloud domain in the LSF tomographic algorithm described by Levis et al. (2020).

The notions of cloud shape or cloud boundary, while intuitively well-understood, are, strictly-speaking, quite ambiguous both in physical and optical sense. Physically, cloud is a spatially distributed collection of water droplets with no hard-defined surface. It can be artificially bounded by a level surface of the droplet number concentration corresponding to a threshold value chosen by the observer. Similarly, “optical surface” of the cloud can be defined using an arbitrarily chosen brightness threshold. It is affected by the solar and viewing geometries. For example, “optical surfaces” are systematically shifted relative “physical” ones towards the side of the cloud directly illuminated by the Sun. The higher is the brightness threshold, the smaller is the cloud-shape cutout derived using the technique described above.

However, while it is difficult to assign a precise physical meaning to the cloud shape corresponding to a single brightness threshold, the collection of such shapes derived for a range of thresholds may carry information about the internal structure of the cloud.

In this study we assume that the reflectance measured by the RSP at a certain viewing angle is quantitatively one-to-one related to the directional COT (dCOT) along the corresponding view ray. This assumption allows us to use the measured reflectance to build dCOT tomogram, which will be then converted (using inverse Radon transform) into a 2D distribution of the cloud extinction coefficient. This assumption is, of course, a simplification in our essentially 3D setting, where multiple scattering within the whole cloud contributes to the reflectance at any particular direction. However, we will show below that this simplified assumption is sufficient for our purposes. When reflectance becomes a proxy for dCOT, it loses its optical properties and gains physical ones. Having this in mind, we introduce “reflectance-proxy” (RP) which coincides with the measured reflectance where it is present, while can be computed for other (virtual) view rays using inter- or extrapolation (in the same way as dCOT can be extended). Some of such view rays (e.g., horizontal ones) are inconsistent with the observation geometry and do not exist in real observation datasets. Also, the extended RP values are not expected to comply with the RT laws, so the RP is a rather abstract dataset that gets its physical meaning only after conversion to dCOT. This conversion will be discussed in Section 7.

Following the above assumptions we create an intermediate reflectance-proxy tomogram by interpolating the reflectances observed at the actual view rays to all possible view rays from a high-resolution grid of angles and offsets (relative to a chosen cloud center). This construction utilizes the collection of cloud shapes corresponding to a range of brightness thresholds. In a realistic analysis, the brightness threshold array consists of discrete values. The corresponding cloud shapes also form a discrete family, which is expected to be nested. Figure 6A presents such a cloud-shape family constructed for our simulated cloud using an array of 5 brightness thresholds. In addition to them, we may also use the maximum observed reflectance value that we associate with the cloud center depicted by white diamond in Figure 6A. Then, for each threshold we assign its value to all points of the corresponding cloud-shape curve. This operation converts the cloud shapes into level curves of an abstract 2D “reflectance-proxy distribution” (RPD), which then is extended to the entire 2D domain by interpolation between the RP values at the level-curve points. Specifically, for each point P in this domain we first select the two level curves (let us identify them as “1” and “2”) between which P is located. Then we determine the shortest distances from P to these two curves: d ₁ and d ₂ respectively. These distances are then used to compute the weights w ₁ = d ₂/(d ₁ + d ₂) and w ₂ = d ₁/(d ₁ + d ₂) with which the respective level-curves’ RP values contribute to the linearly interpolated value assigned to P. Moving-average smoothing takes care of interpolation irregularities. Our RPD is shown in Figure 6B. This RPD can be used for creation of cloud shapes corresponding to threshold values in-between of those from the initial discrete array, however, we use it in a different way.

Given the way the cloud shape corresponding to a certain brightness threshold is constructed, all the actual view rays tangent to it correspond to the same reflectance value equal to this of the threshold. We extend this property to all possible (not only actual) view rays tangent to this cloud shape. This, of course, is an abstraction not necessarily consistent with the nature of realistic light scattering within the cloud. This extension allows us to assign a RP value to any view ray (which from now on we will also call “chord”) crossing the RPD domain. The RP value assigned to a chord is this of the RPD’s level curve to which this chord is tangent. For the functional shape of the RPD in our case (which monotonically increases towards the cloud center) it is easy to see that this value is simply the maximum of the RPD along the chord.

In more complicated cases, such as e.g., two partially merged clouds, a chord may be tangent to two level curves with different values (each surrounding a different cloud center). In this case the along-chord RPD would have two local maxima, and the value of the larger of them should be assigned to the chord. This follows from the general structure of RPD with higher-value level curves being inside the interiors of the lower-value ones (since clouds are brighter at their center(s) than at their edges). So a chord coming from an observation with a lower threshold b cannot be tangent to a contour curve corresponding to a higher threshold a > b, since to do so it would have to cross another b-value level curve encompassing the a-value one (instead of being tangent to it). On the other hand, chords corresponding to a higher threshold a can cross level curves of any lower value since they are invisible for the a-based cloud masking.

While the RPD can be defined on a 2D grid with very fine resolution (we use 1 m, which certainly is an excess), the chords’ RP values can be also combined into a function of angles and offsets on a very fine grid (we use 1° resolution in angle and 1 m—in offset). This function is the RP tomogram R _tom (ψ, ρ). Figure 7A presents the RP tomogram computed for our cloud. Alongside determining the chords’ RPD maxima for the RP tomogram, we also collected the lengths of these chords when bounded by the cloud shape corresponding to the lowest brightness threshold (here 0.0015). These lengths also form a tomogram L _tom (ψ, ρ) defined on the same grid as R _tom (ψ, ρ); it is shown in Figure 7B.

On the next step of our analysis we need to find out how to relate the dCOT tomogram to the RP- and the chord-length-tomograms from Figure 7. Note that the dCOT tomogram can be determined up to an arbitrary constant factor, since the inverse Radon transform of it is also defined up to such a factor. The actual dCOT tomogram which we want to estimate from the RSP measurements is shown in Figure 4B. Note that using L _tom alone as the dCOT tomogram would result in a constant extinction value within the implied cloud domain. Indeed, having a constant extinction coefficient makes the dCOTs to be proportional to the corresponding chord lengths. This suggests that inclusion of L _tom into analysis together with R _tom would provide us with a better representation of the spatial configuration of the cloud.

Finding a good empirical proxy formula relating R _tom and L _tom to the dCOT tomogram τ _tom is a trial and error process, in which the LES data plays an active role (since we know what the result should be). After several trials we have stopped at the formula relating the reflectance R to the dCOT τ in the case of single scattering and normal incidence of a wide beam on a plane-parallel medium: R = 1 2 1 − e − 2 τ b , (8) where b is a fitting parameter playing the role of the backscattering coefficient. It is assumed to be constant throughout the medium and not related to other physical parameters of the cloud; we use b = 0.1. Equation 8 translates into the following proxy formula for the dCOT tomogram τ t o m ψ , ρ = − ln 1 − 2 b R t o m ψ , ρ L t o m ψ , ρ 2 max L t o m . (9)

Here the denominator 2 max L _tom is actually irrelevant (since the inverse Radon transform of τ _tom is defined up to a constant factor) and is included only to make τ _tom dimensionless as dCOT should be. The dCOT tomogram for our cloud (with an arbitrary value scale) is presented in Figure 8A.

As we mentioned above, the extinction distribution derived from the dCOT tomogram using inverse Radon transform is defined up to an arbitrary constant factor and has to be calibrated. The calibration factor can be determined by iteratively computing the 3D-RT-simulated reflectances for scaled extinction fields until they match the RSP measurements. This way is computationally expensive and also involves uncertainties associated with the 2D nature of our retrievals that do not provide the full extent of the 3D cloud structure. A better alternative to this is to rely on additional measurements such as the COT derived from the RSP’s nadir reflectances using a LUT, or cloud-top extinction coefficient retrieved from the measurements made by the High Spectral Resolution Lidar (HSRL) which is often deployed on the same airborne platform as the RSP during field campaigns.

The algorithm for COT retrievals from the RSP’s nadir total reflectances uses a LUT computed assuming a plane-parallel geometry. Thus, while working fairly well for stratiform clouds, the retrievals for small popcorn Cu can be significantly affected by 3D effects such as leaking of scattered light through the cloud’s sides and shadowing (see e.g., Marshak et al., 1999; Várnai and Marshak, 2001; Várnai and Marshak, 2002). Both effects reduce the nadir reflectance (by about a factor of two in our case) compared to that from a stratiform cloud with the same COT. This causes significant underestimation of COT retrievals as long as the measurements are interpreted using plane-parallel LUT. Figure 9A illustrates this problem. Here the initial COT derived from the LES microphysics is shown in blue and has the maximum of about 35. The red curve represents the COT computed using the extinction field resulted from our tomographic procedure (Figure 8A) and rescaled to have the same maximum as the LES-derived COT. It is appeared to be shifted left by about 50 m towards the bright side of the cloud relative the blue curve. The solid green curve depicts the COT derived using the RSP’s nadir-view algorithm applied to the virtual RSP measurements. While being shifted to the left even further than the red curve, the value of the optically derived COT underestimates the microphysical one by the staggering factor of 4, and is certainly not suitable for calibration of the tomographic retrievals.

An alternative calibration method can be based on the extinction coefficient values at cloud top presumably known from HSRL retrievals. Figure 9B illustrates such scheme. Here we assume that virtual HSRL (which we do not currently simulate) provides accurate extinction coefficients (coinciding with the LES values) near cloud top (here at 1,400 m altitude). Then we find the calibration factor by scaling the tomographic extinction maximum at this altitude to match the maximum of the HSRL curve (shown in blue). The resulting tomographic extinction is shown in red. Application of this calibration factor to the tomographic COT results in only 16% overestimation (at maximum) of the LES COT from Figure 9A (and the same overestimation of the extinction coefficient).

We continue to explore the potential of the RSP’s own measurements of the polarized reflectance to provide the calibration constant. They are dominated by single scattering and, therefore, are less affected by 3D effects. However, for now, we want to separate calibration issues from assessing the accuracy of the characterization of the 2D cloud structure. For this purpose we adopt the calibration made by matching the tomographic COT maximum to the LES COT maximum. The tomographic COT then corresponds to the red COT-curve in Figure 9A. The extinction field in Figure 8B is calibrated accordingly. The same calibration is used in the comparisons presented below.

Figure 10 presents the comparisons between the tomographic retrievals of the extinction coefficient from Figure 8B (right) with the initial LES extinction field. Plots in the left column of Figure 10 show that the retrieved field appears to be consistently shifted relative the original one by about 50 m to the left (towards the bright side of the cloud). This is a side effect of using brightness thresholds for derivation of cloud shapes: brighter parts of the cloud look optically and, therefore, physically thicker. Nevertheless, both fields show quite a structural similarity. To better demonstrate this in right-column plots we shifted the retrieved field to the right by 50 m, thus, improving the comparisons. Figure 10A shows color contour plot of the LES extinction coefficient (from Figure 4A) with the two level curves of the retrieved distribution being over-plotted. White curve corresponds to k _ext = 0.02 m⁻¹, black—to k _ext = 0.07 m⁻¹ (these values are represented by blue and red respectively in the color plot). We see quite good spatial agreement between the two fields in Figure 10A, left: within about 10 m both vertically and horizontally. The retrieved field also replicates the characteristic tilted-column shape of the initial cloud. For more detailed comparisons we plotted in Figure 10B the horizontal profiles of both fields along the transect at 1,100 m altitude. Similarly, Figure 10C compares vertical profiles of extinction coefficient along the transect at the cloud center.

Figure 11 presents quantitative comparison between the retrieved extinction coefficient with that from the initial LES dataset. First, the retrieved k _ext field (with 1 m × 1 m spatial resolution) is interpolated to the LES data points (located on 100 m × 40 m grid). Then only the data points with positive values in both datasets are taken into comparison (to avoid the influence of zero-value points outside the cloud). The direct (as in Figure 10A, left) scatter-plot comparison of k _ext is presented in Figure 11A. Figure 11B shows the same type of comparison but with the retrieved field shifted to the right (away from the brighter cloud side) by 50 m (as in Figure 10B, right). The two datasets appear to be unbiased (having almost zero-mean difference) that is not surprising given that the LES COT was used to calibrate the RSP-derived extinction field. So the measure of the retrieval accuracy in this case is the standard deviation σ _k of the difference, which is 0.015 m⁻¹ (20.5% of the LES dataset’s maximum) for the un-shifted field and improves with the 50 m shift to 0.011 m⁻¹ (15.1% of the LES dataset’s maximum). The correlation between the two datasets also improves with the shift: from 73 to 84%. The LES extinction values used in this comparison have the mean of 0.023 m⁻¹, the median of 0.017 m⁻¹, and the maximum of 0.073 m⁻¹. The width of the scatter plot in Figure 11B does not seem to depend on the k _ext magnitude, so the best overall accuracy assessment is that 65% of the retrieval datapoints have extinction values within σ _k = ±0.01 m⁻¹ from their LES counterparts, and 96%—within 2σ _k = ±0.022 m⁻¹ (the corridor shown by dotted lines in Figure 11). For the un-shifted dataset 97% of points lie within its 2σ _k -corridor with the half-width of 0.030 m⁻¹.

Alexandrov et al. (2020) described a method for derivation of DSD vertical profiles from RSP-measured polarized reflectances. This method is a further development of our operational algorithms (Alexandrov et al., 2012a; Alexandrov et al., 2012b) based on the polarized rainbow analyses. The only new feature introduced by Alexandrov et al. (2020) was the aggregation of the RSP measurements to a range of points at the cloud’s side instead of the stereo-derived cloud top. The “cloud’s side” here is a part of the cloud-shape curve corresponding to some brightness threshold. Only points from the bright side of the cloud are considered. Note that in the case of a cloud tilted away from the Sun (as in Alexandrov et al., 2020) the stereo CTH on the bright side is usually close to a cloud-shape curve. Thus, the standard RSP DSD retrievals can be used for derivation of the profile without performing an alternative aggregation. This, however, is not the case for the cloud considered in this study which is slightly tilted toward the Sun.

The RSP-measured polarized reflectances are dominated by single scattering and are representative of the cloud layer of unit COT or about 50 m into the cloud (Alexandrov et al., 2018). Thus, the cloud shape chosen for derivation of a representative DSD profile should be somehow deeper into the cloud then the very first detectible curve corresponding to the smallest brightness threshold. On the other hand, the vertical extent of cloud shape rapidly decreases with the brightness threshold (as seen in Figure 6), thus, shortening the range of the corresponding DSD profile. In order to balance the range and representativeness of the profile we chose the cloud shape corresponding to two different thresholds one for the bright and the other for the shadowy sides of the cloud: 0.030 and 0.005 respectively (Figure 5B).

Figure 12A schematically illustrates aggregation of RSP measurements to a single point on the cloud surface from Figure 5B. The green lines represent the view rays (from different aircraft positions) corresponding to the view angles within the rainbow scattering range (135°–165°). Figure 12B shows the part of the cloud shape selected for derivation of the effective radius (Figure 12C) and variance (Figure 12D) profiles using the parametric fit method (Alexandrov et al., 2012a). Note that MYSTIC 3D RT code does not use bin-resolved DSDs (instead assuming gamma-distribution functional shapes with LES-derived r _eff and v _eff = 0.1). Hence, the non-parametric RFT retrieval technique (Alexandrov et al., 2012b) does not provide any additional information in this case.

Cloud-side profiles of DSD parameters appear to provide a good estimate of their inside-cloud counterparts both in simulated and real datasets. The 2D cross section of the r _eff distribution in our LES output presented in Figure 2B shows much stronger variation in the vertical dimension than in the horizontal one. Also, our statistical analysis of the data collected during CAMP²Ex demonstrated very good agreement (on average) between the DSD parameters derived from the RSP measurements at cloud tops and those measured in situ inside clouds when both datasets were plotted vs. their respective altitudes (these results will be published elsewhere). The LES model used by Levis et al. (2020) for demonstration of their LSF tomographic algorithm assumed the same r _eff profile and the same v _eff value of 0.1 throughout the cloud. While any LSF algorithm generally works as a “black box” not revealing dominant sources of specific retrievals, the tests performed by Martin and Hasekamp (2018) on simulated layer clouds suggest that the vertical profile information likely comes from the cloud side. They show that when the horizontal coverage of the cloud layer achieves 100% of the measurement domain (so its side is no longer seen) the retrievals loose any sensitivity to vertical profile.

In our case r _eff in Figure 12C and v _eff in Figure 12D both show little variability with altitude within the profile range. While the LES r _eff (blue curve) has a modest increase with height from 14 to 16 μm, this trend is not captured by the RSP retrievals (red curve). Thus, within our retrieval accuracy we can simply assume the constant values r _eff = 15.7 μm and v _eff = 0.14 throughout the cloud.

The effective radius and variance profiles derived from the polarized reflectances at the cloud side are assumed to be good representations of the corresponding profiles inside the cloud (as we discussed in the previous section). Then the vertical profiles r _eff(h) and v _eff(h) can be used to obtain the 2D cross section of the droplet number concentration N _c (x, h) from that of the extinction coefficient k _ext (x, h) using Eq. 2: N c x , h = k e x t x , h 2 π r e ff 2 h 1 − v e ff h 1 − 2 v e ff h × 1 0 6 . (10)

Here x and h are the horizontal and vertical coordinates respectively; N _c is measured in cm⁻³, k _ext—in m⁻¹.

In our case droplet size profiles are considered to be altitude-independent (r _eff = 15.7 μm, v _eff = 0.1), so N _c (x, h) is proportional to k _ext (x, h) and its plot in Figure 13A is a rescaled version of that in Figure 8B. Figure 13B presents comparison of the retrieved N _c with the initial LES field (shown in color). Over-plotted white curve is the level curve of the retrieved N _c (Figure 13A) corresponding to 20 cm⁻³, black curve corresponds to the 60 cm⁻³ level (these level values are represented by blue and red respectively in the color plot). As in Figure 10A (right) the retrieval curves here are shifted to the right by 50 m. Comparisons between retrieved (red) and initial LES (blue) droplet number concentrations for the horizontal transect at 1,100 m altitude and the vertical transect at the cloud center are presented in Figures 13C,D, respectively. Retrieved N _c field in both bottom plots is also shifted by 50 m to the right.

Figure 14 presents quantitative comparison of the retrieved N _c values with those in the LES dataset. As in Figure 11 the retrievals are interpolated to the LES data points. Only the points with N _c > 1 cm⁻³ and positive k _ext (i.e., positive r _eff) in both datasets were selected for the comparison. The LES N _c values have the mean of 28.75 cm⁻³, median of 28.36 cm⁻³, and maximum of 71.74 cm⁻³. Figure 14A shows the original scatter plot, while in the right panel the retrieved field is shifted to the right by 50 m (as in Figure 13B). The droplet number concentration comparisons are not as good as those of the extinction coefficient due to limited accuracy of our retrieved droplet size profiles. However, they demonstrate the possibility of quantitative estimation of 2D N _c field. There are very small negative biases in both plots (less than 0.5 cm⁻³), while the standard deviation σ _n provides more definitive accuracy estimate. As expected, the accuracy improves with the shift from 17.6 to 12.8 cm⁻³, while the correlation increases from 65 to 81%. In Figure 14A 96.5% of points lie within the 2σ _n -corridor (depicted by two dotted lines). The shift does not change this number much (97.7%) while the corridor is visibly narrower in Figure 14B. Thus, we can estimate our accuracy of N _c retrievals as ±13 cm⁻³, regardless of the N _c value (within the range of our model), which in this case constitutes 17% of the maximal N _c and 43% of its mean or median.

The ability to derive 2D droplet number concentration field opens possibility to calibrate tomographic retrievals using in situ measurements inside the cloud.

We presented a new remote sensing tomographic technique allowing for retrieval of cloud internal structure from external measurements made by the Research Scanning Polarimeter. While tomography (in narrow sense) is an active technique incompatible with the geometry of atmospheric remote sensing, we developed a “semi-tomographic” approach in which the tomogram of the cloud is estimated from passive measurements instead of being measured directly. This tomogram is then converted into spatial distribution of extinction coefficient using inverse Radon transform (filtered backprojection), which is the standard procedure in the actual tomography (used e.g., in medical CT scans). This procedure is computationally inexpensive compared to approaches relying on multi-dimensional least-square fitting since it does not require iterative 3D RT simulations.

We illustrated and validated the proposed technique using simulated RSP observations of a LES-generated Cu cloud. The radiation field was computed using MYSTIC 3D RT code. On the first step of our algorithm we use the RSP’s view rays grazing the cloud “surface” to create a nested family of cloud shapes corresponding to a range of radiometric thresholds. This family is used for estimation of the “reflectance tomogram” which is then converted into COT tomogram using the proxy relation Eq. 9. While this empirical relation is the result of an educated guess, its validity has been proven using the simulated data.

In the second step Radon transform is applied to the COT tomogram resulting in 2D distribution of the extinction coefficient defined up to an unknown constant factor. To determine this factor a calibration procedure based on additional data should be used. While the RSP’s own measurements would be preferable for the calibration, this appeared to be problematic. Our first choice of the calibration dataset was the RSP-derived nadir COT (given that the LES-derived COT was successfully used for self-test of Radon transform on the LES data). However, unfortunately, the COT derived from the RSP’s nadir measurements assuming plane-parallel geometry appeared to have a significant low bias due to 3D effects (e.g., leaking of light through the cloud’s sides) and could not be used for calibration. We consider the RSP’s polarized reflectance as a potential calibration source since it is dominated by single scattering and is less sensitive to 3D effects. Alternatively, the values of extinction coefficient at cloud top derived from correlative lidar measurements can be also used (this is a relatively new data product developed by the HSRL team). In our simulations we found that this type of calibration leads to a rather modest 16% overestimation of the overall 2D extinction values and the COT. While more work is needed to resolve the calibration issues, in this study we wanted to estimate the accuracy of the method apart from the calibration uncertainty. To do this we used the “ideal” calibration based on the LES-derived COT (assuming that it is known from a hypothetical independent measurement).

The reflectance-based cloud shapes at the first step of our algorithm are naturally shifted towards the bright side of the actual cloud and this bias propagates to the retrieved 2D extinction field. Thus, a point-by-point comparison of this field with the initial LES-derived distribution showed notable improvements after 50 m shifting away from the bright side of the cloud. We understand that at different cloud and illumination conditions the optimal shift may be also different. However, in real remote sensing the retrievals of cloud’s physical structure are more important than that of its exact location, so an error in the latter can be tolerated (see Supplementary Material for other potential sources of the retrieval uncertainties). The retrievals of the extinction coefficient appeared to be unbiased in value, so we measured the retrieval uncertainty using the standard deviation σ _k of the difference with the LES values. We found that for the shifted dataset σ _k = 0.01 m⁻¹ (13.7% of the LES dataset’s maximum) and 98% of the retrieved values of k _ext lie within 2σ _k = 0.02 m⁻¹ from their LES counterparts.

In the next step we converted the extinction retrievals into these of cloud droplet number concentration N _c . In order to do this we first estimated the droplet size distribution parameters (effective radius and variance) within the cloud. We assume (based by the LES data) that the DSD-parameters profiles do not significantly vary in the horizontal direction (which may not be the case for other real or simulated clouds). This allows us to use throughout the cloud the profiles of r _eff and v _eff derived from the RSP’s polarized reflectance measurements made along the cloud side. These retrievals were made using the rainbow (cloud bow) fitting technique (Alexandrov et al., 2012a; Alexandrov et al., 2018). In our example both DSD parameters appeared to be nearly constant within the retrieval accuracy and vertical range making the conversion of the 2D k _ext field into the spatial distribution of N _c particularly straightforward (just a rescaling). As in the extinction case the comparison of the N _c retrievals with the corresponding LES values improves with the 50 m shift achieving the accuracy (the standard deviation of the difference) σ _n = 12.3 cm⁻³ (17% of the maximal N _c ) while 95.3% of the points in the scatter plot lie within the 2σ _n -corridor. The correlation between the two datasets is 80%.

In the upcoming continuation of this series we will present examples of tomographic analyzes of real clouds observed during NASA’s Cloud, Aerosol and Monsoon Processes Philippines Experiment (CAMP²Ex) conducted in 2019. During this campaign the RSP was deployed on-board of NASA’s P3-B aircraft together with a number of other remote-sensing and in situ instruments, measurements of which we can use to both improve and validate our retrievals. In particular, we plan to use the extinction coefficients at cloud top derived from the measurements made by the High Spectral Resolution Lidar (HSRL2) to calibrate our tomographic extinction field, while the HSRL cloud-top height product will be used to validate our retrievals of the cloud geometry. We will also look at the measurements from the Airborne Precipitation Radar (APR-3) to validate our retrievals of cloud shape and position.

A detailed analysis of the requirements to the sensor’s angular range and resolution will be also presented in our future publications (with different RSP platforms in mind, as well as other airborne and satellite instruments). Some preliminary estimates of the tomographic retrieval accuracy are presented in Supplementary Material accompanying this paper.