Clouds cover roughly two thirds of the globe. Reflecting and absorbing solar and thermal radiation, they strongly affect the Earth’s energy (Rossow et al., 2002). The way this budget is modulated depends, among other factors, on the cloud particle size distribution. Since clouds contribute the largest uncertainty to estimates of the Earth’s changing energy budget, accurate retrievals of particle size from satellite observations is one of the most critical goals in cloud remote sensing. Indeed, cloud particle size has not only a significant influence on cloud response to aerosol modification (Várnai and Marshak, 2002; Oreopoulos and Platnick, 2008; Spencer et al., 2019), but it is also a key parameter for understanding aerosol-cloud-precipitation interactions (e.g., Albrecht, 1989; Twomey, 2007).

The current approach to simultaneously estimating cloud optical depth τ (vertical integral of cloud extinction coefficient over cloud physical thickness) and the effective particle radius r _eff (the ratio of the 3rd to the 2nd moment of particle size distribution) goes back to Twomey and Seton (1980). Their method was further developed and popularized by Nakajima and King (1990). Since then, this approach has been thoroughly studied in both theory and practice (see Platnick et al., 2003; Platnick et al., 2017 and references therein) and has been applied to a variety of observations, for example to Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) data. The approach uses reflectances from two narrow spectral bands: one in the visible (or near-infrared) region where cloud absorption is negligible, and the other in a near-infrared region where solar radiation is slightly absorbed by cloud particles. When the two reflectance measurements are combined, both τ and r _eff can be determined (Nakajima and King, 1990).

The retrievals are based on library calculations for plane-parallel homogeneous clouds (Platnick et al., 2003). The underlying assumption of plane-parallel clouds and the use of a one-dimensional (1D) radiative transfer forward model can introduce substantial errors for retrievals applied to inhomogeneous cloud scenes (e.g., Loeb and Coakley, 1998; Zuidema and Evans, 1998; Stephens and Kummerow, 2007; Evans et al., 2008; Di Girolamo et al., 2010; Wolters et al., 2010). Left and middle panels in Figure 1, reproduced from Zhang and Platnick (2011), show the global monthly means of effective radii retrieved using a visible band in combination with the 2.1 μm and 3.7 μm MODIS bands, respectively.

Figure 1 shows the dramatic difference between r _eff retrieved from 2.1 μm to 3.7 μm. As demonstrated in Zhang and Platnick (2011), 2.1 µm is affected more by inhomogeneous cloud structure than 3.7 μm, which remains relatively more stable. The difference is well correlated with cloud fraction (CF): lower CF leads to higher values of r _eff being retrieved using 2.1 µm. Moreover, the error in r _eff retrieved from 2.1 μm strongly increases with the cloud horizontal heterogeneity index (H _σ) introduced by Liang et al. (2009). This suggests the need for physically meaningful structural variables to account for cloud heterogeneity in new retrieval techniques (e.g., Fu et al., 2019).

The use of the plane-parallel assumption in droplet effective radius retrievals can result in substantial errors in inhomogeneous or broken cloud regions. Here we propose a new and completely different retrieval approach that has the potential of filling the gap in current retrieval capabilities for inhomogeneous clouds. This technique comes from nuclear reactor physics and was developed in 1960 s (e.g., Bell and Glasstone, 1970). We propose to adapt the technique for space-based cloud remote sensing. It is not based on the 1D assumption and can be applied to any inhomogeneous and broken cloud fields. Obviously, the proposed technique has its own strengths and weaknesses that will be also addressed here.

The outline of the paper is as follows. Section 2 discusses the operationally retrieved cloud droplet effective radius and cloud optical depth and then Section 3 examines the ratio of cloud reflectance over single scattering albedo. Next, Section 4 states the two main hypotheses for interpreting satellite observations. Section 5 discusses retrievals of single scattering albedo and number of scatterings. Finally, Sections 6, 7 summarize our theoretical research questions and provide a few concluding remarks.

Figure 2 illustrates the relationship between the operationally retrieved values of droplet effective radius r _eff and cloud optical depth τ for all 29,148 data points of water clouds in a MODIS Aqua cloud product granule over the South Pacific Ocean. After considering all data points (grey dots), we selected from them the 1794 data points with the fixed value of ∼0.48 for the ratio of 0.87 μm and 2.13 μm reflectances (blue dots). Out of those, we selected 103 data points where the scattering angle between solar and viewing directions is ∼99° (red dots). We note that the values of 0.48° and 99° were chosen arbitrarily, as representative values that occur abundantly in the sample MODIS granule shown in Figure 2. As we see, the red points lie along a curve that has a roughly hyperbolic shape, suggesting that the product of r _eff and τ is nearly constant.

For a sample cloudy pixel in the left panel of Figure 2, Figure 3 illustrates the ratio of the observed MODIS Top-of-Atmosphere (TOA) reflectance I _λ over the spectral single scattering albedo ω_0λ (calculated based on the retrieved effective radius r _eff of 13 μm) plotted vs. the reflectance I _λ(Ω, Ω₀) at four MODIS bands (λ = 0.86, 1.65, 2.13, and 3.75 μm). Here Ω and Ω₀ are solar and viewing directions, respectively. A linear fit indicates that the relationship I λ Ω , Ω 0   /   ω 0 λ = p   I λ Ω , Ω 0 + q Ω , Ω 0 (1) holds with a slope p = 0.943 and an offset q = 0.011.

Equation 1 can be rearranged as I λ Ω , Ω 0 = ω 0 λ   q Ω , Ω 0 / 1 − p ω 0 λ (2)

Analyses of successive order of scattering suggest that Eq. 2 represents the exact solution of the 3D radiative transfer equation with non-reflecting boundaries such as a black underlying surface (Huang et al., 2007; Knyazikhin et al., 2011; Yang et al., 2017). Here p is the recollision probability, defined as the probability that a photon scattered in the medium will collide in the medium again. The variable q Ω , Ω 0 = ρ Ω i 0 Ω 0 quantifies escape events: ρ Ω is the probability that a scattered photon will escape the medium in a given (up- or downward) direction Ω and i 0 Ω 0 is the interceptance, defined as the portion of photons in the incident beam coming from direction Ω₀ that collide with cloud droplets for the first time. The interceptance i ₀ is strongly sensitive to 3D medium structure and is related to direct transmittance: their sum is equal to 1. The fraction of intercepted photons initiates the process of photon–medium interactions while the recollision probability determines the number of scattering events as results of multiple interactions: the higher p, the more interactions the photons undergo. On average, an intercepted photon will have N scat = 1 / 1 − p ω 0 λ (3) interactions (Marshak et al., 2011).

Indeed, in the literature of vegetation canopy remote sensing, Eq. 2 is usually written as (e.g., see Knyazikhin et al., 2013; Yang et al., 2017) I λ Ω , Ω 0 = D A S F Ω , Ω 0 W λ = ρ Ω i 0 Ω 0 1 − p   [ ω 0 λ 1 − p 1 − p ω 0 λ ] , (4)

Where DASF(Ω,Ω₀) = ρ(Ω)i ₀(Ω₀)/(1 − p) is the Directional Area Scattering Factor and W _λ = ω_0λ(1 − p)/(1 − pω_0λ) is the scattering coefficient of a radiative active layer (Knyazikhin et al., 2013). Note, that in Eq. (4) the extinction coefficient and the scattering phase function, normalized by single scattering albedo, do not depend on λ . As a result, DASF becomes independent of wavelength while W _λ is almost independent of the sun-view geometry but varies with wavelength.

To complete our retrieval approach, we need to estimate the slope p from Eq. 1. For this, we combine radiances at MODIS non-absorbing (say, 0.87 μm) and water absorbing channels (1.64, 2.13, and 3.75 μm) using Twomey and Bohren’s (1980) approximation (4), which relates 1-ω_0λ at these channels to r _eff as shown in Figure 5. For any given p, the number of scatterings N _scat is related to the single scattering albedo ω₀ through Eq. 3 (see also Eq. 9 in Marshak et al., 2011). If we then add to Figure 4 the curve of N _scat as a function of ω₀, the two independent sets of curves will intersect (see Figure 7). For each spectral reflectance I _λ of a cloudy pixel, this intersection provides the two fundamental radiative transfer parameters: the number of scatterings N _scat and the single scattering albedo ω₀ averaged over all photon paths contributing to the observations. Naturally, both N _scat and ω₀ are functions of wavelength λ. Note that using ω₀ instead of r _eff will allow a more flexible interpretation of droplet size if we need to adjust the shape of the particle size distribution (whereas r _eff is tied to a specific distribution-shape).

Here we briefly summarize Figure 7, which is the key plot explaining the proposed approach. For each cloudy pixel that is associated with a fixed solar-viewing geometry, the applicable thin line is selected based on the ratio as A = 1 − I _λ=2.13/I _λ=0.87 (see Eq. 5), while the applicable thick line is specified by Eq. 1 for the spectrally-invariant slope p (the recollision probability) obtained by performing a regression based on Eq. 2 to determine what p and q values work best for a set of radiances observed at multiple wavelengths at the same location, with the ω₀ values at different wavelengths being related through Eq. 8.

First, let us point out that the plane-parallel assumption has been neither stated nor used in the above approach. Indeed, the fundamental radiative transfer theory behind this approach is based on Eqs 1, 2, which are valid for any 3D scattering and absorbing medium (Várnai and Marshak, 2002; Várnai and Marshak, 2009; Cahalan et al., 2005; Marshak et al., 2006). Even so, both the number of scatterings N _scat and the single scattering albedo ω₀ (averaged over the photon path) can be obtained for each wavelength λ (note that passive 1D reflectance-based cloud retrievals are fundamentally retrieving a 1D ω₀). Our main assumption here is that the “absorbing wavelength photons” follow the same path as the non-absorbing ones. In other words, we assumed that the solution of the radiative transfer equation, I _λ, depends on {τ_λ, P _λ} (where P _λ is the spectral phase function) in a way that does not change with wavelength, and so the wavelength dependence of I _λ comes only from the ω_0λ spectra (Marshak et al., 2012).

We note that a such separation of variables is natural for radiative transfer in leaf canopies (Knyazikhin et al., 2011), where the scattering objects are much larger than the wavelength of solar radiation—and the dependence on τ and p is determined by canopy structure (Schull et al., 2007), while the dependence on ω₀ comes entirely from leaf physiology (Knyazikhin et al., 2013). For atmospheric radiative transfer this assumption is not met, since the size of scattering objects (air molecules, aerosol and cloud particles) is comparable to (or smaller than) the wavelength of solar radiation. However, in cloudy atmospheres these assumptions can be met approximately for a wide range of wavelengths (Marshak et al., 2011).

As a proof of concept, Figure 8 illustrates for several cloud optical depths the impact of neglecting the wavelength-dependence of scattering phase functions. As expected, the impact decreases as the optical depth increases: For example, for λ = 2.13 μm and τ = 5, the errors are about 8%–12%, while for τ = 40 they are reduced to 3%–8%, depending on the Solar Zenith Angle (SZA).

Finally, we note that 1D radiative transfer theory predicts that the reflectance functionally depends on the scaled optical depth τ (1-g) (e.g., van de Hulst, 1980). Assuming that g is known a priori, the impact of phase function variations can be mitigated and the assumption of wavelength-independent phase functions can be relaxed to a certain degree. Quick preliminary calculations showed that the difference between the pair of diamond symbols for τ = 5 in the left panel of Figure 8 corresponds to a roughly 2 μm change in r _eff. The difference is even smaller in case of larger optical depths (middle and right panels in Figure 8). This bias can be much smaller than the one due to the assumption of plane-parallel geometry in case of strongly inhomogeneous clouds.

The above approach has been partially developed only for a black surface. However, if we know the reflectivity of the underlying surface, the surface albedo effects can be adequately removed (Stephens and Kummerow, 2007). We are familiar with the decomposition of the reflected radiation into the “black surface” problem and the additional radiative field due to the interaction between the underlying surface and the medium (Knyazikhin and Marshak, 2000; Knyazikhin et al., 2005).

While the value of retrieving ω₀ (and hence r _eff) is clear and may warrant using the proposed approach even for just this purpose, one may wonder, what information about cloud properties does the retrieved number of scatterings convey? First of all, for non-absorbing scattering (ω₀ = 1) in a plane-parallel medium, one can relate the average number of scatterings N _scat to cloud optical depth τ (and scattering asymmetry factor g) using Monte Carlo simulations or the diffusion approximation (Eqs 6, 7; see also Figure 3 in Kato and Marshak (2009), which relates N _scat to τ for different viewing zenith angles). For absorbing wavelengths (ω₀ < 1), the more complex formulae also include the single scattering albedo ω₀ (Platnick, 2001a; Marshak and Davis, 2005, pg. 559). Thus, under the common plane-parallel assumption, we get cloud optical depth τ and, combining τ with r _eff, we can also get liquid water path as done by current operational retrievals of, for example, MODIS (e.g., Platnick et al., 2003). Second, the average number of scatterings N _scat helps estimating horizontal photon transport; as horizontal transport can be expressed as a function of N _scat, ω₀ and g (Platnick, 2001b). Indeed, the root-mean-square displacement of reflected and transmitted photons can be estimated from diffusion theory if the average number of scatterings is known (Platnick, 2001a; Platnick, 2001b; Marshak and Davis, 2005, pg. 560–561). We note that inconsistencies between the retrieved N _scat values and the optical depths retrieved independently using the 1D Nakajima-King approach may help identify the cases where the current operational cloud products are significantly impacted by heterogeneity effects. Finally, the average number of scatterings is closely related to the average photon path length thus it is a useful quantity in its own right. We note that the average number of scatterings obey reciprocity (Platnick, 2001b), thus all results are valid upon exchange of solar and viewing directions.

In this paper we propose to retrieve two fundamental radiative transfer parameters: single scattering albedo, ω₀, and number of scatterings, N _scat. This retrieval does not require the assumption of plane-parallel geometry. The traditional cloud microphysical properties like optical depth, τ, and droplet size, r _eff, can be easily obtained from these parameters.

Eliminating the 1D assumption does not come for free here. We assume that the photons’ trajectory is the same in absorbing and non-absorbing channels. This assumption is strong and creates a bias in τ and r _eff. However, in many cases, the bias appears likely to be (much) smaller than the biases that the assumption of 1D geometry brings for strongly inhomogeneous clouds.

In addition, the retrieved number of scatterings and single scattering albedo convey broad information about cloud properties. For example, ω₀ allows a more flexible interpretation of droplet size in cases of different shapes of particle size distribution, while N _scat allows to infer the impact of photon horizontal transport.