Remote Sensing of Cloud Droplet Radius Profiles using solar reflectance from cloud sides. Part I: Retrieval development and characterization

Abstract. Convective clouds play an essential role for Earth's climate as well as for regional weather events since they have a large influence on the radiation budget and the water cycle. In particular, cloud albedo and the formation of precipitation are influenced by aerosol particles within clouds. In order to improve the understanding of processes from aerosol activation, over cloud droplet growth to changes in cloud radiative properties, remote sensing techniques become more and more important. While passive retrievals for spaceborne observations have become sophisticated and commonplace to infer cloud optical thickness and droplet size from cloud tops, profiles of droplet size have remained largely uncharted territory for passive remote sensing. In principle they could be derived from observations of cloud sides, but faced with with the small-scale structure of cloud sides, classical passive remote sensing techniques are rendered inappropriate. In this work the feasibility is demonstrated to gain new insights into the vertical evolution of cloud droplet effective radius by using reflected solar radiation from cloud sides. Central aspect of this work on its path to a working cloud side retrieval is the analysis of the impact unknown cloud surface geometry has on effective radius retrievals. Using extensive 3D radiative transfer calculations on the basis of realistic droplet size resolving cloud simulations, the sensitivity of reflected solar radiation to cloud droplet size is examined. Sensitivity is enhanced by considering the pixel surrounding to resolve ambiguities caused by illumination and cloud geometry. Based on these findings, a statistical approach is used to provide an effective radius retrieval. An in-depth sensitivity study of the presented approach on the basis of a wide range of radiative transfer test cases demonstrates the feasibility to retrieve cloud particle size profiles from cloud sides.


1 Introduction 1.1 Current state of passive remote sensing of clouds There exist various methods to infer optical properties (e.g. optical thickness and cloud droplet effective radius) from observation of cloud tops from above using information about the scattered and absorbed radiation in the solar spectrum (e.g., Plass 20 and Kattawar (1968); King (1987)). Phase detection is the first step for every cloud property retrieval. Spectral absorption differences in the near-infrared or brightness temperature differences in the thermal infrared are commonly used to distinguish 1 Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2018-234 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 17 September 2018 c Author(s) 2018. CC BY 4.0 License. between liquid water and ice (e.g., Nakajima and King (1990)). Various operational techniques exist to retrieve microphysical cloud properties like cloud thermodynamic phase and effective particle size (e.g., Han et al. (1994); Platnick et al. (2001); Roebeling et al. (2006)).
Remote sensing of cloud and aerosol parameters is mostly done by use of multi-spectral sensors, i.e., using only a limited number of spectral bands. Common examples of spaceborne imagers are the Advanced Very High Resolution Radiometer 5 (AVHRR), the Moderate Resolution Imaging Spectroradiometer (MODIS), and the Spinning Enhanced Visible Infrared Imager (SEVIRI). However, there are concerns about measurement artifacts influencing retrievals of aerosol and cloud properties caused by small-scale cloud inhomogeneity which are unresolved by the coarse spatial resolution of spaceborne platforms (Zinner and Mayer, 2006;Marshak et al., 2006b;Varnai and Marshak, 2007).
Non-imaging systems like the Solar Spectral Flux Radiometer (SSFR, Pilewskie et al., 2003) or the Spectral Modular Air- 10 borne Radiation measurement sysTem (SMART, Wendisch et al., 2001;Wendisch and Mayer, 2003) were used for cloud remote sensing from ground (McBride et al., 2011;Jäkel et al., 2013) or aircraft (Ehrlich et al., 2008;Eichler et al., 2009;Schmidt et al., 2007). The imaging spectrometer of the Munich Aerosol Cloud Scanner (specMACS, Ewald et al., 2016), is the instrument for which the retrieval in this manuscript has been developed. Marshak et al. (2006a); Martins et al. (2011) proposed cloud side scanning measurements and Zinner et al. (2008); Ewald 15 et al. (2013) presented steps towards a cloud side retrieval for profiles of phase and particle size of convective clouds in order to observe the vertical development of cloud microphysics. Similar to earlier satellite retrievals they propose to use solar radiation in the near-visible to near-infrared spectral regions reflected by cloud sides. Especially the vertical dimension of these observations should reflect many aspects of cloud-aerosol-interaction as well as mixing of cloudy and ambient air (Martins et al., 2011;Rosenfeld et al., 2012). In order to provide a resolved vertical profile a fairly high spatial resolution of the order 20 of 100 m or better is needed. This means that a method has to consider 3D radiative transfer effects.
Albeit sophisticated, the studies of Zinner et al. (2008) and Ewald et al. (2013) are limited to an idealized geometry and simplified cloud microphysics. First, they focus on a space-like perspective for a fixed viewing zenith and scattering angle above the cloud field where sun and sensor have the same azimuth. Therefore, their studies lack the varying geometries of an airborne perspective and avoid the challenge to identify suitable observation positions within the cloud field. Moreover, the 25 spatial resolution of their model cloud fields of 250 m is still rather coarse with respect to an airborne perspective of cloud sides.
Second, the effective radius is only parameterized in their studies. For all cloud fields, the effective radius profile is calculated by using a sub-adiabatic ascent of one air parcel in the context of a fixed cloud condensation nuclei concentration. Finally, the approach was not tested for the potential bias to always detect larger effective radii with increasing cloud height; a potential pitfall that could be caused by the prior information contained in the forward calculations. 30 Since the diverse perspectives and the high spatial resolution of airborne cloud side measurements hampered the application of the approach presented by Zinner et al. (2008) and Ewald et al. (2013) until now, the present work will extend and test their ideas in the context of an airborne perspective. In the course of this part 1, following scientific objectives will be addressed: 1. Extend the existing approach to realistic airborne perspectives and develop methods to test the sensitivity of reflected radiances from cloud sides to cloud droplet radius with the observer position within the cloud field. 35 2 Atmos. Meas. Tech. Discuss., https://doi.org /10.5194/amt-2018-234 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 17 September 2018 c Author(s) 2018. CC BY 4.0 License.

Monte Carlo approximation
When no analytical expression for the likelihood probability is available, Monte Carlo sampling from the joint distribution can be used to approximate the likelihood and posterior probability (Mosegaard and Tarantola, 1995). Using a radiative transfer model, this sampling yields a histogram n(L 0.87 , L 2.10 , r eff ) of the frequency of observed radiances L 0.87 and L 2.10 and the corresponding effective radius r eff . With the histogram n as a very simple non-parametric density estimator (Scott et al., 5 1977), the following relation between histogram n and the joined probability p fwd (L 0.87 , L 2.10 , r eff ) and marginal probability p fwd (r eff ) can be made: p fwd (L 0.87 , L 2.10 , r eff ) ∝ n(L 0.87 , L 2.10 , r eff ) (1) p fwd (r eff ) ∝ n(r eff ) = L0.87 L2.10 n. (2) For a successful estimation of these two probabilities, the number of radiative transfer results contained in histogram n needs 10 to be large enough. Simultaneously, the forward simulation has to cover all values expected in the real world application.
With the likelihood probability p(L 0.87 , L 2.10 | r eff ) as a conditional probability, it can be written as the quotient of the joined probability p fwd (L 0.87 , L 2.10 , r eff ) and p fwd (r eff ) describing the ensemble of forward calculations contained in n. In this work, the likelihood probability p(L 0.87 , L 2.10 | r eff ) is approximated by the histogram n(L 0.87 , L 2.10 , r eff ) which is normalized with n(r eff ) using the relations in Equation (1) and Equation (2): 15 p(L 0.87 , L 2.10 | r eff ) = p fwd (L 0.87 , L 2.10 , r eff ) p fwd (r eff ) p(r eff | L 0.87 , L 2.10 ) = p(L 0.87 , L 2.10 | r eff ) p pr (r eff ) p(L 0.87 , L 2.10 | r eff ) p pr (r eff ) dr eff (4) In Equation (3), the prior distribution of r eff in the radiative transfer ensemble is removed by the normalization with the marginal probability p fwd (r eff ). In a final step, the likelihood probability can be used with an arbitrary prior p pr (r eff ) to get 20 the posterior probability p(r eff | L 0.87 , L 2.10 ) given measurements L 0.87 and L 2.10 . Hereby, the arbitrary prior p pr (m i ) must be included within the bounds of the implicit prior p(r eff ) in the forward calculations. Values of r eff that are not included in the forward calculations cannot be retrieved since the likelihood probability p(L 0.87 , L 2.10 | r eff ) is not defined for them. Figure 1 shows how such a Monte Carlo approximation (blue histogram) of a posterior distribution (red line) could look like for given radiance measurements L 0.87 and L 2. 10 . Besides an estimated mean effective radius r eff , the standard deviation of 25 the posterior distribution also yields the uncertainty connected with this estimate.

Radiation transport model
The analysis of radiative transfer effects in one-dimensional clouds is done using DISORT (Stamnes et al., 1988). The representation of 3D radiative transfer in realistic cloud ensembles is done using the Monte Carlo approach with the Monte Carlo code for the physically correct tracing of photons in cloudy atmospheres (MYSTIC; Mayer,2009  with the Monte Carlo sampling of posterior distributions mentioned above, this method will be termed "3D radiative transfer forward modeling" in the following. Both codes are embedded in the radiative transfer library libRadtran (Mayer et al., 2005;Emde et al., 2016) which provides prerequisites and tools needed for the radiative transfer modeling. The atmospheric absorption is described by the representative wavelengths absorption parametrization (REPTRAN; Gasteiger et al., 2014). This parametrization is based on the HITRAN absorption database (Rothman et al., 2005) and provides spectral bands of different 5 resolution (1 cm −1 , 5 cm −1 , and 15 cm −1 ). Calculations have shown that the spectral resolution of 15 cm −1 (e.g. ∆λ = 1.1 nm at 870 nm, ∆λ = 6.6 nm at 2100 nm) best suits the spectral resolution of common hyperspectral imagers. The extraterrestrial solar spectrum is based on data from Kurucz (1994) which is averaged over 1.0 nm. In order to include vertical profiles of gaseous constituents, the standard summer mid-latitude profiles by Anderson et al. (1986) are used throughout this work. Precomputations of the cloud scattering phase function and single scattering albedo are done using the Mie tool MIEV0 from 10 Wiscombe and Warren (1980). When not mentioned otherwise, a Gamma size distribution with α = 7 was used for the Mie calculations. The high computational costs of the 3D Monte Carlo radiative transfer method for tracing large numbers of photons are reduced using the Variance Reduction Optimal Option Method (VROOM) (Buras and Mayer, 2011), a collection of various variance reduction techniques.

15
In order to calculate realistic posterior probability distributions p(r eff | L 0.87 , L 2.10 ), likelihood probabilities, produced by a sophisticated forward model, have to be combined with a realistic prior. While Marshak et al. (2006a) used statistical models to obtain this prior of 3D cloud fields, the physical consistency of cloud structures and cloud microphysics are an advantage of the explicit simulation of cloud dynamics and droplet interactions. Following Zinner et al. (2008), this work applies the threedimensional radiative transfer model MYSTIC to realistic cloud fields which were generated with a large eddy simulation 20 (LES) model on a cloud resolving scale. While Zinner et al. (2008) uses realistic cloud structures combined with a bulk microphysics parametrization, this work extends their approach by including explicit simulations of fully consistent, spectral cloud microphysics. In order to cover clean as well as polluted atmospheric environments, LES model outputs with different background cloud condensation nuclei (CCN) concentrations will be used.  (Rauber et al., 2007). The simulations use an adapted version of the Regional Atmospheric Modeling System (RAMS) coupled to a microphysical model Feingold et al. (1996) and described in more detail in Jiang and Li (2009). In addition to the high spatial resolution, cloud microphysics are explicitly represented by size-resolved simulations of droplet growth within each grid box. The cloud droplet distributions 5 cover radii between 1.56 µm to 2540 µm which are divided into 33 size bins with mass doubling between bins. All warm cloud processes, such as collision-coalescence, sedimentation, and condensation/evaporation are handled by the method of moments developed by Tzivion et al. (1987Tzivion et al. ( , 1989. Droplet activation is included by using the calculated supersaturation field and a given cloud condensation nucleus concentration in two versions with N CCN = 100 cm −1 and N CCN = 1000 cm −1 . The LES simulations (dx25-100 and dx25-1000; Jiang and Li, 2009) have a domain size of 6.4 × 6.4 × 4 km with a spatial resolution of 10 10 m in the vertical and a spatial resolution of 25× 25 m in the horizontal with periodic boundary conditions. As initial forcing, thermodynamic profiles collected during the RICO campaign (Rauber et al., 2007) were used. With condensation starting at a cloud base temperature of around 293 K at 600 m, the cloud depth of the warm cumuli varies over a large range from 40 m to a maximum of 1700 m (Jiang and Li, 2009).
In order to sample a representative prior from this cumulus cloud simulations, a 2 hour (12 h − 14 h LT) model output is 15 sampled every 10 min for both background CCN concentrations. As input for the following radiative transfer calculations, microphysical moments are derived from the simulated cloud droplet spectra. Effective radius r eff , liquid water content LWC and total cloud droplet concentration N d can be calculated from mass mixing ratios m i in g kg −1 and cloud droplet mixing ratios n i in kg −1 given for the 33 LES size bins. With 1067 g m −2 , the LWP maximum is found co-located with a LWC maximum of over 2 g m −3 inside the strongest convective core ( Figure 2a). The growth of cloud droplets with height can clearly be seen in the cloud vertical profile ( Figure 2b).
As the main object of interest of the retrieval developed herein is thie vertical droplet growth, the range of variation of this process in the cloud data sets is explored in more detail. Figure 3 shows so called contoured frequency by altitude diagrams 25 (CFADs; Yuter and Houze, 1995) for effective radius r eff , liquid water content LWC and total cloud droplet number concentration N d . Their typical profile is summarized with black lines over all sampled time steps for N CCN = 1000 cm −1 and red lines for N CCN = 100 cm −1 . The frequently occurring low values of N d and LWC are associated with grid boxes at cloud edges while the wide spectrum of larger values can be found within the cloud cores. As intended the two cloud ensembles cover a wide range of possible values for r eff and N d between low ("clean") and high CCN concentration ("polluted"). Small droplet 30 r eff and slow growth during ascent characterize the N CCN = 1000 cm −1 case, while quick growth to much higher values is present for smaller N CCN = 100 cm −1 . LWC and cloud lower and upper boundaries show only small differences.
Based on the cloud base droplet number N cb , temperature T cb , pressure and a saturation adiabatic lapse rate (here we assume 4 Kkm), "adiabatic" reference values can be calculated for an ensemble of droplets growing by condensation during ascent neglecting entertainment of dry environmental air (broken lines). The existence of other effects (e.g. entrainment, coalecence) becomes evident in comparison with the modeled LWC and N d profiles, as the adiabatic theory provides only an upper limit to their values. In contrast, r eff follows the adiabatic limit more closely with sub-adiabatic values between 60 − 80% which is in accordance with in-situ aircraft observations during the RICO campaign (Arabas et al., 2009) and other studies (Martin et al., 1994).

Selection of suitable cloud sides
A key component of the Bayesian approach is the selection of a suitable sampling strategy to explore the likelihood distribution p fwd (L 0.87 , L 2.10 | r eff ). This is especially true, if the sampling of the observation parameter space is done with a computational expensive 3D radiative transfer method. Following Mosegaard and Tarantola (1995), the sampling of the model space can be 5 improved when the model space is sampled with the intended measurements in mind. Instead of sampling the radiative transfer in 3D cloud fields completely at random, the indented measurement location and perspective should be taken into account.
To that end, the we introduce a technique to select suitable locations within the LES model output for which cloud sides are visible from the airborne perspective. Cloud side measurements are intended for clouds within several kilometers from the instrument location. With the sun in the back, azimuthal positions of ±45°around the principal plane will be accepted 10 for an airborne field of view which is centered slightly below the horizon. To ensure reproducibility, an analytical method is chosen to select observation locations to sample the likelihood distribution. The field-of-view is modeled by an observation kernel k FOV with an azimuthal opening angle of ∆ϕ = 45°and a zenithal opening angle of ∆ϑ = 40°, centered around 5°b elow the horizon. As a function of radial distance, the observation kernel comprises a scalar weighting to curtail the location where clouds are desired. Figure 4 shows a vertical cross-section of this kernel. The arbitrary score is strongly negative in the 15 vicinity of the observer to penalize locations where clouds are too close. Observation distances of 3 km to 5 km turned out to maximize the likelihood to observe a complete cloud side in the used LES model output. For a distance of 2 km and onward, the weighting score becomes thus positive with a maximum at 3.5 km to favor locations with clouds in this region. For all LES cloud fields on average, this method positions the observer at a distance of around 4 km from cloud sides.
Subsequently, the field of cloudy grid boxes is convolved with the observation kernel at an observation altitude of h = 1.7 km, creating a two-dimensional score field s obs . For every cloud field and chosen azimuthal orientation, the observation position is then placed where s obs has its global maximum. In Figure 4a, the already introduced LES cloud field (12 h 40 min LT) is 5 shown in combination with the corresponding score field s obs obtained for a viewing azimuth of φ = 315°. The observation position is indicated by the yellow dot, where s obs has its global maximum as recognizable by the green color. Also depicted is the field of view towards the largest cloud in the center of the domain. The red region in s obs would be unfavorable for a cloud side perspective since it would be too close to the cloud. For the selected perspective shown in Figure 4a, a simulated truecolor image is shown in Figure 5a.  method has its limitations when it comes to highly structured cloud sides with horizontally inhomogeneous microphysics. By neglecting the penetration depth of photons, reflection from deeper within the cloud are disregarded. Radiance observations at given wavelengths can not carry the information from this first grid box alone, but from the full multi-scattering path. The found effective radius r eff therefore becomes biased towards droplet sizes found directly at cloud edges. However, due to very low LWCs, these grid boxes have only a marginal contribution to the overall reflectance. 5 As Platnick (2000) showed, the penetration depth of reflected photons in the visible spectrum lies within some hundred meters while in the near-infrared spectrum the penetration depth is only a few dozen meters. The co-registration of responsible cloud droplet sizes with modeled radiances is essential. Besides the observation perspective, this apparent effective radius r eff app also depends on the observed wavelength since different scattering and absorption coefficients lead to different cloud penetration depths. 10 As discussed by Platnick (2000), there exist analytical as well as statistical methods to consider the contribution of each cloud layer to the apparent effective radius r eff app . Advancing the one-dimensional weighting procedures of Platnick (2000) and Yang et al. (2003), the 3D tracing of photons in MYSTIC is utilized to calculate the optical properties of inhomogeneous, mixed-phase clouds. The apparent effective radius r eff app for a photon can be described as a weighted, linear combination (Equation (5)) of the individual effective radii r eff the photon encounters on its path through the cloud: In Equation (5), the effective radii are weighted with the corresponding extinction coefficient k ext of the cloud droplets along the path length in each grid box. Subsequently, the mean over all photons traced for one forward simulated pixel leads to the apparent effective radius r eff mc of this pixel: In the summation Equation (6), the photon weight p w,n is used to account for the different photon path probabilities. As the photon weights p w,n are also used in the calculation of L 0.87 and L 2.10 , the apparent effective radius r eff mc can be derived simultaneously. This method was integrated within the MYSTIC 3D code and will therefore be referred to as the MYStic method To Infer the Cloud droplet EFFective Radius (MYSTIC REFF). For the cloud scene shown in Figure 5a, Figure 5b shows the apparent effective radius r eff mc obtained with MYSTIC REFF. Compared with the effective radius found at the 5 cloud edge, r eff mc appears much smoother. The vertical gradient of r eff mc also compares better with the vertical gradient of effective radii r eff shown in Figure 3b. The method shows very good agreement with the analytical solution of Yang et al. (2003) for homogeneous mixed-phase clouds and Platnick (2000) for one-dimensional clouds with a vertical effective radius profile.
4 The cloud geometry effect and its mitigation 10 Reflected radiance at non-absorbing wavelengths is mainly influenced by optical thickness and by the amount of radiation incident on the cloud surface. For the latter, cloud surface orientation with respect to the sun is decisive. This is a problem for all retrievals using this radiance to derive τ c or, in combination with an absorbing wavelength, τ c and r eff (Nakajima and King, 1990) in all situations with unknown (non-plane-parallel) clouds. In such a situation, e.g. for cloud side observations, the limitation to optically thick clouds can be a way out. With increasing optical thickness τ c , the solar cloud reflectance 15 becomes less sensitive to variations of τ c . In contrast to the typical observation geometry from above where a plane-parallel cloud is assumed, cloud surface orientation is unknown and only the scattering angle ϑ s is known. In the following study, the ambiguity caused by the unknown cloud surface orientation and the remaining sensitivity to the effective radius will be explored. Molecular absorption and scattering will be neglected in the following idealized study.

Ambiguities of reflected radiances
20 Figure 6 shows the basic geometry for cloud side remote sensing. The normal which defines the cloud surface isn, the vector pointing from sun into direction of light propagation is denoted withŝ and the vector pointing into direction of the observer is denoted withv. The viewing zenith angle ϑ and the sun zenith angle ϑ 0 are still referenced within the global coordinate system which is defined parallel to the ground. Corresponding to these two angles, two additional angles exist which describe the inclination ofŝ andv on the oriented cloud surface: the local illumination angle ϑ * 0 and the local viewing angle ϑ * with 25 respect to the cloud surface. In a first step, all vectors are assumed to be within the principal plane (the plane spanned byŝ and n). Figure 6 shows two different viewing geometries of a vertically-oriented cloud surface under the same local illumination angle ϑ * 0 . For both cases, Figure 7 shows spectral radiances at λ = 870 nm and λ = 2100 nm for a clockwise rotation of the cloud surface. The radiative transfer calculations were done with DISORT for an optically thick (τ = 500) water cloud with a fixed r eff of 9 µm. The arrows in Figure 7 indicate the progression of radiance values, as it could be observed during a 30 cloud surface rotation within the principal plane as shown in Figure 6. The figure uses a typical two-channel diagram with the absorbing channel on the x-axis and the non-absorbing on the y-axis. E.g. Nakajima and King (1990) Figure 7 the dependence of reflected radiance in both channels on the systematic variation of τ and r eff values for plane-parallel clouds (hereafter denoted as "2-wavelength retrieval" and "2-wavelength diagram"). The similarity of these lines to the isolines for fixed r eff and varying τ in their diagrams is striking. Nakajima and King (1990) use the radiance in the non-absorbing wavelength to infer the optical thickness τ . Numerous studies Cahalan et al. (1994); Varnai and Marshak (2002); Zinner and Mayer (2006); Vant-Hull et al. (2007) pointed out, that 5 tilted and therefore more shadowed or illuminated cloud sides have a huge impact on the retrieval of optical thickness.The radiance similarity of cloud surface rotation and optical thickness variation further underlines the necessity to restrict the retrieval to optically thick clouds when the cloud surface orientation is unknown. Therefore, the following analysis will determine the remaining information content in L 0.87 and L 2.10 about r eff for optically thick clouds with an unknown cloud surface orientation. 10 A difference between the direct backscatter case with a scattering angle of ϑ s = 180°and the case with a scattering angle of ϑ s = 150°becomes already evident in Figure 7. While the spectral radiance first increases at both wavelengths as the local illumination and viewing angle becomes smaller, it is only in case of direct backscatter that spectral radiances decrease the same way as they increased when the illumination angle becomes more oblique again. In case of scattering angle ϑ s = 150°, as long as the local viewing angle ϑ * is smaller than the local illumination angle ϑ * 0 , spectral radiances at λ = 2100 nm are lower 15 than for the remaining part of the rotation when ϑ * > ϑ * 0 . Within the principal plane, a rotation of the cloud surface produces a characteristic bow structure for scattering angles ϑ s < 180°(evident in Figure 7). For optically thick clouds with unknown cloud surface orientation, this introduces the ambiguity to be dealt with into the relation between L 0.87 , L 2.10 and cloud optical properties. For the oblique viewing geometry, radiances for r eff = 9 µm within the upper branch of the bow structure coincide with radiance values for r eff = 7 µm. Nevertheless, towards higher values of L 0.87 there remain unambiguous regions where  Figure 7. Calculations of spectral reflection were done for an optically thick (τ = 500) water cloud with a fixed effective radius r eff = 9 µm with a fixed scattering angle of ϑs = 180°(orange line, 9 µm) and three different effective radii with a fixed scattering angle of ϑs = 150°(blue lines, 7 µm, 9 µm and 13 µm).

3D
Next, the analysis is extended from principal plane considerations to a full 3D setup. To this end, 3D MYSTIC radiance simulations were done for a spherical, optically thick cloud and different scattering regimes of ϑ s = 180°and 150°. Figure Figure 8  shows the results in 2-wavelength diagrams (center) and as images for the direct backscatter direction (left) and for a scattering angle of 150°(right). In the 2-wavelength diagrams the radiance pairs (L 0.87 , L 2.10 ) from the 3D radiative transfer forward simulation are shown as scattered points, the results from the one-dimensional DISORT simulations for different effective radii are shown as isolines. Within the principal plane, the large green and red dots indicate cloud surfaces with same local illumination angle ϑ * 0 = 30°. Figure 9 illustrates the viewing and illumination setup at the green and red dot. In the following, the 5 green dot will mark the cloud surface with smaller local viewing angle (ϑ * <ϑ * 0 ), the red dot the cloud surface with larger local viewing angle (ϑ * >ϑ * 0 ). In the 2-wavelength diagrams, the gray isolines mark radiance pairs with same local illumination angle ϑ * 0 . The corresponding images show radiances L 0.87 , L 2.10 and the colored ratio L2.10 L0.87 to aid mapping radiance pairs with their spatial location.
This ratio reflects the radial symmetry of the local illumination angle for ϑ s = 180° (Figure 8 left part). For this direct 10 backscatter geometry in Figure 8, 3D results for r eff = 9 µm match the 1D DISORT results for r eff = 9 µm very closely. Due to the radial symmetry of the local illumination angles, radiance values also decrease radially symmetric with more oblique cloud surfaces. Albeit restricted to airborne or spaceborne platforms, this perspective minimizes the 3D effect on radiance ambiguities caused by unknown cloud surface orientations. The picture changes when the observer leaves the backscatter geometry as shown for a scattering angle of ϑ s = 150°in Figure 8 on the right. As already shown with the DISORT results in 15 Figure 7, the radiance pairs form a bow-like pattern with higher L 2.10 values at more oblique surface orientations. Furthermore, the red and green dots with same local illumination angle now become separated since ϑ * = ϑ * 0 . While radiance at the nonabsorbing wavelength drops considerably with a more oblique local viewing angle (ϑ * = 60°, red dot), radiance at the absorbing wavelength changes only slightly at very steep local viewing angle (ϑ * = 0°, green dot). Consequently, droplets at the green dot with r eff = 9 µm could be miss-interpreted as effective radius r eff = 7 µm or even 5 µm . 20 For a deeper insight into the different observation at the green and red dot, we analyze the angular distribution of cloud reflectance at the absorbing and non-absorbing wavelength. For a fixed illumination angle ϑ * 0 = 30°, Figure 10 shows the DISORT radiances for different effective radii and all possible local viewing angles. The two locations considered in Figure 8  (right) are, again, marked with green and red dots. Obviously, the angular characteristic differs between the absorbing and non-absorbing wavelength. For a wide range of scattering angles ϑ s < 180°, the radiance at the absorbing wavelength remains quite constant between steep (green dot) and oblique (red dot) viewing perspective while the radiance at the non-absorbing wavelength is considerably smaller for the oblique viewing perspective. This asymmetric behavior becomes less pronounced for scattering angles near ϑ s = 180°, when the red and green dots both move symmetrically towards the backscatter peak in 5 Figure 10. The reason for this different angular reflectance is connected with different photon penetration depths at the two wavelengths. Smaller penetration depths lead to a more uniform reflection, while larger penetration depths lead to a stronger reflection perpendicular to the cloud surface. Marshak et al. (2006a) and Zinner et al. (2008) do not investigate these reasons for the cloud geometry caused scatter of reflected radiance in detail. Nonetheless, they suggest to limit the influence of missing geometry information by additional 10 consideration of vertical thermal radiation temperature gradients (containing part of the geometry information). In the following a more systematic use of available geometry information in the visible, near-infrared spectrum itself is presented.

Additional information from surrounding pixels
Here a technique is presented that uses information from surrounding pixels to classify the geometrical environment of the considered pixel. Already Varnai and Marshak (2003)  As discussed in the preceding section, ambiguous radiances are mainly caused by the lower L 0.87 radiances from oblique cloud surfaces when compared to the more steeper cloud surfaces. With the sun behind the observer and convex cloud structures, cloud regions with steeper surfaces should generally be brighter than their surroundings while oblique surfaces should be darker than their surroundings. In line with recent developments of multi-pixel retrievals, the method should therefore determine if the pixel is surrounded by darker pixels or surrounded by brighter pixels. At the same time, the method should be robust 5 regarding instrument or radiative Monte Carlo noise between adjacent pixels. To this end, a 2D Gaussian band-pass filter is used to classify different illumination regimes in simulated or measured radiance images. As a 2D band-pass filter, it compares the brightness of each pixel with the brightness of other pixels in the periphery. Pixels are classified according to their positive or negative radiance deviation compared to their surrounding pixels. The band-pass filter consists of two 2D Gaussian functions H LP (x, y) and H HP (x, y) which specify the inner and outer search radius for this comparison. By neglecting directly adjacent pixels, the filter is insensitive to pixel-to-pixel noise of CMOS sensors or Monte Carlo radiative transfer calculations. With D(x, y) as the distance (Equation (7)) from the origin pixel, D H and D L limit the inner and outer search radius respectively: The absolute deviation values obtained from H BP (x, y) can vary from scene to scene. For the classification into steep or oblique perspectives, we are more interested in the relative deviation. To constrain the gradient classifier g class (x, y) into a fixed interval, 20 the arcus tangent function is used (Equation (11)).
To demonstrate the method, the cloud field illustrated in Figure 5 was used again for radiance calculations, but with a fixed effective radius of r eff = 8 µm. For this fixed effective radius, the broad radiance distribution of L 0.87 and L 2.10 shown in Figure 11a is mainly caused by the different cloud surface orientations discussed in the previous section. In order to identify the regions leading to the upper part of the radiance scatter cloud in Figure 11a, an exponential function was fitted (black line) 25 to the data points to determine the positive (blue) or negative (red) deviation ∆L 2.10 from the typical best fit (black) line for each radiance pair.
In the following, this deviation ∆L 2.10 is taken as a reference for a perfect separation of 3D radiance ambiguities. A method that would yield the same separation from observable parameters could be used as a proxy to mitigate the problem of ambiguous radiances. For the gradient classifier g class , a maximal correlation with ∆L 2.10 is found when the band-pass operates between   Figure 11b shows the separation provided by the gradient classifier (Equation (11)). Reference and proxy are also shown as images in Figure 12, where the deviation ∆L 2.10 from the fit in the scatter plot is shown as reference on the left (Figure 12a) and the result of the gradient classifier on the right ( Figure 12b).
Apparently, the gradient classifier g class is able to reproduce the general appearance of the gradient classifier g class used in 5 Figure 12. It is able to separate the radiance distribution into positive and negative radiance deviations ∆L 2.10 at high as well as at low radiance values. Large g class values are more likely to be associated with a steeper illumination angle compared to the viewing angle, while smaller values are more likely to be associated with a more oblique illumination angle compared to the viewing angle. For two pixels with same illumination angle, g class > 0 therefore marks the upper radiance branch in Figure 7 for ϑ s = 150°, while g class < 0 marks the lower radiance branch. Based on this feature, the band-pass classification g class can 10 17  be used as a proxy to determine the location of a pixel within the radiance distribution. In the following, the gradient classifier g class is used for the 3D forward calculation ensemble as well as for real measurements.

Exclusion of cloud shadows
Depending on the illumination, cloud surfaces can also be in direct shadow if the local solar zenith angle onto the cloud surface ϑ ′ 0 is larger than 90°. Illuminated cloud parts can also cast shadows onto other cloud parts. Without direct illumination, reflected 5 photons from these parts originate from previous scattering events and are affected by those. For this reason, shadowed cloud parts have to be filtered out before applying any retrieval based on direct illumination.
Usually radiation from shadow regions encountered more absorption compared to directly reflected light (Vant-Hull et al., 2007). This enhanced absorption is visible in Figure 13, where the reflectivity at 2.1 µm drops considerable ( Figure 13a). As a proxy of enhanced absorption, the reflectivity ratio R 0.87 /R 2.10 ( Figure 13b) increases in this region. In the following, this ratio 10 will be used as shadow index f shad to exclude pixel for which light has likely undergone multiple diffuse reflections: The manual inspection of many scenes confirmed f shad > 3.5 as viable threshold for most shadowed cloud regions. In addition, a threshold of L 0.87 > 75 mWm −2 nm −1 sr −1 is used to focus the retrieval on optically thicker clouds and to filter out clearsky regions.

Retrieval
In this section, the Monte Carlo sampled posterior distributions p(r eff | L 0.87 , L 2.10 ) will be used to infer droplet size profiles from convective cloud sides. As mentioned in Section 2.2, the posterior p(r eff | L 0.87 , L 2.10 ) can be derived from Bayes' theorem by solving the easier forward problem p(L 0.87 , L 2.10 |r eff ) for all values of r eff .  The three-dimensional radiative transfer code MYSTIC is applied to LES model clouds to obtain simulations of realistic specMACS measurements. A whole ensemble of these MYSTIC forward simulations of cloud sides will then be incorporated within the statistical framework introduced in Section 2.1. Subsequently, the sampled statistics of reflected radiances are analyzed for their sensitivity to the effective cloud droplet radius.
5.1 Implementation of the 3D forward radiative transfer ensemble 5 In the following, an ensemble of 3D radiative transfer simulations is created to sample the posterior probability distribution p(r eff | L 0.87 , L 2.10 ). The ensemble of simulated cloud side measurements is set up by using the method to select suitable observation perspectives introduced in Section 3.1. During the radiative transfer calculations, the MYSTIC REFF method (Section 3.1) determines the apparent effective radius which links the simulated radiances with the corresponding cloud droplet sizes. Despite the variance reduction methods in the MYSTIC code itself (Buras and Mayer, 2011), the time-consuming 3D 10 technique still limits the number of model runs. Figure 14 illustrates the different illumination setups and the viewing geometry included within the 3D forward simulation ensemble. With LES cloud tops between 1.5 and 2.0 km, the airborne perspective is set to an altitude of h = 1.7 km. Since the retrieval should also be applicable in tropical regions, solar zenith angles were chosen at ϑ 0 = 7, 27, 47 and 67°. impact on radiances at L 0.87 and L 2.1 (analysis not shown), the scattering properties were derived according to Mie theory using modified gamma size distributions with a fixed width of α = 7 and the effective radius as simulated by RAMS. For the ensemble, the surface albedo was set to zero since the influence of radiation reflected by vegetation on the ground is masked 5 out in measurements. This technique will be described in the following part 2 paper.
Atmospheric aerosol was included by using the continental average mixture from the Optical Properties of Aerosols and Clouds (OPAC) package (Hess et al., 1998). The aerosol optical thickness (AOT) at 550 nm is around τ 550 a = 0.15 for this profile. This aerosol profile is typical for anthropogenically influenced continental areas and contains soot and an increased amount of insoluble (e.g. soil) as well as water-soluble (e.g. sulfates, nitrates and organic) components.  20 In the next step, simulated radiances were binned into a multidimensional histogram with equidistant steps in L 0.87 , L 2.10 , r eff , ϑ and g class . Table 1 shows the range and step sizes of this histogram. During this discretization, radiances were counted in adjoining bins by linear interpolation. In the following, the histogram will be normalized to yield the posterior probability p(r eff | L 0.87 , L 2.10 ).

Biased and unbiased priors 25
In our input cloud field, the effective radius always increases with height which might impact the retrieval. The retrieval should not exhibit any trend towards a specific profile. Instead of an unbiased result, the retrieval would otherwise only reflect a-priori knowledge about the vertical profile of cloud microphysics. For this reason, the assumed prior is a key element to be considered in the sampling of the posterior and the subsequent Bayesian inference. In the case of cloud side remote sensing, two possible priors p pr (r eff ) come into mind: a uniform prior or the LES model provided prior. For the LES model the prior is a function of 30 viewing geometry, as some r eff are more likely to be observed under certain viewing directions. In particular, relative frequency of r eff for different scattering angles ϑ s and gradient classes should be the same. For this problem, the prior probability should be uniform in r eff to avoid the introduction of any bias.
Another important prerequisite of the Monte Carlo based Bayesian approach is the sufficient sampling of the likelihood probability p(L 0.87 , L 2.10 , r eff ). Naturally, effective radii not included in the ensemble of forward calculations cannot be retrieved using Bayesian inference. Furthermore, it should be kept in mind that sparsely sampled likelihood regions are probably not representative for the whole distribution. This is especially true for the smallest and largest effective radii contained in the LES model.
These values are chosen to be at least 4 µm larger than the maximum values found in the respective ensembles to ensure positive and realistic values for r flip eff . In order to preserve the optical thickness τ orig of the original cloud field (Equation (16)), the original liquid water content LWC orig is changed (Equation (17)

Radiance and posterior distributions
The following section will present the radiance histograms n(L 0.87 , L 2.10 , r eff ) and the corresponding posterior distributions p(r eff | L 0.87 , L 2.10 ). Analogous to the likelihood distribution, the first gives the spread of radiances for a given effective radius r eff , while the latter describes the spread of effective radii r eff for a given radiance pair L 0.87 and L 2. 10 . Figure 15 shows a 2D histogram of simulated radiance combinations L 0.87 and L 2.10 for the airborne perspective. The histogram shows the results for 5 the effective radius bin centered at r eff = 10 µm, the scattering angle bin between ϑ s = 130°and 140°and the gradient class bin g class = 4 which holds pixels that are brighter as their surroundings. The spread of the radiance from three-dimensional cloud sides, for the most part, can be explained by the one-dimensional DISORT results (dashed line for variable cloud surface inclination).
After normalization of the histograms in Equation (3) and after the application of the uniform prior in Equation (4), the 10 posterior probabilities p(r eff | L 0.87 , L 2.10 ) can be examined. Figure 16 shows posterior probabilities as a function of r eff for different radiances L 2.10 at the absorbing wavelength corresponding to the colored dots in the histogram panels ( Figure 15).
The vertical lines indicate the corresponding mean effective radius for each posterior distribution which were derived using Equation (18). The descending order of mean effective radii with ascending radiance L 2.10 demonstrates the general feasibility to discriminate different effective radii in cloud side measurements. Albeit the relatively large spread in r eff , the measurement of a radiance pair (L 0.87 , L 2.10 ) can still narrow down r eff to ±1.5 µm around the most likely value. In the following, the tabulated set of posterior distributions is used as lookup table for the effective radius retrieval.

Bayesian inference of the effective radius
Based on this lookup table of posterior probabilities p(r eff | L 0.87 , L 2.10 ), the actual retrieval of effective radii can now be introduced. After a set of spectral radiance pairs L 0.87 and L 2.10 has been measured, the band-pass filter (Section 4.2) is 5 applied to the L 2.10 image to derive the gradient classifier g class . Scattering angles are calculated from the orientation and navigation data of the aircraft. With the four parameters, L 0.87 , L 2.10 , g class and ϑ s defined for each pixel, the corresponding posterior is retrieved from the lookup table by linear interpolation between posteriors defined at the bin centers of the lookup table. Finally, the mean effective radius r eff and the corresponding standard deviation σ(r eff ) can be derived as first and second moments of the posterior distribution: This 1-sigma standard deviation σ(r eff ) will be referred to as the statistical retrieval uncertainty.
6 Numerical analysis of the retrieval The next section will examine the stability of the statistical relationship between reflected radiance and cloud droplet size.

15
How well can we retrieve the cloud droplet size after different viewing directions and cloud surface orientations have been combined within one lookup table? To answer this, the statistical retrieval is applied to simulated cloud side measurements for which the underlying effective radius is known. First, this is done for scenes that have already been included in the lookup table. Using scenes with normal and flipped effective radius profile, the lookup table is tested for an inherent bias towards a specific effective radius profile that could be caused by the chosen forward sampling strategy. In addition, tests are repeated for the same scenes with a fixed effective radius of 8 µm; a case which is not included in the lookup table.
6.1 Analysis of the sampling bias 5 By design the retrieval should not exhibit any trend towards a specific profile. The polluted scene (N CCN = 1000 cm −1 ), already introduced in Figure 5a, will be used as a first case study. Figure 17 shows result of the statistical effective radius retrieval with the normal effective radius profile on top, the flipped profile in the center and a fixed effective radius profile at the bottom,. The retrieved mean effective radius is shown in the right panels (Figure 17c,f,i), the apparent effective radius r eff mc ("the truth") is shown in the left panels (Figure 17a,d,g). The center panels compare the mean vertical profile (lines) and its 10 spatial standard deviation (shaded areas) of the apparent (black) and the retrieved (green) effective radius. Furthermore, the green error bars show the error estimate σ(r eff ) provided by the retrieval.
Overall, the retrieval reproduces all three profiles quite well. However, there are also large differences (up to ±3 µm) for specific cloud regions. For the normal profile as well as for the flipped profile, the retrieved effective radius agrees well with its true value for the upper half of the cloud side. In the lower half of the cloud side near cloud base, the retrieval underestimates 15 the large droplets of the flipped profile by up to 1 µm.
Altogether, the statistical relationship between reflected radiance and cloud droplet size seems stable enough to be used for highly complex cloud sides. Moreover, these first results indicate that the retrieval seems to be resilient to the unrealistic, flipped cloud profiles included in its lookup table. Although the retrieval showed minor problems to retrieve the flipped profile, no substantial bias towards a specific effective radius profile could be detected. Mean values for all heights agree within the 20 natural variability in the LES data; the retrieval error estimate σ(r eff ) seems to overestimate the uncertainty. Since these results were only obtained for a single cloud side scene, the following section will investigate these findings for a representative number of scenes.

Statistic stability for included scenes
In a first step, the retrieval will be tested for perspectives which are already included in the lookup table. This is done to test the 25 retrieval for biases and to obtain a robust measure of correlation between the retrieval and the cloud side scenes it is composed of. By comparing this correlation with the correlation for cloud side scenes that are not included in the lookup table, this analysis will also be used to detect a potential over-fitting. There is the risk that the lookup table only reflects 3D effects that are specific for the included cloud side scenes.
In total, 9 cloud side perspectives were randomly chosen from the polluted as well as the clean dataset. For each perspective, 30 the normal, the flipped as well as the fixed effective radius profile were tested. This amounts to 2 × 9 × 3 = 54 test cases. For this statistical comparison, only reliable results with an retrieval error estimate σ(r eff ) of less than 2.5 µm were included.    Counts Figure 19. 2D histograms and linear regressions to determine the correlation between the true effective radius r eff and the retrieved effective radius reff,retr for (a) the included polluted cases (NCCN = 100 cm −1 ) with normal and flipped effective radius profile and (b) for the not included cases with normal and flipped effective radius profile.
slightly better ability to detect the cases with a normal effective radius profile. Nevertheless, comparable linear regressions show no substantial bias towards the normal or the flipped cases. This confirms the observation made in the case study shown in Figure 17.
For the fixed effective radius profiles, the histogram in Figure 20 shows the deviation of the statistical retrieval with two distinct modes. With most likely values between 0 µm and −1 µm, the retrieval underestimates the effective radius slightly.

5
A second mode is found where the retrieval overestimates r eff with values between 1 µm and 2 µm. In combination, there is  only a slight overestimation of 0.10 µm with a larger standard deviation of 1.17 µm. A further investigation showed that the overestimation peak is connected with and found around undetected cloud shadows.

Statistic stability for unknown scenes
To check the retrieval for potential over-fitting, the retrieval was applied to unknown cloud side scenes. Nine new cloud side perspectives were selected from the polluted and the clean LES runs. While the forward ensemble contains viewing azimuths 5 of 45°, 135°, 225°and 315°, these new cloud side perspectives were chosen for new viewing azimuths of 0°, 90°, 180°and 270°and only normal effective radius profiles were used this time. Figure 18 shows one of these new cloud side scenes with a normal effective radius profile. In contrast to Figure 17, this clean (N CCN = 100 cm −1 ) scene features a much larger range of cloud droplet sizes. Again, the statistical retrieval detects the effective radius profile well. The true values for r eff remain within the retrieval error estimate σ(r eff ).The comparison for all 10 not included cloud sides is shown in Figure 19b. The correlation for the not included cases is 0.93, where around 339, 000 pixel are compared in total.
With nearly the same correlation coefficient, the retrieval performance remains the same when faced with unknown cloud side scenes. It can therefore be concluded that the retrieval is not trained only for the included cloud side scenes. Rather, it represents the statistical relationship between reflected radiance and cloud droplet size for this cloud ensemble.  size was analyzed in the context of a unknown cloud geometry.
(2) In this course, 3D radiative effects caused by the unknown cloud surface orientation were investigated. A technique was proposed to mitigate their impact on cloud droplet size retrievals by putting pixel in context with their surrounding. (3) Finally, an effective radius retrieval for the cloud side perspective was developed and tested for cloud scenes which were used during the retrieval development as well as for unknown scenes. The scope of this work was limited to the liquid part of convective liquid water clouds, e.g. Cumulus and Cumulus congestus, which 5 exhibit well-developed cloud sides. In principle, the proposed technique could also be extended to ice clouds.
In a first step, this work introduced a statistical framework for the proposed remote sensing of cloud sides following Marshak et al. (2006a). A statistical relationship between reflected sunlight in a near-visible and near-infrared wavelength and droplet size is found following the classical approach by (Nakajima and King, 1990). By simulating the three-dimensional radiative transfer for high-resolution LES model clouds using the 3D radiative transfer model MYSTIC (Monte Carlo code for the 10 physically correct tracing of photons in cloudy atmospheres), probability distributions for this relationship were sampled.
These distributions describe the probability to find a specific droplet size after a specific solar reflectance pair of values has been measured. In contrast to many other effective radius retrievals, this work thereby provides essential information about the retrieval uncertainties which are intrinsically linked with the reflectance ambiguities caused by three-dimensional radiative effects. Furthermore, this work developed a technique (Section 4.2) to reduce 3D radiance ambiguities when no information 15 about the cloud surface orientation is available. More precisely, additional information from surrounding pixels was used to classify the environment of the considered pixel. The numerical analysis of the statistical retrieval showed a RMSE between retrieved and true r eff of around 1 µm to 1.5 µm.
For the airborne measurement perspective aimed for (see part 2 of this work), the statistical retrieval reliably detects the present effective radius profile, while sanity checks showed no prior bias of the retrieval towards specific cloud droplet size profiles.
This is an essential prerequisite for all consecutive interpretation of the retrieval results. Furthermore, the retrieval performance remained the same when faced with unknown cloud side scenes not included in the ensemble used for the retrieval. It can 5 therefore be concluded that the retrieval is not over-fitted and that it represents the statistical relationship between reflected radiance and cloud droplet size for this cloud side perspective.
Limited to optically thick water clouds, this work investigated the main reason for reflectance ambiguities from cloud sides.
Within the same 3D cloud side scene, viewing perspectives onto cloud surfaces can be steeper or more oblique as the illumination angle. As a consequence, the correlation between reflected solar radiance pairs and droplet sizes becomes ambiguous. 10 The next important step is the application of the proposed retrieval technique to real measurements. In combination with simultaneous in-situ measurements, airborne cloud side measurements have been acquired with the hyperspectral cloud and sky imager specMACS. In a follow-up paper part 2, the proposed retrieval will be validated with this independent in-situ data.
A further important point is the development of a distance mapping for the retrieval. The height and location assignment of retrieval results is not just of uttermost importance for the comparison with in situ measurements and models, but also essential 15 to estimate the cloud distance for a potential aerosol correction. Here, first promising results could be achieved by exploiting the oxygen A-band absorption at λ = 762 nm presented in Zinner et al. (2018). In conclusion, the present work developed a working effective radius retrieval to measurements of clouds sides applicable to real measurements and thus paved the way for further research on this topic.