Multi-View Polarimetric Scattering Cloud Tomography and Retrieval of Droplet Size

Abstract

Tomography aims to recover a three-dimensional (3D) density map of a medium or an object. In medical imaging, it is extensively used for diagnostics via X-ray computed tomography (CT). We define and derive a tomography of cloud droplet distributions via passive remote sensing. We use multi-view polarimetric images to fit a 3D polarized radiative transfer (RT) forward model. Our motivation is 3D volumetric probing of vertically-developed convectively-driven clouds that are ill-served by current methods in operational passive remote sensing. Current techniques are based on strictly 1D RT modeling and applied to a single cloudy pixel, where cloud geometry defaults to that of a plane-parallel slab. Incident unpolarized sunlight, once scattered by cloud-droplets, changes its polarization state according to droplet size. Therefore, polarimetric measurements in the rainbow and glory angular regions can be used to infer the droplet size distribution. This work defines and derives a framework for a full 3D tomography of cloud droplets for both their mass concentration in space and their distribution across a range of sizes. This 3D retrieval of key microphysical properties is made tractable by our novel approach that involves a restructuring and differentiation of an open-source polarized 3D RT code to accommodate a special two-step optimization technique. Physically-realistic synthetic clouds are used to demonstrate the methodology with rigorous uncertainty quantification.

1. Motivation, Context & Outline

Clouds play a significant role at local and global scales, affecting weather, the water cycle, solar power generation, and impacting the Earth’s energy balance [1]. Moreover, uncertainties in global climate models are significantly affected by our limited understanding, and therefore modeling, of cloud dynamics and microphysics [2]. Thus, understanding, modeling, and predicting cloud properties is a key issue with worldwide socio-economic implications that is in the center of many research studies [3]. Much of the current understanding relies on routine remote sensing of cloud properties such as by the MODerate resolution Imaging Spectrometer (MODIS) [4]. In practice, global-scale retrievals have so far been based on an individual pixel basis, using an approximation that clouds are plane-parallel slabs. This approximation uses a 1D radiative transfer (RT) model, which leads to biases in many retrievals [5] while other retrievals simply fail [6]. Convective clouds are therefore a blind spot due to their inherently 3D nature.

In its 2018 Decadal Strategy for Earth Observation from Space [7], the National Academies of Sciences, Engineering, and Medicine have indeed identified “Clouds, Convection, and Precipitation” as one of its five top-priority Targeted Variables for NASA’s next generation of satellite missions. To bridge this gap, new technology is needed to study clouds as 3D volumetric objects, on a global scale. The CloudCT [8] space mission, by the European Research Council (ERC), is specifically destined to provide data and products for this goal. It will involve 10 nano-satellites orbiting in formation, thus acquiring simultaneously unique multi-view measurements of such vertically-developed 3D clouds (Figure 1).

Moreover, common retrieval of cloud droplet characteristics use two optical bands simultaneously [9]: a visible band, where reflected radiance increases with cloud optical thickness, and a shortwave infra-red (SWIR) band, where absorption by condensed water depends on cloud droplet size. To sense droplet size in 3D by CloudCT or other future missions, sensors will need to have either SWIR or polarization capability.

1.1. Why Polarized Light?

There is an additional caveat in common retrievals, which rely on SWIR absorption [9]. In addition to absorption, light undergoes multiple scattering in clouds. Multiple scattering diminishes sensitivity to droplet microphysics. High sensitivity to microphysics is embedded in single-scattering events. It is thus beneficial to pick-up single-scatter signals, out of the strong multiply-scattered background radiance. Polarization signals of scattered light are dominated by single-scattering events, and are thus highly sensitive to the type and size specifications of scatters. Polarization therefore provides a significant signal for retrieval of droplet size distributions. In contrast, the intensity signal, which undergoes multiple-scattering events before reaching the sensor, is insensitive to the droplet size and provides complementary information about optical densities within the cloud.

In recent years, there is growing interest in polarimetric imagers for remote sensing of clouds and aerosols [10,11,12,13,14,15,16]. In turn, increased interest in polarimetric sensing capabilities has led to the development of 1D and 3D polarized (or “vector”) RT codes [17,18] with an aim of improving retrieval algorithms. Motivated by the CloudCT mission formulation—only the first of many to come in innovative passive cloud remote sensing—we develop herein a novel framework for 3D remote sensing of cloud properties using multi-view polarimetric measurements.

1.2. Why Passive Tomography?

From its etymology, the word “tomography” means a slice-by-slice recovery of an object’s 3D internal structure using 2D projections of cumulative density. In the computer age, this task is termed Computed Tomography (CT) [19]. Common medical CT approaches are transmission-based X-ray CT or single-photon emission computed tomography (SPECT). There, 2D projections represent straight line-of-sight (LOS) integrals of the local X-ray opacity or nuclear marker density, respectively. In both imaging modalities, the inverse problem of recovering the medium content is linear [20].

Medical CT modalities generally use active radiation. Active methods are also used for atmospheric sensing or scatterers by radar and lidar. There, a transmitter and receiver are generally collocated and signals are based on backscattering and time-resolved two-way transmission. Probing is solved per LOS using methods which are computationally relatively simple. However, the technology is expensive, horizontal sampling is generally very limited, and irradiance decays fast from the transmitter. Passive sensing is less expensive, uses minimal power, and can image wide swaths of Earth. Thus global coverage mandates passive imaging from space. Consequently, this paper focuses on the derivation of 3D passive tomography of scatterer fields. In passive imaging of scatterers, the light source irradiating the atmosphere is the sun: uncontrolled, steady and mono-directional. Passive remote sensing does not benefit from pulsed sources for echo-location. It should rely on multi-angular data. In general, however, retrieving atmospheric scatterer fields in 3D requires a full forward model of scattering in 3D. The model generalizes both the direct transmission model of linear CT and the diffusion limit. It can handle atmospheric regimes where neither a direct transmission model nor the diffusion limit are valid. 1.3. Context

In one sense or another, without tomography per se, three-dimensional cloud analysis using remote sensing data has been a topic of research for over two decades. There are well-known biases in optical thickness retrievals using the independent pixel approximation (IPA). To remove these biases, Marshak, et al. [21] invokes an inverse nonlocal IPA (NIPA). Marshak et al. [21] further demonstrates this approach on a single realization of a stochastic model for marine stratocumulus (Sc). NIPA is based on the RT Green’s function for horizontal (cross-pixel) transport of solar radiation between irradiation and escape to the above-head sensor, which can be computed numerically and conveniently parameterized. Using many instances of the same stochastic cloud model, neural networks (NNs) were successfully trained to retrieve domain-average cloud properties from nadir radiance fields [22,23,24]. Zinner et al. [25] generalized NIPA to multi-spectral retrievals of both optical thickness and column-averaged effective radius. Iwabuchi and Hayasaka [26] creatively combined NIPA-like physics-based and NN-like statistical regression methods to retrieve spatially varying properties of Sc clouds without suffering from IPA-induced biases.

Multi-angle Imaging Spectro-Radiometer (MISR)/Terra [27] was launched in 1999 and soon inspired research into cloud or cloud field reconstruction based on its unique multi-angle imaging capability. Marchand and Ackerman [28] used the original scalar version of the Spherical Harmonics Discrete Ordinates Method (SHDOM) iteratively to fit MISR observations of marine Sc by modifying the column optical thickness under each pixel. Other early pioneers were Seiz and Davies [29] who reconstructed the outer shell of a 2D transect through a large convective cloud by applying labor-intensive stereographic methods to MISR’s multi-angle imagery. Cornet and Davies [30] then estimated the mean optical thickness of that same cloud by performing repeated Monte Carlo 3D RT simulations, holding the outer shape constant, to fit the MISR radiometry. Evans et al. [31] trained a neural network to improve on the IPA retrieval of mean optical depth and its standard deviation using one or more MISR angles. To that effect, these last authors used SHDOM to compute multi-angle radiance fields from very many cloud fields generated with a Large Eddy Simulation (LES). This started a new trend in uncertainty quantification (UQ) in cloud remote sensing, using LES for ground truth that has by now become the standard, and used here as well.

There is a renewed interest in determining a cloud’s outer surface in its own right. On the one hand, very high resolution imagery from ground [32] or space [33] has been processed by stereophotogrammetric methods to generate dense “point clouds” representing many features seen from two cameras. On the other hand, the multi-view technique of “space carving”—used here for initializing the iterative tomography—defines the 3D region that contains a certain cloud by intersecting threshold-based cloud masks from multiple directions. Space carving methods have been applied to determine conservatively the outer shape of clouds observed with the airborne hyperangular Research Spectro-Polarimeter (RSP) [34], MISR [35] (same cloud as Reference [29]), and other sensors (e.g., References [36,37]). Note that methods based on cloud-masking, feature detection-and-tracking and/or stereoscopic location do not require radiometrically calibrated sensors. The loosened radiomeric use seems sub-optimal, since most sensors are in fact painstakingly calibrated to stringent (NIST-traceable) standards. Zinner et al. [38] have, to some extent, reconciled the notion of cloud surface sensing and RT-based cloud property estimation. Their approach scans the illuminated sides of vertically-developed convective clouds using a multispectral sensor mounted on an airborne platform. Then, a standard bi-spectral technique infers the vertical profile of cloud particle size along the side of the cloud. Alexandrov et al. [39] used LES clouds and 3D vector RT to demonstrate cloud-side microphysical profiling with RSP’s scanning polarimetric capability.

Remote sensing of clouds in complex 3D scenarios has already benefited from the advent of convolutional neural networks (CNNs), which excel in mapping one or more (e.g., multi-spectral) images, to inherent properties of the objects in the scene. The CNN is trained on a large database where the multi-spectral radiances are simulated with 3D RT and the pixel-scale cloud properties (typically, optical thickness and column-mean effective droplet radius) are known from the LES and assumed microphysics. Deep learning has thus been applied to both nadir satellite imagery [40] and to an array of ground-based sky cameras [41].

To close this non-exhaustive survey of cloud or cloud scene reconstruction techniques, we mention data fusion, which generally involves a combination of active and passive sensors. Liou et al. [42] blazed this path toward the inference of 3D cloud structure by creatively using a combination of ground-based mm-wave radar and multi-spectral imagery from a NOAA satellite overflying the same region. The radar delivers only a vertical 2D “curtain” of radar reflectivities from an extensive cirrus layer, while the multi-spectral imager provides raw information about the horizontal structure of the same cirrus layer. Barker et al. [43] used A-train satellite data (CloudSat, CALISPO, MODIS) to construct a plausible 3D cloud scene by populating off-nadir columns (where actively-probed profiles are lacking) by looking for a sub-satellite (hence actively profiled) pixel that is close in multi-spectral space. This “donor–receiver” algorithm is then iterated until a broadband 3D RT estimation of up-welling flux at the top of the atmosphere (TOA) agrees with the CERES measurement of it, thus satisfying a form of radiative closure. Fielding et al. [44] formalize cloud scene reconstruction with fused ground-based passive/active observations using an ensemble Kalman filter technique. Ensemble Kalman filtering is common in data assimilation for numerical weather prediction. Fielding et al. [44] assimilate observations into a suite of forward models based on 1D or 3D RT, however, rather than into the primitive equations of atmospheric dynamics.

1.4. Outline

Linear CT models (analogous to medical Xray CT and SPECT) were used to study gas emission and absorption in 3D plumes in the vicinity of pollution sources [45,46] or volcanoes [47,48]. There, Rayleigh-scattered sunlight was transmitted through the gas to a spectrometer on a platform flying around the plume. Following the vision of Werner et al. [49], Huang et al. [50,51] used scanning microwave radiometers to reconstruct 2D slices of particle density in clouds based on its impact on local emissivity. Linear CT was also adapted by Garay et al. [52] to characterize a smoke plume over water emanating from a coastal wild fire. There, the signal is sunlight scattered to space and detected by the Multi-angle Imaging Spectro-Radiometer (MISR) sensor [27] at nine viewing angles. The analysis in Reference [52] yields the direct transmission through the plume per LOS, from which linear CT analysis yields the plume density without using solar radiometers under the plume. Aides et al. [53] formulated CT based on single-scattered light. Their forward model is based on sets of broken-ray paths, where light changes direction once from the sun to a sensor. All the above atmospheric tomography methods assume the medium to be optically thin enough for direct and once-scattered radiation to dominate the measured radiance.

Biomedical imaging also involves CT modalities which are not based on linear models. A notable case in scattering-based tomography of X-ray CT is in Reference [54]. It exploits scattered radiation to obtain reconstruction with significantly lower dose and simpler acquisition hardware relative to linear CT systems, while enhancing estimation of chemical decomposition per volume element. Another example is Optical Diffusion Tomography (ODT) [55,56,57], which uses non-ionizing near-infrared light. It is worth noting the work by Che et al. [58] which departs from physics-based approaches into the realm of machine learning.

In contrast to ODT and the linear models described, which rely on either very high or very low orders of scattering, we use a steady-state polarized 3D RT forward model which generalizes far beyond these assumptions. This approach can seamlessly extract information about droplet sizes and densities from all orders of scattering.

Specifically, we adopt the vector Spherical Harmonics Discrete Ordinates Method (vSHDOM) [59,60], a popular computational 3D polarized RT model that, importantly here, is open source. We then formulate an inverse 3D RT problem for cloud tomography utilizing multi-view multi-spectral polarimetric images. In contrast to linear CT, the image formation model is nonlinear in the cloud’s microphysical and density variables. Our approach seeks an optimal fit of observed data to a physical model based on droplet microphysical and density parameters. Our work here generalizes our prior inversion approaches, which have been based on intensity-only observations [61,62,63,64,65]. Our generalization takes advantage of polarimetric measurements, thereby providing the benefits described in Section 1.1.

In Section 2, we cover basic cloud droplet optics using Mie scattering theory and the fundamentals of polarized 3D RT. The latter yields radiance which has a clear decomposition into single- and multiply-scattered light. This decomposition supports the solution to the inverse problem at hand. We then lay out our 3D cloud tomography method where we target three basic microphysical properties, volumetrically. Necessary but tedious mathematical details are presented in the Appendix A. Subsequently, the new 3D cloud tomographic capability is demonstrated on realistic synthetic clouds from an LES that provide ground truth for unambiguous retrieval error quantification. We conclude with a summary of our results and an outline of future developments, mostly looking toward CloudCT and other future space-based uses.

2. Theoretical Background

This section describes bulk microphysical parameterization of scattering media, the polarimetric radiative transfer image formation (forward) model and the relation between them. The section also describes the coordinate systems in use (per-scatterer, imager and Earth frames). We further decompose the polarized radiance into single-scattered and high-order scattered components. These foundations are used in subsequent sections, to formulate tomographic recovery. For concise reading, this section does not elaborate on some definitions. An extended, self-contained theoretical background section is given in Appendix A, which the readers are encouraged to follow to ease reproduction of our results.

2.1. Scatterer Microphysical Properties

In the lower atmosphere, cloud particles are droplets of liquid water that are very nearly spherical, having radius r. They are however polydisperse, with a droplet size distribution denotedn(r). For most remote-sensing purposes,n(r)is parameterized using an effective radius,_re(inμm) and a dimensionless variance,_ve [66]. A commonly used parametric size distribution with empirical support [66] is the Gamma-distribution (Figure 2), detailed in Appendix A.1. Let_ρwbe the density of liquid water. An important cloud characteristic is the water mass density or Liquid Water Content (LWC) per unit volume. The parameters_reand_ve dictate the shape of the droplet size distribution. LWC relates to droplet size and concentration, thus to optical opacity. In Section 3 the parametersLWC,_re,_veare unknown quantities, which we seek to estimate on a 3D grid.

2.2. Single Scattering of Polarized Light

In the absence of molecular absorption, material optical properties can be approximated using a single wavelength for each spectral band. This approximation, commonly used in multi-spectral remote sensing, relies on a single rendering of spectrally-averaged optical properties. It is valid if wavelength dependencies within a spectral band are weak, a condition met when narrow bands are considered. Macroscopic optical cross-sections are then expressed as averages weighted by the size-distributionn(r). We respectively denote the extinction, scattering and absorption cross-sections as

_σt(λ),_σs(λ),_σa(λ).

Note that throughout the text, dependency onλis generally omitted for simplicity; however, it is used at specific points as needed.

Scattering, as a fraction of the overall interaction [67], is expressed by the dimensionless single scattering albedo

ϖ=_{σs σt}.

The extinction coefficient (or optical density), denoted byβ , is expressed in terms of the LWC as [68]

β=LWC·_σ˜t.

Here,_σ˜tis the mass extinction coefficient (in units of^m2/g).

Letωand^ω′ be the unitary incident and scattered ray direction vectors respectively in Figure 2. Single-scattering geometry is defined by the local coordinate system of the incoming beam’s electric fields. The electric field of incoming light is decomposed into components along orthogonal directions. We set them as

_E1∝ω×^ω′,_E2∝_E1×ω.

The scattering angle isθ=^cos−1(ω·^ω′). The angular redistribution of singly-scatterred light by a volume element is defined by the macroscopic phase matrixP(θ). For spherical (or just randomly-oriented) particles, the phase-matrixP(θ) takes the following symmetric form [66]

Pθ=_p11θ_p21θ00_p21θ_p22θ0000_p33θ−_p43θ00_p43θ_p44θ,

where_p11is the (unpolarized) scattering phase-function.

It is convenient to define the state of polarized light in terms of the Stokes [66] vector I=^I,Q,U,V⊤, where I is the unpolarized intensity. The Stokes components relate to statistics of the fields_E1,_E2 mentioned in Equation (4), as described in Appendix A.2. The degrees of polarization (DOP) and linear polarization (DoLP) are respectively defined as the ratios^Q2+^U2+^V2/I,^Q2+^U2/I. The angle of linear polarization (AoLP) is1/2^tan−1(U/Q) . Referring to the matrix elements in Equation (5), in single-scattering of unpolarized incident sunlight, the DoLP of scattered light amounts to the ratio|_p21|/_p11.

Mie Scattering

Mie theory describes how light interacts with a spherical particle of size comparable toλ [69]. Denoteμ=cosθ. Scattering of the Stokes vectorIis described by the phase matrix_PMie(μ), which is fully defined by six matrix components:

_p11Mie(μ),_p12Mie(μ),_p22Mie(μ),_p33Mie(μ),_p43Mie(μ),_p44Mie(μ).

Expressions of Mie and Rayleigh scattering phase matrices are given in Appendix A.3.

Mie scattering due to water droplets peaks at specific angles. For a single droplet or monodisperse material,^PMiehas sharp scattering lobes at angles that depend on the droplet’sr/λratio. A macroscopic voxel contains droplets in a range of radii r, smoothing the scattering lobes. The smoothing effect depends on_ve (Figure 3). Two angular domains that stand out for remote-sensing purposes are the cloud-bow (θ∈[^135°,^155°]) and glory (θ∈[^175°,^180°]). Both domains have peaks that are sensitive to the droplet microphysical parameters, and are significantly polarized (i.e., peaks are visible in the_p12Miecomponent). The latter fact renders these peaks distinguishable in the presence of a multiply-scattered signal component.

2.3. Multiple Scattering of Polarized Light

The Radiative Transfer Equation (RTE) [70] describes multiple scattering interactions of monochromatic partially polarized light within a medium. Transmittance between two points_x1,_x2is

T_x1→_x2=exp−_∫x1x2 β(x)dx.

An atmospheric domainΩis bounded by∂Ω. The intersection of∂Ωwith a ray originating at pointxin direction−ω (Figure 4) is denoted_x0(x,ω). Denote the Stokes vector field asIx,ω. Then I(_x0,ω)is the Stokes vector of radiation which propagates in directionωat boundary point_x0(x,ω) . The non-emissive forward RT model [70] couplesIx,ωto a vector source fieldJx,ω (Figure 4) by

Ix,ω=I(_x0,ω)T_x0→x+_∫x0xJ(^x′,ω)β(^x′)T^x′→xd^x′,

Jx,ω=ϖ(x)4π_∫4πPx,ω·^ω′Ix,^ω′d^ω′.

Equations (8) and (9) are solved numerically, either directly with an explicit solver [60] or indirectly using a Monte-Carlo path tracer [71]. We use vSHDOM [60] to simulate scattered Stokes components of a realistic atmosphere, having both Mie and Rayleigh scattering due to water droplets and air molecules.

Multiple scattering interactions are defined using two coordinate systems. Local scatterer coordinates are set by(_E^1,_E^2). Stokes measurements in satellites, however, are defined in Meridional coordinates. Letz^denote the zenith direction vector at every point on Earth. In meridian coordinates, the electric field components are defined by direction vectors

_m^1=z^×ω∥z^×ω∥,_m^2=ω×_m^1.

Each pixel-scale Stokes measurement is described by a coordinate system defined by_m^1and_m^2. The transformation between the two coordinate systems amounts to a multiplication ofIby a Mueller rotation matrix.

SamplingIx,ωat the location of each camera and direction of each camera pixel yields the measured Stokes vector. A measurement k is taken at the camera position_xk, LOS direction_ωk, and wavelength_λk (Figure 4). Thus, Equations (8) and (9) yield the pixel measurement model

I[k]=I_x0,_ωkT_x0→_xk+_∫x0xk J^x′,_ωkβ(^x′)T^x′→_xkd^x′.

2.4. Single-Scattering Separation

It is often convenient to separate the single-scattering contribution from the rest of the radiance field [72]. The solar irradiance at the top of the atmosphere (TOA) is_FSun. It is unoplarized, thus corresponds to a Stokes vector_FSun=^{_FSun,0,0,0⊤}. The Sun is modeled as an ideal directional source with direction_ωSun. A solar ray heading to pointxintersects the TOA at point_xSun. The solar transmittance is given byT_xSun→x. Letδdenote Dirac’s delta. Thus,Ican be written as a sum of the diffuse component_Id, and direct solar component:

Ix,ω=_Idx,ω+δω−_ωSunFSunT_xSun→x.

Inserting (12) into (9) yields

Jx,ω=_Jdx,ω+ϖ(x)4πPx,ω·_ωSunFSunT_xSun→x,

where

_Jdx,ω=ϖ(x)4π_∫4πPx,ω·^ω′_Idx,^ω′d^ω′.

Consider Figure 4. It shows a broken-ray path of direct sunlight, which undergoes single scattering at^x′, then reaches the camera at_xk. Denote this broken ray by

_xSun→^x′→_xk.

The direction vector between^x′and_xkis_ωk. This broken ray projects to a pixel corresponding to_ωkat camera location_xkthus contributing to the measurementI(_xk,_ωk) . Using Equations (8) and (13), the single-scattered contribution from^x′is

_ISinglexSun→^x′→_xk=ϖ(^x′)4πβ(^x′)P^x′,_ωk·_ωSunFSunT_xSun→^x′T^x′→_xk.

Thus, the entire single-scattered signal accumulates contributions along the LOS

_ISinglexk=_{∫x0xk ISinglexSun}→^x′→_xkd^x′.

2.5. Ray Tracing

Ray tracing computes a function over a straight line through a 3D domain. A common operation is path-integration (e.g., Equations (7) and (8)). Leth(x)be a continuous field. Define a grid of discrete points_xg,whereg=1,2,…,_Ngrid. Denoteh[g]=h(_xg). A path-integral overh(x)is numerically computed using an interpolation kernel K

_∫x1x2 h(x)dx=∑g=1_Ngridh[g]_∫x1x2 Kx−_xgdx.

For zero-order interpolation (i.e., voxel grid), (18) degenerates to

_∫x1x2 h(x)dx=∑g=1_Ngridh[g]_ℓgx1→_x2,

where_ℓgx1→_x2 is the intersection of the path with voxel g (Figure 4). For voxel indices g that do not intersect the path_x1→_x2, the value of_ℓgx1→_x2is 0.

3. Cloud Tomography

So far, we described the forward (image-formation) model, that is, how images are formed, given cloud properties. In this work, we formulate a novel inverse tomographic problem of recovering the unknown cloud microphysical properties, volumetrically. In voxel g, the vector of unknown parameters isLWC[g],_re[g],_ve[g]. The unknown microphysical parameters are concatenated to a vector of length3_Ngrid

Θ=^{…,LWC[g],_re[g],_ve[g],…⊤},1≤g≤_Ngrid.

Neglecting circular polarization, each pixel measures a Stokes vector,_yI=_yI,_yQ,_yUat_Nλwavelengths. Let_Nviewsand_Npixdenote the number of view points and camera pixels. The total number of Stokes measurements is thus_Nmeas=_{Nλ Nviews Npix}. The measurement vector of length3_Nmeasis expressed as

y=^{_yI[1],…,_yI[_Nmeas]⊤}.

In this section, we formulate the use of measurementsy(multi-view, multi-pixel, multi-spectral, polarimetric measurements) for tomographic retrieval ofΘ (3D volumetric cloud density and microphysics). It is worth mentioning at this point that the Stokes components are not measured directly. Rather, they are computationally retrieved from measurements taken using different states of polarization filters (see Appendix C for the AirMSPI measurement model).

3.1. Polarimetric Information

To make an initial assessment of the sensitivity of polarimetric measurements, we simulate a simple homogeneous cubic cloud (Figure 5), parameterized by two microphysical parameters:(LWC,_re) . Back-scattered Stokes measurements are taken at the TOA for the same angles and wavelengths as sampled by the Airborne Multi-angle Spectro-Polarimetric Imager (AirMSPI) [14]. AirMSPI is used here as a proxy instrument to demonstrate our approach. However, the approach is applicable to other polarimetric imagers.

LetI[k],U[k],Q[k]be simulated Stokes components at measurement index k, and let_vebe held constant. Define a cost (data-misfit) function for each of the Stokes components

_DILWC,_re=∑k=1_Nmeas^{I[k]−_yI[k]2},

_DQLWC,_re=∑k=1_Nmeas^{Q[k]−_yQ[k]2},

_DULWC,_re=∑k=1_Nmeas^{U[k]−_yU[k]2}.

Equations (22)–(24) are 2D manifolds. Figure 6 plots the cost manifolds for different solar azimuth angles,_ϕ0. While there is an ambiguity betweenLWCand_rewhen relying on_DI, there are better defined minima for_DQand_DU. This indicates that polarization measurements carry valuable information.

3.2. Inverse Problem Formulation

Denote_IΘas the image formation model. Tomography can be formulated as minimization of a data-fit function. We perform

Θ^=argminΘD_IΘ,y=argminΘ^{_IΘ−y⊤ Σ−1}_IΘ−y,

HereΣis related to the co-variance matrix of the measurement noise. For brevity, we omit the subscriptΘbut remember that

I≡_IΘ,J≡_JΘ,β≡_βΘ,ϖ≡_ϖΘ,P≡_PΘ,T≡_TΘ.

Assuming noise in different pixels, wavelengths and angles is uncorrelated, (25) degenerates to

Θ^=argminΘ∑k=1_Nmeas^{I[k]−_yI[k]⊤ R−1}I[k]−_yI[k].

The matrixRdepends on the particular sensor technology. Description ofR , tailored to the AirMSPI sensor, is detailed in Appendix C.

We solve (27) by a gradient-based approach. The gradient with respect to the unknown parametersΘis

_∇ΘD_IΘ,y=2∑k=1_Nmeas^{I[k]−_yI[k]⊤ R−1} _∇ΘI[k].

The term_∇ΘI[k] is the Jacobian of the sensing model. Equation (28) is used to formulate an update rule for an iterative optimization algorithm

_Θb+1=_Θb−_{χb ∇Θ}D_IΘ,y,

where b denotes the iteration index and_χb is a scalar. We use L-BFGS [73] for numerical optimization that, in particular, determines the value of_χbadaptively. One approach to computing the gradient_∇ΘD is the Adjoint RTE [74,75]. Due to the recursive nature of the RTE, computing the gradient through the exact Jacobian_∇ΘI[k]is computationally expensive. In the following sections, we derive a method to make the computation of the gradient tractable and efficient. We do that by approximating the Jacobian_∇ΘI in a tractable way, using a two-step iterative algorithm [61,63].

3.3. Iterative Solution Approach

We formulate an iterative algorithm which alternates between two steps (See the diagram in Figure 7).

Starting with an initial guess,_Θ0 , Step 1 uses vSHDOM to compute the forward (recursive) RT equations. This renders synthetic images according to the multi-view geometry, spectral bands and spatial samples of the cameras. Keeping_Id fixed, Step 2 efficiently computes an approximate gradient with respect toΘ. The approximate gradient is fed into an L-BFGS step to update the current estimate _Θb.

Step 1: RTE Forward Model

The first step in the estimation approach is running the forward model in Equations (8) and (9) using a numerical RTE solver. This requires transforming microphysical to optical properties at every voxel (g) and spectral band (λ):

LWC[g],_re[g],_ve[g]⟶_βλ[g],_ϖλ[g],_Pλ[g].

This transformation can be implemented as described in Appendix A, using specifically Equations (A8)–(A13). However, using this implementation during each optimization iteration can be time-consuming. Therefore, we use pre-computed lookup tables,_σ˜λre,_ve,_ϖλre,_ve,_Pλre,_ve, over the grids_re∈_remin,…,_remaxand_ve∈_vemin,…,_vemax. With these pre-computed tables and{LWC[g],_re[g],_ve[g]} , vSHDOM [60] renders the Stokes vector at each 3D voxel and direction. This is the forward modeling procedure. The result is the set of fieldsIx,ω,_Idx,ω,Jx,ω,_ISinglexSun→^x′→_xk.

Step 2: Approximate Jacobian Computation

The forward vRTE model in (11) depends on optical propertiesβ,ϖ,P , which themselves depend on the sought microphysics. The Jacobian at voxel g is expressed by applying the chain-rule to (11). For example, the derivative with respect to the effective radius is

∂I[k]∂_re[g]=∂I[k]∂β[g]∂β[g]∂_re[g]+∂I[k]∂ϖ[g]∂ϖ[g]∂_re[g]+∂I[k]∂P[g]∂P[g]∂_re[g].

Analogously, replacing_re in (31) withLWCor_veyields the respective microphysical derivatives. We proceed by expressing the derivatives∂{β,ϖ,P}/∂{LWC,_re,_ve}. Afterwards, we expand and combine the derivatives∂I/∂{β,ϖ,P} to express (31).

For each voxel, the derivatives ofβ,ϖ,Pwith respect to the microphysics are calculated using pre-computed tables

∂β∂LWC=σ˜(_re,_ve),∂β∂_re=σ˜(_re+_εre ,_ve)−σ˜(_re,_ve)_εre ,

∂ϖ∂LWC=0,∂ϖ∂_re=ϖ(_re+_εre ,_ve)−ϖ(_re,_ve)_εre ,

∂P∂LWC=0,∂P∂_re=P(_re+_εre ,_ve)−P(_re,_ve)_εre ,

where_vederivatives are computed analogously to the_rederivatives. Using the shorthand notation_∂g≡{∂∂LWC[g],∂∂_re[g],∂∂_ve[g]}, the overall Jacobian is given by a sum of terms

_∂gI[k]=_A1+_A2+_A3+_A4+_A5+_A6.

The full expression for each term in Equation (35) is given in Appendix B. For example,

_A1=−_ℓgx0→_xkI_x0,_ωkT_x0→_xk∂gβ.

Let us focus on the term

_A4=_∫x0xk ϖ(^x′)4π_∫4πP^x′,_ωk·^ω′_∂gI^x′,^ω′d^ω′β(^x′)T^x′→_xkd^x′.

This Jacobian term stands out, because it is the only term which requires computing the derivative ofI. This derivative is computationally expensive becauseIis computed recursively through the RTE (Equations (8) and (9)). In principle, a change in the microphysics of one voxel can recursively affect the radiance at every other voxel. We decompose_∂gI using the diffuse-direct decomposition of (12)

_∂gI^x′,^ω′=_{∂g Id}^x′,^ω′+δ^ω′−_{ωSunFSun∂g}T_xSun→^x′.

At the core our approach for computational efficiency is the assumption that the diffuse light_Idis less sensitive to slight changes in the microphysical properties of any single voxel g. Rather,_Id is impacted mainly by bulk changes to the over-all volume. Thus, we approximate (37) by keeping_Idindependent ofΘfor a single iteration of the gradient computation, that is,

_{∂g Id}≈0.

This bypasses the complexity of recursively computing_{∂g Id}.

It is important to note that at every iteration, the Jacobian_∇ΘI[k]still is impacted by_Id. This is because_IdaffectsI through Equation (12), andIappears in the terms_A1,…,_A6. As the estimated medium properties evolve through iterations, so does_Id (in Step 1, above). We just assume during Step 2 that_{∂g Id} is negligible compared to other terms in Equation (35).

Contrary to_Id , the single-scattered component is highly sensitive to changes in the micro-physical properties of a single voxel. We therefore include an exact treatment of single-scattering in the gradient computation (Appendix B). This is the essence of our numerical optimization approach. It enables tackling multiple-scattering tomography, in practice. Simulation results presented in the following section rely on additional numerical considerations (e.g., initialization, preconditioning, convergence criteria), which are all described in the accompanying Appendix D.

4. Computational Efficiency

Section 3.2 and Section 3.3 describe our iterative gradient-based approach to Equations (28) and (29). The approach relies on an approximate Jacobian. We argue that our gradient computation is efficient, relative to the state of the art. Recall that the goal is estimation of three unknown cloud-droplet parameters per voxel. Each Jacobian component expresses differentiation relative to each parameter in each voxel. Accounting for the need to run the forward model in Equation (28) and the Jacobian, brute force computation requires3_Ngrid+1full 3D vRTE computations, for each iteration of the optimization.

An alternative tractable approach relies on adjoint differentiation of 3D vRTE. To the best of our knowledge, only Martin et al. [75] have at least formulated a potentially efficient atmosphere/surface tomography based on the adjoint 3D vector RTE. However, the method was only demonstrated using 2D scalar RT, showing promise in retrieving cloud extinction [76]. The adjoint approach requires computing one full 3D vRTE and one full adjoint-3D vRTE, per gradient computation. Due to symmetry, the overall complexity of each iteration of the optimization would be that of two 3D-vRTE computations.

Our method to approximate the gradient per iteration requires lower resources: one forward full 3D vRTE (for Equation (28)) computation and a simple single-scattering computation. Intuitively, our approach compromises gradient accuracy at voxels where the adjoint Stokes vector has undergone multiple scattering events. These voxels will tend to be at larger optical distances from the sensors. The scattering events cause a decorrelation of the adjoint Stokes vector from the optical properties and therefore loss of sensitivity to the properties of individual voxels. The measurement sensitivity to cloud properties is therefore expected to diminish. Thus, while the accuracy of the gradient in these voxels is reduced, so is the magnitude of the gradient. This region of diminished sensitivity coincides with the “hidden core” of the cloud. We refer the reader to Forster et al. [77] for more details on the hidden core concept.

5. Simulations

In the context of tomography, simultaneity is defined by the evolution timescales of targeted clouds and by the time it takes to acquire multi-angular data. Becasue updraft timescales of shallow cumulus are∼15min [78], our target clouds significantly evolve over a typical single-platform observation time span (e.g., 7 min for Terra’s MISR). Hence, single-platform systems do not satisfy simultaneity in out context. As mentioned, real data of simultaneous spaceborne multi-angular polarimetric images of clouds does not yet exist, but a mission to supply this data is in the works. Therefore, we use careful simulations to test the approach. We simulate an atmosphere with molecular Rayleigh scattering and liquid water clouds. Rayleigh scattering is calculated from the Air Force Geophysics Laboratory (AFGL) database [79] for a summer mid-latitude atmosphere. Mie tables are pre-computed for_re∈[4,25]μm and_ve=0.1with_Nre =100. The surface is Lambertian with a water-like albedo of0.05 . For realistic complexity, a Large Eddy Simulation (LES) model [80] was used to generate a cloud field. Each voxel is of size 20 × 20 × 40 m³ . The LES outputs [80] are clouds with 3D variableLWCand 1D (vertically) variable_re . A typical value [81] of_ve=0.1was chosen. Consequently, the present recovery demonstrations recoverLWCand_reon their respective native LES grid. On the other hand,_ve=0.1is excluded from the unknowns. The assumption of fixed_ve is not fundamental to our approach or implementation, however, relaxing this assumption warrants an investigation which will be done in future work. The spatial variability of the recovered parameters is discussed in more detail in Section 6.

From the generated cloud field, two isolated cloudy regions are taken for reconstruction:

• Scene A: An atmospheric domain of dimensions 0.64 × 0.72 × 20 km³ with an isolated cloud (see synthetic AirMSPI nadir view in Figure 8).

• Scene B: An atmospheric domain of dimensions 2.42 × 2.1 × 8 km³ with several clouds of varying optical thickness (see synthetic AirMSPI nadir view in Figure 9).

Synthetic measurements rendered with the spatial resolution and angular sampling of AirMSPI [14], namely, 10 m pixels and 9 viewing angles:±^70.5°,±^60°,±^45.6°,±^26.1°, and^0°from zenith, where ± indicates fore- and aft-views along the northbound flight path. The solar zenith angle is^15°from nadir in the measurement plane, that is,^0°solar azimuth. We simulate measurements at AirMSPI’s three polarized spectral bands:λ=0.47,0.66,0.865μm. The bandwidths are narrow enough (≈46 nm) to render images using a single representative wavelength per band.

Single scattering albedos for these wavelengths equal unity (up to 4 significant digits). In other words, and in sharp contrast with the operational Nakajima–King [9] two-wavelengths non-tomographic retrieval, absorption by droplets plays no role in this demonstration of the tomography of cloud microphysics. The measurements are synthesized with realistic noise, according to the AirMSPI data acquisition model (see Appendix C).

Qualitative volumetric results of the recovered LWC for Scene A are shown in Figure 10. Scatter plot of the recovered LWC and the recovery results of_re for Scene A are given in Figure 11. Qualitative volumetric results of the recovered LWC for Scene B are shown in Figure 12. A scatter plot of the recovered LWC and the recovery results of_re for Scene B are given in Figure 13.

For quantitative assessment of the recovery, we use local mean errorϵ , and global bias measures [53]ϑ:

_ϵLWC=∥LWC^_{−LWC∥1∥LWC∥1},_ϑLWC=∥LWC^_∥1−_{∥LWC∥1∥LWC∥1},_ϵre =∥_re^−_{re ∥1}∥_{re ∥1}.

We further define the Root Mean Square (RMS) error and bias of the retrieved LWC as:

RMS=1_{Ngrid∥LWC^−LWC∥21},bias=1_Ngrid∑g=1_Ngrid(LWC^[g]−LWC¯[g]).

The quantitative error measures upon convergence for the two scenes are:

Scene A:_ϵre ≈11%,_ϵLWC≈30%,_ϑLWC≈−4%,RMS≈0.065g/^m3,bias≈−0.002g/^m3,

Scene B:_ϵre ≈13%,_ϵLWC≈29%,_ϑLWC≈−5%,RMS≈0.085g/^m3,bias≈−0.0035g/^m3.

Using a 2.50 GHz CPU, the recovery run-time of cloud properties in Scenes A,B was∼13h and∼10days, respectively.

Multi-angular tomographic retrieval enables vertical resolution of the droplet effective radius. By contrast, a homogeneous droplet radius is typically retrieved by mono-angular observations fitted to a plane-parallel homogeneous cloud model. The retrieval errors of droplet radii in the demonstrations above are significantly smaller than retrieval errors of a homogeneous droplet radius. The latter can easily exceed 50% in similar conditions to our study that is, shallow cumuli and illumination conditions (see e.g., Reference [82]).

6. Discussion

In the formulation of our retrieval, and in nature, the droplet effective radius_reand droplet effective variance_vevary in 3D. This allows us to, in principle, localize the variations in cloud microphysical properties which is a significant advance over other methods. Operational remote sensing algorithms, which rely on 1D RT and plane-parallel cloud models, retrieve a single value for_re(and for_ve ), for each cloudy pixel. This occurs both in bi-spectral [9,83] and polarimetric [11,84] techniques. The physical interpretation is complicated by the difficulty in determining which portion of the cloud the retrieved quantities corresponds to.

In our current demonstrations we assume 1D variable_reand a completely homogeneous_ve. This is much more general than operational methods, but still a degenerate representation compared to full 3D variability. We now discuss the realism of such a representation.

Vertically varying droplet size distributions arise as a result of moist adiabatic thermodynamics in vertically stratified environments (e.g., Reference [81]). Further 3D variability forms due to the diabatic mixing of air parcels in the turbulent processes of entrainment and detrainment. Despite the complexity of these processes, there is evidence that the horizontal variability_reis indeed small over a cloud scale.

This evidence comes from in-situ aircraft observations of shallow cumulus [85,86,87], modelling studies [88] and theory [89]. However, there are also select observations of monsoonal clouds [90] and theoretical arguments [89] that suggest there is a sharp gradient in the droplet effective radius in the very outer shell of the clouds. If this is the case, then a representation having vertical-only variation of_reloses validity at the outer shell.

The value of_ve can also vary significantly across different environmental conditions. This is seen in research flights including in-situ measurements [91,92,93,94,95]. Moreover, in LES simulations of shallow cumulus clouds with bin microphysics,_ve might range from 0.01 to 0.26 [96]. The core of a cloud tends to have a low effective variance as condensation is the dominant process there [96,97]. Cloud edges, in contrast, experience also evaporation and entrainment mixing, as the cloud is diluted by environmental air [98]. This tends to increase_ve. If the cloud has precipitation, spatial variability of_ve increases [99].

While horizontal variations reduce the accuracy of our demonstrated retrieval, these points show that the approximations used in our demonstrations are often reasonable and a significant advance. To further push the potential of the general framework proposed here and retrieve additional degrees of freedom describing microphysical variations in 3D, we may include more information. One option is to include additional measurements, for example, by using a combination of AirMSPI [14] and Research Spectro-Polarimeter (RSP) [100]. Another option is to introduce regularization schemes, which mathematically express the natural trends of horizontal variability mentioned above. The 3D tomographic approach presented in the paper offers a lot more flexibility than current operational analyses and enables the development and application of retrievals that might provide new insight into cloud microphysics.

7. Summary & Outlook

We derive the tomography of cloud microphysics based on multi-view/multi-spectral polarimetric measurements of scattered sunlight. This novel type of tomography uses, for the first time, 3D polarized RT as the image formation model. The recovery of droplet sizes is made possible due to the crucial phase-function information carried by the polarized signal. We define a model-fitting error function and compute approximate gradients of this function to make the recovery tractable. Demonstration are done on synthetic 3D clouds, based on a Large Eddy Simulation with the effective radius assumed to vary only vertically. Our open-source code is publicly available [101].

Future work will address the extent to which polarimetric measurements penetrate the cloud and the relation between_rein the outer shell and_re in the cloud core, as defined by Forster et al. [77]. Furthermore, we will relax the fixed_ve assumption that was used in the demonstrations, and thus assess full microphysical retrieval capabilities of polarization measurements. Moreover, future plans include experimental demonstration and use, while the CloudCT formation orbits. Research is ongoing [77] about the adaptation of our cloud tomography methodology for satellite data as can be obtained from the multi-view imaging from MISR on Terra and a SWIR view from the collocated MODIS, but with significantly lower spatial resolution than used so far.

Lastly, we note that our atmospheric tomography approach described herein can be adapted to aerosols, including dense plumes of wild fire smoke, volcanic ash, and dust near their sources. The main difference in this extended application is that the forward model must be further linearized with respect to the complex refractive index of the varying aerosol material, which has been done for 1D polarized RT (e.g., Reference [102]), even for non-sphericity effects (e.g., Reference [103]).

Author Contributions

Conceptualization, A.L. and Y.Y.S.; methodology, A.L., Y.Y.S., A.B.D.; software, A.L., J.L.; validation, A.L. and J.L.; formal analysis, A.L.; investigation, A.L.; resources, Y.Y.S.; data curation, A.L. and J.L.; writing-original draft preparation, A.L.; writing-review and editing, A.L., Y.Y.S., A.B.D. and J.L.; visualization, A.L.; supervision, Y.Y.S.; project administration, A.L. and Y.Y.S.; funding acquisition, Y.Y.S. All authors have read and agreed to the published version of the manuscript.

Funding

The authors are grateful to the US-Israel Binational Science Foundation (BSF grant 2016325) for continuously facilitating our international collaboration. Aviad Levis' work was partially supported as a Zuckerman Foundation STEM Leadership Fellow. Yoav Schechner is a Landau Fellow supported by the Taub Foundation. His work was conducted in the Ollendorff Minerva Center (BMBF). Anthony Davis' work was carried out at JPL/Caltech, supported by NASA's SMD/ESD/(RST+TASNPP) and ESTO/AIST programs. That research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). Support for Jesse Loveridge's work from JPL under contract #147871 is gratefully acknowledged. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 810370: CloudCT).

Acknowledgments

We thank I. Koren, D. Rosenfeld, A. Aides, D. Diner, L. Di Girolamo, and G. Matheou for support and fruitful discussions. We acknowledge F. Evans and A. Doicu for the online vSHDOM code.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Extended Theoretical Background

Appendix A.1. Scatterer Microphysical Properties

In the lower atmosphere, cloud particles are droplets of liquid water that are very nearly spherical, having radius r. They are however polydisperse, with a droplet size distribution denotedn(r). For most remote-sensing purposes,n(r)is parameterized using an effective radius inμ m and a dimensionless variance [66]:

_re=_∫0∞ ^r3n(r)dr_∫0∞ ^r2n(r)dr,_ve=_∫0∞ ^{r-_re2 r2}n(r)dr_{re2 ∫0∞} ^r2n(r)dr.

A commonly used parametric size distribution with empirical support [66] is the Gamma-distribution (Figure 2):

n(r)=Nc^r(_ve-1-3)exp[-r/(_{re ve})],

where we require0<_ve<1/2. Herec=^{(_{re ve})(2-_ve-1)}/Γ(_ve-1-2)is a normalization constant and

N=_∫0∞n(r)dr

is the droplet number concentration. Let_ρwbe the density of liquid water. An important cloud characteristic is the water mass density or Liquid Water Content (LWC) per unit volume:

LWC=43π_{ρw ∫0∞} ^r3n(r)dr.

It is expressed as LWC =4/3π_{ρw re3}(1-_ve)(1-2_ve) for the Gamma distribution in (A2). In Section 3 these microphysical parameters,LWC,_re,_ve, are unknown quantities we seek to estimate on a 3D grid.

Appendix A.2. Polarized Light

A light wave is associated with orthogonal components of a random electric wave,_E1(t)and_E2(t), where t is time. The components' direction unit vectors are respectively_E^1and_E^2. The wave propagates in directionω=_E^1×_E^2 . It is convenient to define the polarized light state in terms of the Stokes [66] vectorI=^I,Q,U,V⊤. Each component ofIexpresses temporal integration:

I=_{〈E1 E1*+E2 E2*〉t},Q=_{〈E1 E1*-E2 E2*〉t},U=_{〈E1 E2*+E2 E1*〉t},V=i_{〈E1 E2*-E2 E1*〉t},

wherei=-1. Unpolarized intensity is I. The degrees of polarization (DOP) and linear polarization (DoLP) are respectively defined as the ratios^Q2+^U2+^V2/I,^Q2+^U2/I. The angle of linear polarization (AoLP) is1/2^tan-1(U/Q).

Appendix A.3. Single Scattering of Polarized Light

Light interaction with a single particle is described by the total extinction cross-section_st(r,λ), decomposed into scattering and absorption cross-sections, respectively:

_st(r,λ)=_ss(r,λ)+_sa(r,λ).

In Mie scattering by spheres, introduced further on, we have

_st(r,λ)=π^r2_Qt(2πr/λ),_ss(r,λ)=π^r2Qs(2πr/λ),_sa(r,λ)=π^r2Qa(2πr/λ),

whereQt,Qs,Qaare dimensionless efficiency factors, which depend on the normalized size parameter2πr/λ. In the limitr≫λ,Qt≈2. Furthermore, when_ss(r,λ)≫_sa(r,λ), thenQs≈2andQa≈0.

When molecular absorption can be neglected, material optical properties can be approximated using a single wavelength for each narrow spectral band. This approximation, commonly used in multi-spectral remote sensing, relies on a single rendering of spectrally-averaged optical properties. It is valid if wavelength dependencies within a spectral band are weak, a condition met when narrow bands are considered. Macroscopic optical cross-sections are then expressed as weighted averages

_σt(λ)=_{〈st(r,λ)〉r},_σs(λ)=_{〈ss(r,λ)〉r},_σa(λ)=_{〈sa(r,λ)〉r}.

We note that (A8) aggregates scattering properties rather than electric fields. This aggregation holds for scatterer populations that are in each other's far field (i.e., are≫λ apart) [66]. Here the size-weighted average over some functiona(r)is defined by

_〈a〉r=1N_∫0∞a(r)n(r)dr.

The size integral of (A9) is terminated in practice at_rmax= 70μm. Note that throughout the text, dependency onλis generally omitted for simplicity; however, it is used at specific points as needed.

Scattering, as a fraction of the overall interaction [67], is expressed by the dimensionless single scattering albedo

ϖ=_{σs σt}.

The extinction coefficient (or optical density) is denoted byβ . Following Equations (A3), (A4) and (A8),β=N_σt is expressed in terms of the LWC as [68]

β=LWC43π_{ρw 〈^r3〉rσt}=LWC·_σ˜t.

Here,_σ˜tis the mass extinction coefficient (in units of^m2/g).

Letωand^ω' be the unitary incident and scattered ray direction vectors respectively in Figure 2. Single-scattering geometry is defined by the local coordinate system of the incoming beam's electric fields. As stated above, the electric field of incoming light is decomposed into components along orthogonal directions. We set them as

_E1∝ω×^ω',_E2∝_E1×ω.

The scattering angle isθ=^cos-1(ω·^ω'). The angular redistribution of singly-scattering light from a sphere of is defined by the4×4dimensionless Mueller matrix_Ps(θ,r). The macroscopic phase matrix is the size-weighted average

P(θ)=_{〈ss(r)Ps(θ,r)〉r σs}.

For spherical (or just randomly-oriented) particles, the phase-matrixP(θ) takes the following symmetric form [66]

Pθ=_p11θ_p21θ00_p21θ_p22θ0000_p33θ-_p43θ00_p43θ_p44θ,

where_p11is the (unpolarized) scattering phase-function. In single-scattering of unpolarized incident sunlight, the DoLP of scattered light amounts to the ratio|_p21|/_p11.

Appendix A.3.1. Rayleigh Scattering

The Rayleigh model describes light scattering by particles much smaller than the wavelength. The Rayleigh phase matrix takes the following form [70]

_PRaylθ=341+^cos2θ-34^sin2θ00-34^sin2θ341+^cos2θ000032cosθ000032cosθ.

The single-scattering DoLP due to air molecules is then

_DoLPRayl(θ)=^sin2θ1+^cos2θ.

According to (A16) a maximum DoLP is attained at single-scattering angleθ=^90°.

Appendix A.3.2. Mie Scattering

Mie theory describes how light interacts with a spherical particle of size comparable toλ [69]. Denoteμ=cosθ. Mie scattering is defined in terms of complex-valued amplitude scattering functions_S1(μ),_S2(μ), which correspond to scattering of the_E1,_E2electric field components. Scattering of the Stokes vectorIis described by the phase matrix_PMie(μ), which is fully defined by six matrix components:

_p11Mie=ϱ2_{S1 S1*}+_{S2 S2*},_p12Mie=ϱ2_{S1 S1*}-_{S2 S2*},_p22Mie=ϱ2_{S1 S1*}+_{S2 S2*},_p33Mie=ϱ2_{S1 S2*}+_{S2 S1*},_p43Mie=ϱ2_{S1 S2*}-_{S2 S1*},_p44Mie=ϱ2_{S1 S2*}+_{S2 S1*}.

Here,ϱis a normalization constant, set to satisfy12_{∫-11 p11Mie}(μ)dμ=1.

Mie scattering due to water droplets peaks at specific angles. For a single droplet or monodisperse material,^PMiehas sharp scattering lobes at angles that depend on the droplet'sr/λratio. A macroscopic voxel contains droplets in a range of radii r, smoothing the scattering lobes. The smoothing effect depends on_ve (Figure 3). Two angular domains that stand out for remote-sensing purposes are the cloud-bow (θ∈[^135°,^155°]) and glory (θ∈[^175°,^180°]). Both domains have peaks that are sensitive to the droplet microphysical parameters, and are significantly polarized (i.e., peaks are visible in the_p12Miecomponent). The latter fact renders these peaks distinguishable in the presence of a multiply-scattered signal component.

Appendix B. Jacobian Derivation

In Equation (35) of the main text, the Jacobian is written as a sum of six terms

_∂gI[k]=_A1+_A2+_A3+_A4+_A5+_A6.

In this section we expand and describe each of these terms. Using Equations (7) and (19), the transmittance derivative is

_∂gT_x1→_x2=-T_x1→_x2ℓgx1→_x2∂gβ.

Then,

_A1=-_ℓgx0→_xkI_x0,_ωkT_x0→_xk∂gβ,

_A2=_ℓgx0→_xk∫_{x0 xk∂g}ϖ4π∫4πP^x',_ωk·^ω'I^x',^ω'd^ω'β(^x')T^x'→_xkd^x',

_A3=_ℓgx0→_xk∫_{x0 xk}ϖ(^x')4π∫4π_∂gP^x',_ωk·^ω'I^x',^ω'd^ω'β(^x')T^x'→_xkd^x',

_A4=∫_{x0 xk}ϖ(^x')4π∫4πP^x',_ωk·^ω'_∂gI^x',^ω'd^ω'β(^x')T^x'→_xkd^x',

_A5=_ℓgx0→_xk∂gβ∫_{x0 xk}J^x',_ωkT^x'→_xkd^x',

_A6=-_ℓgx0→_xk∂gβ∫_{x0 xk}J^x',_ωkβ(^x')T^x'→_xkd^x'.

Note thatIx,ωandJx,ω are computed in Step 1 and are therefor ready for use when computing_A1,_A2,_A3,_A5and_A6. Furthermore,_ℓgx0→_xk=0for any voxel that is not on the LOS of pixel k. Therefore, the terms_A1,_A2,_A3,_A5,_A6are computed using a single path tracing_xk→_x0.

We now give special attention to_A4 in Equation (A23). Using the diffuse-direct decomposition of (12), we decompose (A23) as

_A4=∫_{x0 xk}ϖ(^x')4π∫4πP^x',_ωk·^ω'_{∂g Id}^x',^ω'd^ω'β(^x')T^x'→_xkd^x'+∫_{x0 xk}ϖ(^x')4π∫4πP^x',_ωk·^ω'δ^ω'-_{ωSunFSun∂g}T_xSun→^x'd^ω'β(^x')T^x'→_xkd^x'.

The first term in (A26) is based on_{∂g Id}, i.e., a derivative of the diffuse (high order scattering) component. Herein lies a recursive complexity. In principle, a differential change in the microphysics of one voxel can recursively affect the radiance at every other voxel, and this affects all the pixels. To make calculations numerically efficient, we approximate (A26). The approximation assumes that relative to other components in the Jacobian,_Idis less sensitive to a differential changes in the microphysical properties at voxel g. Thus, (A26) is approximated by keeping_Idindependent ofΘfor a single iteration of the gradient computation, that is,

_{∂g Id}≈0.

The second term in (A26) is based on differentiation of the direct component. This is straight-forward to compute using (A19). Consequently, using Equation (A27) and the definition of_ISinglexSun→^x'→_xk in (16), the term_A4in (A26) is approximated by

_A4≈_A4˜=_∂gβ∫_{x0 xkℓgxSun}→^x'_ISinglexSun→^x'→_xkd^x'.

The term_ℓgxSun→^x'in (A28) contributes to voxels outside of the LOS. The integral in_A4˜ is computed with a broken-ray [104] path_xk→^x'→_xSun , as illustrated in Figure 4.

Using Equations (8), (A19), (A20) and (A25),_A1and_A6are combined to

_A1,6=_A1+_A6=-_∂gβI_xg,_ωk.

Overall, in our iterative procedure, we approximate the Jacobian in (35) by

_∂gI[k]=_A1,6+_A2+_A3+_A4˜+_A5.

Equations (A20)-(A28) formulate the Jacobian in terms of a voxel grid (zero-order interpolation). However, in practice we use a trilinear interpolation kernel K in (18), consistent with vSHDOM internal interpolation [59].

Appendix C. Measurement Noise

The inverse problem defined in the main text is formulated in terms of measured Stokes vectors (Equaiton (21)). However, Stokes vectors are not measured directly. Rather, they are derived from intensity measurements taken through filters. The raw intensity measurements are noisy. Noise is dominated by Poisson photon noise, which is independent across different raw measurements. However, the estimation of Stokes components from independent intensity measurements yields noise which is correlated across the components of the Stokes vector, per-pixel. In this section, we describe the synthesis model we employ to generate realistic noise in simulations. Our synthesis is based on the AirMSPI [105] sensor model. Furthermore, we derive the expression forR , which we use in the recovery process (Equation (27) in the main text).

AirMSPI measures a modulated intensity signal at_Nsub=23subframes. Define a normalized frame which spans the unitless integration time intervalψ∈[-0.5,0.5]. Denote the temporal center and span of each subframe as_ψlandΔψ=1/_Nsub , respectively (Figure A1). Based on the sensing process described in Ref. [105], define the following modulation function, whose parameters are given in Table A1:

M[l]=_J0[κ(_ψl)]+13^πΔψ22 _γ02(λ)_J2[κ(_ψl)]-cos[2(π_ψl-η)]_J0[κ(_ψl)],

with

κ(_ψl)=-2_γ0(λ)sin(π_ψl-η)1+^cot2(π_ψl-η).

Remotesensing 12 02831 g0a1 550

Figure A1.A normalized frame spans the interval[-0.5,0.5], evenly divided into_Nsubsubframes.

Figure A1.A normalized frame spans the interval[-0.5,0.5], evenly divided into_Nsubsubframes.

Table

Table A1.Modulation parameters [105] used for synthesis of AirMSPI measurements.

Table A1.Modulation parameters [105] used for synthesis of AirMSPI measurements.

γ0(470nm)	γ0(660nm)	γ0(865nm)	ξ(470nm)	ξ(660nm)	ξ(865nm)	η
4.472	3.081	2.284	1.0	0.27	0.03	0.009

Here_J0,_J2are the Bessel functions of the first kind of order 0 and 2, respectively. Denote byξ(λ) a wavelength-dependent ratio, which is drawn from quantum efficiencies and spectral bandwidths of each AirMSPI band (Table A1). For the exact calculation of the ratio see Equation (24) of [105]. Using simulated Stokes vectors derived by vSHDOM, AirMSPI measurements are synthesised as passing through two polarization analyzing filters [105]. As defined in Equation (11) in the main text,I[k]is the Stokes vector in pixel k. Correspondingly the intensity isI[k], whileQ[k],U[k]are the polarized components. Measurements l through the two filters of AirMSPI are modeled by

_d0[l,k]=ξ(λ)I[k]+M[l]Q[k]

_d45[l,k]=ξ(λ)I[k]+M[l]U[k],

whereM[l]andξ(λ) are given in Equation (A31) and Table A1, respectively. The units of d are Watts. Defined[k]=^{_d0[1,k],...,_d0[_Nsub,k],_d45[1,k],...,_d45[_Nsub,k]⊤}. In matrix from, the transformations by Equations (A33) and (A34) are written using a single46×3modulation matrixM

d[k]=MI[k].

Detection is by a camera which generates photo-electrons in each pixel well. The relation between_d0[l,k]or_d45[l,k]and the expected unit-less number of photo-electrons in the pixel is given by a gain G. The number of photo-electrons is random (Poissonian) around this expected value. The vector of simulated electron counts is thus synthesized by a Poisson process

e[k]∼PoissonroundG·d[k]=PoissonroundG·MI[k].

The gain G is chosen to let the maximum signal at each camera view (i.e. maximum over pixels, wavelengths and subframe measurements) reach the maximum full-well depth of200,000 electrons, consistent with AirMSPI specifications. Equation (A36) synthesizes raw AirMSPI signals including noise (Figure 8). The synthesized AirMSPI signals, including this noise, are now used as inputs to the calculation of measured Stokes vectors in each pixel and viewpoint. The vector of electron countse[k] in each pixel k is transformed into Stokes synthetic data (Equaiton (21)) using a3×46demodulation matrixW

_yI[k]=^{(M⊤M)-1 M⊤}e[k]=We[k].

The vectors_yI[k]form the data for tomographic analysis.

Our tomographic analysis takes into account the noise properties, including noise correlation. As we now show, the measurement model (A37) yields correlated noise of different Stokes components. Thus,^R-1 (in Equation (27)) is not diagonal. Denote the diagonal co-variance matrix of the photo-electron readings by_Ce-1=diage. Let_I46×46denote the Identity matrix. The signal is generally dominated by unpolarized multiply-scattered background light. Relative to it, the magnitude of the modulated polarization signal is small. Thus, per pixel k, the diagonal matrix_Ce-1[k]is approximately constant with a global weight

_Ce-1[k]≈α[k]_I46×46.

Using Equations (A37) and (A38) for each pixel, the Stokes co-variance matrix is

^C-1[k]=^M⊤ _Ce-1[k]M≈α[k]^M⊤M.

A maximum-likelihood estimator corresponding to a Poisson process should have a weightα[k]∝1/_∥e∥1, to account for higher photon noise in brighter pixels. In simulations, however, we found thatα[k]=1 worked better. This is perhaps due to richer information carried by denser cloud regions, i.e. brighter pixels. Overall the expression we minimize in (27) is

Θ^=argminΘ∑k=1_Nmeas^{I[k]-_yI[k]⊤ M⊤}MI[k]-_yI[k],

that is,^R-1=^M⊤M.

Appendix D. Numerical Considerations

In this section we describe numerical considerations that stabilize the recovery.

Appendix D.1. Hyper-Parameters

Our code requires the choice of hyper-parameters for rendering with vSHDOM [106] in Step 1 and optimization with scipy L-BFGS [73,107] in Step 2. Table A2 summarizes the numerical parameters used in our simulations.

Table

Table A2.Numerical parameters. For vSHDOM parameter definitions, see Ref. [106]. For L-BFGS parameter definitions, see Ref. [107].

Table A2.Numerical parameters. For vSHDOM parameter definitions, see Ref. [106]. For L-BFGS parameter definitions, see Ref. [107].

vSHDOM			L-BFGS
Nμ	Nϕ	splitting accuracy	gtol	gtol	maxls
8	16	0.1	1×10-16	1×10-16	30

Appendix D.2. Preconditioning

Multivariate optimization can suffer from ill-conditioning due to different scales of the sought variables. This is expected when recovered variables represent different physical quantities with different units and orders of magnitude. A preconditioning of the update rule in (29) takes the following form

_Θb+1=_Θb-_χb ^Π-1 _∇ΘD_IΘ,y,

where we apply a diagonal scaling matrixΠ(Jacobi preconditioner) to scale the different physical variables(LWC,_re). Thus,Πtakes the form

Π=diag_ΠLWC,_Πre ,...,_ΠLWC,_Πre .

In our tests, we use_ΠLWC= 15 and_Πre = 0.01 to scale the parameters to a similar magnitude and closer to unity upon initialization.

Appendix D.3. Initialization

The recovery is initialized by the estimation of a cloud voxel mask, which bounds the cloud 3D shape. The 3D shape bound of the cloud is estimated using Space-Carving [37]. Space-carving is a geometric approach to estimate a bound to 3D shape via multi-view images. The following steps are preformed in our space-carving algorithm

1. Each image is segmented into potentially cloudy and non-cloudy pixels (we use a simple radiance threshold).

2. From each camera viewpoint, each potentially cloudy pixel back-projects a ray into the 3D domain. Voxels that this ray crosses are voted as potentially cloudy.

3. Voxels which accumulate "cloudy" votes in at least 8 out of the 9 AirMSPI viewpoints are marked as cloudy.

Outside of the shape bound, LWC = 0 throughout iterations. Within the estimated cloud-shape bound, the volume content is initialized as homogeneous withLWC=0.01g/^m3,_re=12μmand_ve= 0.1. Then, inside of the shape-bound, {LWC,_re,_ve} change throughout iterations, possibly diminishing LWC to very small values.

Appendix D.4. Convergence

Our approach alternates between Step 1 (RTE rendering) and Step 2 (approximate gradient) until convergence (Figure 7). The convergence criteria are dictated by the L-BFGS step: at each iteration, the relative change to the forward model and its gradient are compared to the ftol and gtol parameters (see Table A2 for values used). See SciPy documentation [107] for exact description of the L-BFGS stopping criteria.