Using Two-Stream Theory to Capture Fluctuations of Satellite-Perceived TOA SW Radiances Reflected from Clouds over Ocean

Shortwave (SW) fluxes estimated from broadband radiometry rely on empirically gathered and hemispherically resolved fields of outgoing top-of-atmosphere (TOA) radiances. This study aims to provide more accurate and precise fields of TOA SW radiances reflected from clouds over ocean by introducing a novel semi-physical model predicting radiances per narrow sun-observer geometry. Like previous approaches, this model was trained using CERES-measured radiances paired with MODIS-retrieved cloud parameters as well as reanalysis-based geophysical parameters. By using radiative transfer approxi5 mations as a framework to ingest above parameters, the new approach incorporates cloud-top effective radius and above-cloud water vapor in addition to traditionally used cloud optical depth, cloud fraction, cloud phase, and surface wind speed. A twostream cloud albedo—serving as a function of cloud optical thickness and cloud-top effective radius—and Cox-Munk ocean reflectance were used to describe an albedo over each CERES footprint. A simple equation of radiative transfer, with this albedo and attenuating above-cloud water vapor as inputs, was used in its log-linear form to allow for statistical optimization. 10 We identified the two-stream cloud albedo that minimized radiance residuals and outperformed the state-of-the-art approach for most observer-geometries and solar zenith angles between 20◦ and 70◦, reducing median standard deviations of radiance residuals per solar geometry by up to 13.2% for liquid clouds, 1.9% for ice clouds, and 35.8% for footprints containing both cloud phases. Tested for a variety of scenes, we further demonstrated the plausibility of scene-wise predicted radiance fields. This new approach may prove useful when employed in Angular Distribution Models and may result in improved flux estimates, in 15 particular dealing with clouds characterized by small or large droplet/crystal sizes.

T +243.5 (using T in degree Celsius) (Bolton, 1980) and molecular weights of water and dry air mol h2o and mol air respectively: mr(p, T, rh)dp = 1 g 0 p ctop e s (T ) mol h20 mol air rh(p)dp (1) 90 For our analysis, we filtered the extracted dataset for samples with more than 95% water surface, more than 0.1% cloud fraction, and solar zenith angles between 0 • and 82 • . Table 1 lists the resulting subset of 1.7 billion samples.

Methods for Capturing Radiance Fluctuations
In order to provide hemispherically resolved fields of backscattered radiances to radiance-to-flux-converting ADMs, statistical approaches capture observed radiances together with prevalent scene properties per narrow and discretized Sun-observer-95 geometry. Following Su et al. (2015), solar zenith, viewing zenith, relative azimuth angles were discretized into 2 • intervals, referred to as θ ∆ 0 , θ ∆ v , and ϕ ∆ , respectively. Combinations of θ ∆ 0 , θ ∆ v , and ϕ ∆ were denoted as angular bins and observations were sorted into bins for separate treatment. The following subsection presents the state-of-the-art methodology. Subsection 3.2 introduces a novel semi-physical approach that includes additional parameters. proposed log-linear model (b) capture fluctuations of CERES-measured TOA SW radiances. As this angular bin is within Sun-glint region, (a) shows the LUT-approach for x < 6 (as defined in Section 3.1; note that f1 and f2 are taken between 0-100 in (a) and between 0-1 in (b)).
Colors in (b) mark the amount of above-cloud water vapor. Statistics in both panels summarize each approach's number of samples, bias, and standard deviation of radiance residuals as well as relative deviations.
3.1 State-of-the-Art Approach (Su et al., 2015) 100 An analytic sigmoidal function related TOA SW radiance with MODIS-based f andτ .
Where x = log fτ for a single cloud layer or x = log[(f 1 + f 2 )e f 1 logτ 1 +f 2 logτ 2 f 1 +f 2 ] for two layers and I 0 , a, b, c, and x 0 were free parameters. Optimization of sigmoidal parameters relied on mean radiances that were produced per x interval (every 0.02, shown as black dots in Figure 1a).
Selected angular bin in Fig. 1 had a sun-glint angle of about 14 • and shows how tabulated radiances (colors correspond to wind speed intervals) and sigmoidal curve both covered observed radiances.

115
There are several ways one might incorporate additional variables R e and ACW V into a radiance-predicting statistical model.
One could divide each angular bin's samples into classes of R e and ACW V and repeat sigmoidal fitting for each combination of classes (see Section 3.1). Some bins, however, contained too few samples or failed to cover the full spectrum of at least one of the two parameters. As a viable alternative, we explored radiative transfer approximations as a way to ingest scene properties (i.e. MODIS-based cloud properties and geophysical auxiliary parameters), and this allowed incorporating all samples in a 120 continuous manner.
Working with cloudy atmospheres over ocean surfaces, we assumed that radiance fluctuations were mainly driven by the bidirectional reflection of clouds and water surfaces and by directional absorption through water vapor located above (highly reflective) clouds. We initially set out with following simple equation of radiative transfer: 125 with solar influx S o cos θ 0 and the albedo α of an Earth-atmospheric scene covered by the CERES footprint (hereafter referred to as footprint albedo).
In following subsection, we present how footprint albedo was approximated. This then allowed us to use Eq. 3 in its loglinear form and weight the contribution of reflection and absorption via ordinary least square with free parameters A, B and C.
Like the state-of-the-art methodology (Section 3.1), we applied this approach per angular bin (resolved by 2 • in θ 0 , θ v , and ϕ) allowing us to treat S o cos θ 0 as constant. We also separated by cloud phase but choose a different threshold to discriminate phase. As elaborated in more detail below, we rely on pure liquid and ice phases to, then, treat the mixed phase. Therefore, we consider a footprint as liquid phase for φ 1/2 = 1, as ice for φ 1/2 = 2, and as mixed for φ 1 = 1 and φ 1 = 2. φ 1/2 were rounded 135 in case their values were neither 1 or 2.

Approximating CERES Footprint Albedo
To approximate the albedo within each CERES footprint by means of MODIS-based cloud properties and additional geophysical variables (α surf ace , w 10m , ACW V ), we separately handled clear and cloudy portions.
For clear portions within each footprint, we used the surface broadband albedo of underlying water bodies (referred to 140 α ocean = α surf ace , see Section 2). To capture sun-glint, i.e. the specular reflection at ocean's surface that alters as low-level winds perturb the water surface and tilt reflective facets, we used a Cox-Munk reflectance (Cox and Munk, 1954), as formulated in Wald and Monget (1983) with Fresnel reflection factor ρ(ω) for a perfectly smooth surface, and sun-observer-geometry per CERES footprint: where P (θ n , W 10m ) = 1 150 θ n = arccos cos θ v + cos θ 0 2 cos ω (8) cos 2ω = cos θ v cos θ 0 + sin θ v sin θ 0 cos ϕ To describe the albedo of cloudy portions, we explored the application of two-stream cloud albedo that uses cloud optical thickness and cloud micro-physical properties through asymmetry parameter g or backscattering fraction β. This allowed us to 155 ingest MODIS-basedτ and R e through g(R e ), as explained in more detail in the following subsection. The following solutions are thoroughly described in Meador and Weaver (1980), which presents a unifying theoretical framework to a variety of twostream cloud albedos based on coupled differential equations that describe upward and downward directed intensity fields. We considered two cloud albedos that proved useful for a range of cloud optical thicknesses (King and Harshvardhan, 1986): the Eddington approximation (Shettle and Weinman, 1970) and the Coakley-Chylek approximation (using solution I of Coakley 160 and Chylek, 1975).
The Eddington approximation considered an incident flux explicitly in coupled equations and thus described diffuse intensity fields. Assuming conservative scattering (i.e. a single scattering albedo of 1), a perfectly absorbing lower boundary (α bottom = 0), and no further influx at TOA, the analytical solution for cloud albedo was as follows, where µ 0 = cos θ 0 : The Coakley-Chylek approximation excluded the incident flux in differential equations and thus its intensities referred to total radiation fields (i.e. direct and diffuse). Assuming conservative scattering, a perfectly absorbing lower boundary (α bottom = 0) and only a solar influx at TOA, the analytical solution for cloud albedo was: For the same angular bin as in Figure 1, we present details of the proposed model that highlight essential steps aside from log-linear least-square fitting (Eq. 4). (a) shows the search for an optimal g(Re) (as described in Sec. 3.2): we plotted a two-dimensional slice (showing b and c of Eq. 14) through the three-dimensional space (spanned by a, b, and c). Colors show standard deviations of radiance residuals and point size relates to model bias. The star marks the combination of a, b, and c that produced smallest residual standard deviations and is considered optimal for this bin. (b) compares the g(Re) of the determined optimal solution against Mie-calculations. (c) shows final radiance residuals against cloud homogeneity (x-axis) and cloud optical depth (color). As described in Section 3.2.2., only homogeneous (ν > 10) clouds which were well-retrieved (MODIS-reported portion > 90%) -marked as triangles in (c) -were conisidered for opimization of g(Re) and least-square fitting. Statistcs and error metrics throughout the manuscript incorporate all samples.
where β was substituted with (1−g) 2 as done in textbook solutions (e.g. Bohren and Clothiaux, 2008). Using the Coakley-Chylek 170 approximation and a reflective lower boundary with albedo α bottom > 0 (in this study α bottom = α surf ace ), we produced following cloud albedo: Because it was unclear which solution could explain radiance fluctuations over narrow sun-observer-geometries most successfully, we tested a variety of solutions in Section 4.

175
Both ocean and cloud albedos (for up to two cloud layers) were used to calculate the footprint albedo, using clear fraction f 0 and cloud fractions of layer 1 and layer 2, f 1 and f 2 , respectively: where f 0 + f 1 + f 2 = 1.

180
Before comparing different two-stream approximations in Sec. 4, we performed two steps that ensured statistical optimization for each approximation. Finding an optimal g(R e ) was designed to best capture radiance fluctuations per angular bin. Higher Table 2. In search for optimal g(Re), we list the range (Minimum and Maximum) and step size for each parameter in Eq. 14.

Parameter Minimum Maximum
Step weights to a subset of data per angular bin -homogeneous clouds that were well retrieved -was used to facilitate consistency of radiances across bins. Both steps are explained in more detail below.
As shown in the previous subsection, we used two-stream cloud albedo to explain radiance fluctuations for narrow sun-185 observer-geometries. Applied to all angular bins of an upward hemisphere, it was unclear which g(R e ) to choose. Initial tests that used a g(R e ) from Mie theory (see Fig. 2b) for all geometries proved sub-optimal for some angular bins and left radiance residuals correlated to layer mean effective radius (not shown). We therefore decided to optimize g(R e ) for each angular bin and for each cloud phase (liquid and ice). Inspired by the shape of Mie-calculated g(R e ), we approximatedg(R e ) via a quadratic function: and searched a three-dimensional grid, spanned by a, b, and c, for combinations that minimized the standard deviation of radiance residuals. The search covered parameters a, b, and c as listed in Table 2. As shown in Fig. 2a, we usually found a single optimum value that could minimize standard deviation of radiance residuals and that deviated from Mie calculations ( Fig. 2b).

195
A second step aimed at using a subset of data that was consistent across angular bins. Looking into samples of individual angular bins, we observed stark variability in radiances that could be attributed to cloud horizontal heterogeneity (cloud homogeneity was approximated by ν =τ 2 σ(τ ) 2 ; radiance residuals are shown in Fig. 2c). We suspected that clouds' three-dimensional structure caused tilted cloud facets that led to more or less reflective cloud portions (e.g. as accounted for in Scheck et al., 2018). In order to avoid an uncontrollable impact of cloud heterogeneity onto final models, we decided to select homogeneous 200 samples only for statistical optimization. As a threshold of homogeneity, we used ν > 10 (e.g. Barker et al., 1996;Kato et al., 2005). As shown in Tab. 3 per solar geometry, median homogeneity varied considerably across bins as well as cloud phase, and this resulted in ranging portions of data being selected. For optimization, we further limited selection to CERES footprints with quality flags indicating a confident retrieval of 80% or more of all cloudy MODIS pixels within a CERES footprint. This subset of samples served to optimize the above search for g(R e ) and to find weights via least-square (Eq. 4). To compute error 205 metrics, we used all available samples. An example for the application of the log-linear model in shown in Fig. 1b.  Radiance-predicting statistical models that capture narrow sun-observer-geometries form the basis for empirical angular distribution models. And these statistical models fit observations from satellites, typically capturing how TOA SW radiances measured by a broadband radiometer change with scene type (defined by surface conditions as well as cloud and aerosol 210 properties within the radiometer's footprint area) retrieved using a multi-spectral imager (see Sec. 2). To investigate whether a new approach, the proposed semi-physical log-linear model in Sec. 3.2, is a superior way to fit observations compared to the state-of-the-art approach, the sigmoidal fit described in Sec. 3.1, we took CERES Ed4SSF observations (Sec. 2) of liquidphase clouds along the principal plane of an exemplary solar geometry covering major scattering features of clouds and the ocean surface. We applied the sigmoidal fit as well as a variety of log-linear models, each using a different analytic solution 215 of two-stream cloud albedo (Eqs. 10-12) that is used in this study as a framework to ingest MODIS-based cloud properties.
Looking at the standard deviation of radiance residuals per angular bin (in this study used as a measure of model uncertainty), the Coakley-Chylek approximation using a reflective lower boundary (Eq. 12) outperformed the sigmoidal fit for most bins and by up to 1.5 W m −2 sr −1 (shown in Fig. 3). Only the central portion of sun-glint-affected geometries remained best explained by sigmoidal fits (and accompanied look-up-table approach as laid out in Section 3.1). In contrast, the Coakley-Chylek 220 approximation using a perfectly absorbing lower boundary (Eq. 11) or the Eddingtion approximation (Eq. 10) performed only equally well or worse than sigmoidal fits.
To ensure that the Coakley-Chylek approximation using a reflective lower boundary performed well for other Sun-observergeometries, we processed all angular bins that contained more than 100 samples. As Table 3 shows, this covered between 13 and 96% of all angular bins. We found that for liquid clouds (top panels of Fig. 4) and θ 0 ∼ 20 • -70 • more than half of the 225 bins were better explained by the log-linear approach and errors were reduced by up to 13.2%. For solar geometries θ 0 > 40 • , bins in sun-glint-affected geometries (constituting a portion of all bins in a hemisphere between 10% for a θ 0 ∼ 20 • and 1% for a θ 0 ∼ 75 • ) caused higher uncertainties in log-linear models. For solar geometries θ 0 < 20 • , on the other hand, we found bins outside the sun-glint -i.e. mostly slant observation angles -were best treated with the sigmoidal approach. Few footprints (indicated by circle size) of the top row were treated as mixed in the log-linear model and will be evaluated further below. With 230 these limitations in mind, we use the Coakley-Chylek approximation using a reflective lower boundary standard as two-stream cloud albedo for the remainder of this study.
To determine whether the log-linear approach predicted plausible radiance fields, we tested it on a variety of scenarios. When applied to a range of cloud optical thicknesses, we found a similar radiance response compared to sigmoidal fit (Fig. 5b). Setting cloud fraction to zero (f 1 = f 2 = 0) and using a range of 10 m wind speeds, log-linear and sigmoidal models produced again 235 comparable radiance fields (Fig. 5c). This shows that the sensitivities of the state-of-the-art approach were captured by loglinear models. When varying cloud-top effective radius -a newly added sensitivity -we found radiances grow as droplet size increased (leaving cloud optical thickness constant; shown in Fig. 5e). With a focus on single-scattering features, we found the cloud glory (centered around the direct backscatter) to widen and the cloud glory (positioned about 20 • away from the backscatter) to shift towards the direct backscatter as effective radii became smaller. This observation is corroborated by Mie-calculations of scattering phase functions (e.g. Fig. 1 in Tornow et al., 2018). The newly introduced concept of bin-wise optimized asymmetry parameters (Sec. 3.2.2) made changing cloud bow and glory possible and g(R e ) exhibited a symmetry left and right of the direct backscatter between θ v of -50 • and 0 • (Fig. 5f). For a range of above-cloud water vapor (Fig. 5d) -another newly added sensitivity -we observed that smaller loads produced higher radiances and found a slight increase in sensitivity with larger θ v .

245
We also tested log-linear models on observations of ice-phase clouds. We found that model uncertainties outside the sun glint were of similar magnitude as sigmoidal fits (Fig. 4, bottom panels). Possible reasons will be discussed in Sec. 5. Like the liquid-phase, predicted radiances increased with smaller ice crystal radii. However, distinct scattering features were absent (not shown); possibly a result of ice clouds' rich variety of crystal shapes (e.g. Zhang et al., 1999;Baum et al., 2005) that was unaccounted for. The response to above-cloud water vapor was consistent and covered much of the lower levels (0.03-0.17 kg 250 m −2 , see Table 3).
Roughly 50% of all CERES footprints cover both a liquid and an ice cloud and have been treated separately as "mixedphase". The proposed log-linear approach allows us to handle mixed-phase cases fundamentally differently. Instead of a footprint-effective optical depth (as used in Equ. 2), we can produce a footprint-effective albedo (Equ. 13) and account not only for cloud macro-physical (f 1/2 ,τ 1/2 ) but also for microphysical (R e1/2 ) changes. Optimized asymmetry parameters from   pure liquid and pure ice cases (Fig. 6b) were reused to describe the cloud albedo of respective cloud phase within each mixedphase CERES footprint. Hence only A, B, and C from Eq. 4 needed to be estimated. Fig. 6a illustrates the reduction in model uncertainty for many bins and of up to 2.5 W m −2 sr −1 when using the log-linear approach. Once again, the center of sun-glint remained best captured by the sigmoidal approach. Using a cloud-phase-specific albedo allowed us to account for radiance changes with varying amount of liquid versus ice fraction within a footprint. Fig. 6c shows radiance predictions for different 260 liquid-ice-proportions (which could not be captured by the state-of-the-art approach) and that both approaches agree for 50% liquid and 50% ice cloud footprints. Fig. 6d shows the sigmoidal fit's sensitivity to ranging cloud optical depth was captured by the log-linear approach. Looking at all available sun-observer-geometries (Fig. 4, middle panels) for solar geometries between θ 0 ∼ 20 • -70 • , we found model uncertainty of most bins reduced by as far as 35.8%.
In summary, we showed that the proposed log-linear model had the ability to outperform the existing sigmoidal approach in 265 capturing CERES radiance fluctuations per angular bin. It produced lower uncertainties, added new radiance sensitivities, and allowed to treat mixed-phase footprints in a fundamentally different manner.

Conclusions
Statistical models that capture measurements of TOA SW radiances as a function of corresponding scene type for narrow sun-observer-geometries are the basis for Angular Distribution Models. In this study, we introduced a new alternative that 270 incorporated additional parameters -namely cloud-top effective radius and above-cloud water vapor -via a semi-physical loglinear approach. We found this new approach to better explain radiance fluctuations for the majority of observed geometries and to produce plausible radiance fields.
Incorporating additional parameters that help explain radiance fluctuations may have minimized sampling bias. Ranges in effective radius or above-cloud water vapor varied across bins and ignoring this variation can cause a radiance bias in individual 275 angular bins. Even accounting for parameters that may not affect TOA anisotropy, such as cloud horizontal heterogeneity, has the potential to minimize sampling biases. We found varying portions of heterogeneous samples across bins and suspect that their variation in radiance (cf. Fig. 2c) failed to cancel out. Thus, giving higher (or all) weight to homogeneous samples during regression, as done in this study, should eliminate any sampling bias.
The inclusion of cloud-top effective radius and above-cloud water vapor was successful as evidenced by reduced radiance 280 residuals and credible radiance fields. We failed to reduce radiance residuals for ice-phase clouds and made the following observations looking at ice cloud samples. First, among collected observations, we found footprints to mostly contain homogeneous ice clouds. Second, ice clouds had only small loads of water vapor aloft. Lastly, there was an absence of distinct single-scattering features. We suspect that these are characteristics that drive potential reduction of radiance residuals and that liquid clouds samples, having near asymmetric properties (few homogeneous samples, large loads of water vapor aloft, distinct 285 scattering features), benefitted especially from this new approach.
We successfully used a theoretic framework -inspired by radiative transfer approximations designed for hemispheric averages -and applied it to narrow sun-observer-geometries. A derived byproduct, the asymmetry parameter g(R e |θ ∆ 0 , θ ∆ v , ϕ ∆ ), captured observer-specific multi-scattering. Could this byproduct contain information that allows inference on multi-scattering properties? Monte-Carlo radiative transfer simulations may help in answering this. Future work should simulate radiances, de-290 rive simulation-based g(R e |θ ∆ 0 , θ ∆ v , ϕ ∆ ) and extract additional properties, such as photon path length or number of scattering events.
Statistical models allow finding scene properties that produce similar radiative responses (often referred to as similarity conditions). Like the state-of-art-approach, where different combinations of cloud fraction and cloud optical thickness produced similar radiances, the new semi-physical approach added cloud particle size and above-cloud absorber mass to parameter 295 combinations. A similarity condition explaining albedo through adjusted optical thickness, (1 − g)τ , was found earlier using simulations (e.g. van de Hulst, 1996). To our knowledge, this is the first time adjusted optical thickness (here employed in the framework of two-stream albedo) has been used to capture similarities of observed radiances.
The proposed semi-physical approach can easily be applied to land surfaces. Imager-based bidirectional reflectance distribution function (BRDF) products, such as MCD43GF (MODIS BRDF/albedo/nadir BRDF-adjusted reflectance Climate 300 Modeling Grid gap-filled; Moody et al., 2008), could provide land surface albedo and surface bi-directional reflectance in order to determine each observation's footprint albedo. Future efforts should test if this application over land can compete with CERES' separate treatment by latitude-longitude boxes. Recent efforts that demonstrated circumvention of this regional separation for clear-sky ADMs by using MCD43GF instead indicated a positive outcome (Tornow et al., 2019).
Lastly, we hope this new log-linear approach will form the basis of future angular distribution models. In particular, we 305 expect that cloudy scenes of microphysical extremes (i.e. clouds consisting of very small or very large droplets) observed from the backscattering direction will benefit from radiance-to-flux conversion using new models. More accurate estimates of instantaneous fluxes should benefit EarthCARE' studies of cloud-radiative processes regarding both water and energy fluxes.
We are currently examining this impact on instantaneous fluxes as well as the propagation of updated flux estimates into daily and monthly flux products.

310
Author contributions. FT had the idea, designed the experiment, and conducted the research. CD, HB, and RP had major influence on the development of the methodology through discussion. CD and HB further helped revising this manuscript. JF provided essential resources for data processing and evaluation.

Competing interests. No competing interests.
Acknowledgements. We thank Jason N. S. Cole, Almudena Velazquez Blazquez, Tobias Wehr, and all other members of the CLARA team as 315 well as colleagues at Institute for Space Sciences at Freie Universität Berlin for helpful discussion. We further thank the Atmospheric Sciences Data Center at the National Aeronautics and Space Administration Langley Research Center for providing the Clouds and Earth's Radiant