Articles | Volume 13, issue 7
Atmos. Meas. Tech., 13, 3909–3922, 2020
Atmos. Meas. Tech., 13, 3909–3922, 2020

Research article 21 Jul 2020

Research article | 21 Jul 2020

Using two-stream theory to capture fluctuations of satellite-perceived TOA SW radiances reflected from clouds over ocean

Using two-stream theory to capture fluctuations of satellite-perceived TOA SW radiances reflected from clouds over ocean
Florian Tornow1,2,3, Carlos Domenech4, Howard W. Barker5, René Preusker1, and Jürgen Fischer1 Florian Tornow et al.
  • 1Institute for Space Sciences, Freie Universität Berlin, Berlin, Germany
  • 2Center for Climate Systems Research, Earth Institute, Columbia University, New York, NY, USA
  • 3NASA Goddard Institute for Space Studies, New York, NY, USA
  • 4GMV, Madrid, Spain
  • 5Environment and Climate Change Canada, Toronto, Ontario, Canada

Correspondence: Florian Tornow (


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 semiphysical model predicting radiances per narrow sun-observer geometry. This model was statistically trained using CERES-measured radiances paired with MODIS-retrieved cloud parameters as well as reanalysis-based geophysical parameters. By using radiative transfer approximations as a framework to ingest the 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 two-stream cloud albedo – serving to statistically incorporate cloud optical thickness and cloud-top effective radius – and Cox–Munk ocean reflectance were used to describe an albedo over each CERES footprint. Effective-radius-dependent asymmetry parameters were obtained empirically and separately for each viewing-illumination geometry. 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. We identified the two-stream functional form that minimized radiance residuals calculated against CERES observations and outperformed the state-of-the-art approach for most observer geometries outside the sun-glint and solar zenith angles between 20 and 70, reducing the median SD 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. Geometries affected by sun glint (constituting between 10 % and 1 % of the discretized upward hemisphere for solar zenith angles of 20 and 70, respectively), however, often showed weaker performance when handled with the new approach and had increased residuals by as much as 60 % compared to the state-of-the-art approach. Overall, uncertainties were reduced for liquid-phase and mixed-phase footprints by 5.76 % and 10.81 %, respectively, while uncertainties for ice-phase footprints increased by 0.34 %. 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 particular dealing with clouds characterized by small or large droplet/crystal sizes.

1 Introduction

Radiative fluxes at the top of the atmosphere (TOA) inferred from satellite observations serve many purposes. Instantaneous flux estimates paired with properties of underlying clouds, aerosols, atmospheric gases, and Earth's surface may inform us about the radiative effect of each component of Earth's radiation budget (e.g., Loeb and Manalo-Smith2005; Li et al.2011; Thorsen et al.2018). TOA fluxes may also help to constrain uncertainties concerning cloud–aerosol–radiation interactions, which will be tested in the EarthCARE satellite mission (Illingworth et al.2015). In EarthCARE's radiative closure assessment, observation-based fluxes will be used to help continuously assess both active–passive retrievals of cloud and aerosol properties and results from radiative transfer simulations performed on them (Barker et al.2011; Barker and Wehr2012). Integrals of estimated fluxes over large areas and long time spans (Loeb et al.2018) help us understand the Earth–atmosphere system's current radiation budget (e.g., Stephens et al.2012), thus helping to verify global climate models (e.g., Bender et al.2006; Boucher et al.2013; Calisto et al.2014; Nam et al.2012).

Inferring fluxes from satellite-based radiometry involves a number of steps. The key challenge for solar fluxes, the general focus of this paper, is that constituents such as clouds reflect solar radiation unevenly across the upward hemisphere, and we need to assume how measurements from a subset of directions relate to radiances in directions not viewed. The intention is to adequately represent hemispheric distributions of radiances such that when integrated yield accurate flux estimates. The solution to this challenge has been empirical angular distribution models (ADMs) that learn, via statistical approaches, hemispherically resolved radiance fields associated with atmospheric scenes using years of satellite observations. For clouds over ocean, the specific concern of this paper, early efforts (Suttles et al.1988; Smith et al.1986) worked with ERBE (Earth Radiation Budget Experiment) radiometry as well as GOES (Geostationary Operational Environmental Satellite) measurements and defined four scene types ranging in cloud coverage (including “clear ocean”, which used a cloud cover up to 5 %; two cloudy scene types over ocean; and “overcast”, which blended all surface types). Observations were sorted to produce mean radiances per observed angular ranges for each illumination geometry. Using CERES (Clouds and the Earth's Radiant Energy System) and VIRS (Visible and Infrared Scanner) on the TRMM (Tropical Rainfall Measuring Mission) satellite, Loeb et al. (2003) refined this method and sorted observations into combinations of 12 cloud coverage classes and 14 cloud optical thickness groups and treated ice- and liquid-phase clouds separately. Instead of a discrete scene-type definition, Loeb et al. (2005) defined a continuous description of scene type for the Terra mission, using a sigmoidal function to fit cloud optical thickness and cloud fraction based on MODIS (Moderate Resolution Imaging Spectroradiometer) measurements with CERES-measured TOA shortwave (SW) radiances. They treated footprints containing both ice- and liquid-phase clouds (throughout the paper referred to as “mixed phase”) separately from pure and ice and liquid cases. Much of their state-of-the-art methodology was adapted for the Aqua mission (also hosting CERES and MODIS instruments) using improved cloud algorithms and longer data records (Su et al.2015).

A recent case study (Tornow et al.2018) focused on marine stratocumulus-like clouds of optical thickness τ̃10 and identified additional parameters that influence ADMs: above-cloud water vapor (ACWV) and layer mean cloud-top effective radius Re. They showed that ignoring these parameters could cause deviations in instantaneous flux estimates of about 10 Wm−2. This suggests the non-negligible role of solar absorption and single scattering for determination of cloud reflectance patterns. Features of single scattering, such as the cloud bow and glory for liquid clouds or the specular reflection peak for ice clouds, were generally visible in earlier ADMs (e.g., Loeb et al.2005). These features – solely shaped by the particle phase function that largely depends on particle shape and size – can occur for a wide range of cloud optical thicknesses. Using simulated radiance fields, Gao et al. (2013) demonstrated that scattering regimes, ranging from foremost single scattering to Lambertian-like multiscattering mediums, are functions of the cloud optical thickness. For an intermediate regime, which showed single-scattering features, Gao et al. speculated that the uppermost τ≈1 of cloud is responsible for single-scattering contributions.

This study presents a novel semiphysical model that predicts TOA SW radiances for cloudy scenes over ocean for narrow ranges of sun-observer angles. Estimates are sensitive to Re and ACWV and are compared to results from the state-of-the-art methodology. This new approach used the two-stream approximation to statistically ingest MODIS cloud properties and other geophysical auxiliary parameters. We began by finding the framework of approximations that best explained CERES-observed radiance fluctuations and then demonstrated that semiphysical log-linear models produced tenable radiance fields.

Section 2 presents data from Aqua and Terra satellites used in the current study. Section 3 explains both the state-of-the-art methodology for radiance estimation and the new approach. Section 4 identifies optimal solutions and assesses their properties. Section 5 discusses results and conclusions.

2 Data

Measured TOA SW radiances paired with scene properties – including imager-based cloud properties and further geophysical auxiliary parameters – were obtained from the CERES Ed4SSF (Edition 4.0 Single Scanner Footprint) dataset of Aqua and Terra satellite missions, primarily from days during years 2000–2005 when CERES instruments were measuring in rotating azimuth plane scan mode to provide angular coverage for ADM construction.

We extracted parameters concerning CERES broadband radiometry. Apart from upwelling unfiltered TOA SW radiances I*, covering the spectral range of 0.4–4.5 µm, and their angular geometry (i.e., solar zenith angle θ0, viewing zenith angle θv, and relative azimuth angle φ), we collected downwelling TOA SW fluxes F that incorporate each measurements' prevalent Sun–Earth distance, which allowed normalization of gathered radiances via I=I*S0cosθ0F, with solar constant S0=1361.0Wm−2.

Collocated to each CERES footprint, the SSF dataset summarizes cloud property retrievals (Sun-Mack et al.2018) on the MODIS pixel level taking into account the CERES point spread function (PSF) (Wielicki et al.1996) and reports properties for up to two cloud layers per footprint (given that both layers' cloud-top pressure differed by 50 hPa or more; Loeb et al.2003). We extracted layer cloud fraction f and several statistics on the retrieved field of cloud optical depth τ (layer average of its logarithm τ̃=elogτ, layer average τ¯, and layer SD σ(τ)), as well as layer mean values of cloud condensate phase ϕ (i.e., liquid, ice, or a mixture of both; involving MODIS band 3.7 µm), effective radii of water or ice particles Re (using band 3.7 µm), and cloud-top pressure pctop. A quality flag summarizing the retrieval confidence (using the parameter “Note for cloud layer” from the SSF dataset) was also collected.

Additional geophysical auxiliary parameters provided in the SSF dataset were extracted. We obtained a surface broadband albedo αsurface, surface IGBP (International Geosphere-Biosphere Programme) types, and 10 m surface wind speed w10 m. The wind speed parameter stemmed from GEOS data assimilation version 5.4.1.

Lastly, to incorporate above-cloud water vapor (ACWV) into our analysis, we used layer mean cloud-top pressure (of the layer with larger cloud fraction) and extracted from ERA-20C (ECMWF 20th century) reanalysis (Poli et al.2016) four-dimensional fields whose vertical profiles of relative humidity RH(p) and temperature T(p) that were nearest in time and geolocation to the footprint center. For each CERES footprint we collocated the following vertical integral of mixing ratio mr(p), with saturation vapor pressure es=6.112e17.67TT+243.5 (using T in C(Bolton1980), gravitational acceleration g, and molecular weights of water and dry air molh2o and molair respectively:


For footprints consisting of multiple cloud layers, relying on a single cloud-top pressure may introduce uncertainty, especially for mixed-phase footprints (see Sect. 3) where pressure difference between ice- and liquid-phase layers is exceptionally large.

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.

Table 1Number of CERES footprints obtained after screening for marine clouds. Number are shown in millions (in total 1 711 937 663 footprints).

Download Print Version | Download XLSX

3 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 geometry. Following Su et al. (2015), solar zenith, viewing zenith, and 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. Section 3.2 introduces a novel semiphysical approach that includes additional parameters.

3.1 State-of-the-art approach (Su et al.2015)

An analytic sigmoidal function related TOA SW radiance with MODIS-based f and τ̃.

(2) I θ 0 Δ , θ v Δ , φ Δ = I 0 + a 1 + e - ( x - x 0 ) b c ,

where x=logfτ̃ for a single cloud layer or x=log[(f1+f2)ef1logτ̃1+f2logτ̃2f1+f2] for two layers, and I0, a, b, c, and x0 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 Fig. 1a).

Figure 1For an example angular bin (θ0[20,22], θv[6,8], φ[12,14]), we show how a state-of-the-art sigmoidal fit (a) and proposed log-linear model (b) capture fluctuations of CERES-measured TOA SW radiances. As this angular bin is within the sun-glint region, panel (a) shows the lookup-table approach for x<6 (as defined in Sect. 3.1; note that f1 and f2 are taken between 0 and 100 in panel a and between 0 and 1 in panel 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 SD of radiance residuals as well as relative deviations.


Models were generated separately per cloud phase. A footprint's cloud phase was determined via an effective phase, defined as ϕeff=f1ϕ1+f2ϕ2f1+f2 for two layers, and the following thresholds: liquid for 1<ϕeff<1.01, mixed for 1.01ϕeff1.75, and ice for 1.75<ϕeff2.

To handle radiance fluctuations caused by sun glint, a glint region was defined (sun-glint angles <20). Observations with x>6 in affected geometries remained captured by a sigmoid fit. For x≤6, on the other hand, a lookup-table approach stored mean radiances per wind speed interval (0–2, 2–4, 4–6, 6–8, 8–10, and >10m s−1) and per x interval (<3.5, 3.5–4.5, 4.5–5.5, 5.5–6).

The 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 the sigmoidal curve both covered observed radiances.

3.2 Novel semiphysical approach

There are several ways one might incorporate additional variables Re and ACWV into a radiance-predicting statistical model. One could divide each angular bin's samples into classes of Re and ACWV and repeat sigmoidal fitting for each combination of classes (see Sect. 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 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 the following simple equation of radiative transfer:

(3) I θ 0 Δ , θ v Δ , φ Δ S o cos θ 0 α e - 2 ACWV ,

with solar influx Socos θ0 and the albedo α of an Earth–atmospheric scene covered by the CERES footprint (hereafter referred to as footprint albedo).

In the following subsection we present how footprint albedo was approximated. This then allowed us to use Eq. (3) in its log-linear form and weight the contribution of reflection and absorption via ordinary least-square fitting with free parameters A, B, and C.

(4) log I θ 0 Δ , θ v Δ , φ Δ A + B log α + C ACWV

Like the state-of-the-art methodology (Sect. 3.1), we applied this approach per angular bin (resolved by 2 in θ0, θv, and φ), allowing us to treat Socos θ0 as constant. We also separated by cloud phase but chose 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 ϕ2=2. ϕ1∕2's were rounded if their values were neither 1 nor 2.

3.2.1 Approximating CERES footprint albedo

To approximate the albedo within each CERES footprint by means of MODIS-based cloud properties and additional geophysical variables (αsurface, w10 m, ACWV), 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 as αocean=αsurface; see Sect. 2). To capture sun glint, i.e., the specular reflection at the ocean's surface that is altered as low-level winds perturb the water surface and tilt reflective facets, we used a Cox–Munk reflectance (Cox and Munk1954), as formulated in Wald and Monget (1983) with Fresnel reflection factor ρ(ω) for a perfectly smooth surface, and sun-observer geometry per CERES footprint:

(5) r SunGlint = π ρ ( ω ) P ( θ n , W 10 m ) 4 cos θ 0 cos θ v cos 4 θ n ,



To describe the albedo of cloudy portions, we explored the application of two-stream equations as a function to ingest MODIS-retrieved cloud optical thickness τ̃ and cloud-top effective radius Re through asymmetry parameter g(Re) or backscattering fraction β(Re), 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 for a variety of two-stream 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 Harshvardhan1986): the Eddington approximation (Shettle and Weinman1970) and the Coakley–Chylek approximation (using solution I of Coakley and Chylek1975).

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:

(10) α TwoStream = 3 4 ( 1 - g ) τ ̃ - 1 4 ( 1 - 3 μ 0 ) 1 - e - τ ̃ μ 0 1 + 3 4 ( 1 - g ) τ ̃ .

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

(11) α TwoStream = ( 1 - g ) τ ̃ 2 1 + ( 1 - g ) τ ̃ 2 ,

where β was substituted with (1-g)2 as done in textbook solutions (e.g., Bohren and Clothiaux2008). Using the Coakley–Chylek approximation and a reflective lower boundary with albedo αbottom>0 (in this study αbottom=αsurface), we produced following cloud albedo:

(12) α TwoStream = α ocean + ( 1 - α ocean ) ( 1 - g ) τ ̃ 2 1 + ( 1 - α ocean ) ( 1 - g ) τ ̃ 2 .

Because it was unclear which solution could explain radiance fluctuations over narrow sun-observer geometries most successfully, we tested a variety of solutions in Sect. 4.

Both ocean and cloud albedos (for up to two cloud layers) were used to calculate the footprint albedo, using clear fraction f0 and cloud fractions of layer 1 and layer 2, f1 and f2, respectively:

(13) α = f 0 α ocean + r SunGlint + f 1 α 1 TwoStream + f 2 α 2 TwoStream ,

where f0+f1+f2=1.

3.2.2 Statistical optimization

Before comparing different two-stream approximations in Sect. 4, we performed two steps that ensured statistical optimization for each approximation. Finding an optimal g(Re) was shown to best capture radiance fluctuations per angular bin. Higher weights to a subset of data per angular bin – homogeneous clouds that were well retrieved – were used to facilitate consistency of radiances across bins. Both steps are explained in more detail below.

Figure 2For the same angular bin as in Fig. 1, we present details of the proposed model that highlight essential steps aside from log-linear least-square fitting (Eq. 4). Panel (a) shows the search for an optimal g(Re) (as described in Sect. 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 the SD of radiance residuals and point size relates to model bias. The star marks the combination of a, b, and c that produced the smallest residual SD and is considered optimal for this bin. Panel (b) compares the g(Re) of the determined optimal solution against Mie calculations. Panel (c) shows final radiance residuals against cloud homogeneity (x axis) and cloud optical depth (color). As described in Sect. 3.2.2., only homogeneous (ν>10) clouds which were well retrieved (MODIS-reported portion ≥80 %) – marked as triangles in panel (c) – were conisidered for optimization of g(Re) and least-square fitting. Statistics and error metrics throughout the paper incorporate all samples.


As shown in the previous subsection, we used two-stream cloud albedo to explain radiance fluctuations for narrow sun-observer geometries. Applied to all angular bins of an upward hemisphere, it was unclear which g(Re) to choose. Initial tests that used a g(Re) from Mie theory (see Fig. 2b) for all geometries proved suboptimal for some angular bins and left radiance residuals correlated to layer mean effective radius (not shown). We therefore decided to optimize g(Re) for each angular bin and for each cloud phase (liquid and ice). Inspired by the shape of Mie-calculated g(Re), we approximated g(Re) via a quadratic function,

(14) g ( R e ) = a + b R e + c R e 2 ,

and searched a three-dimensional grid, spanned by a, b, and c, for combinations that minimized the SD 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 the SD of radiance residuals and that deviated from Mie calculations (Fig. 2b).

Table 2In search for optimal g(Re), we list the range (minimum and maximum) and step size for each parameter in Eq. (14).

Download Print Version | Download XLSX

A second step aimed at using a subset of data that were 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 on final models, we decided to select homogeneous 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 Table 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(Re) and to find weights via least-square fitting (Eq. 4). To compute error metrics, we used all available samples. An example for the application of the log-linear model in shown in Fig. 1b.

Table 3Per solar geometry θ0 and per cloud phase (L – liquid, I – ice, M – mixed) as defined in Sect. 3.2, we show what portions of the upward hemisphere were covered with observations and how large the range of cloud homogeneity, above-cloud water vapor, and cloud-top effective radius was. The range lists minima and maxima of median values computed per angular bin within a hemisphere.

Download Print Version | Download XLSX

4 Results

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 properties within the radiometer's footprint area) retrieved using a multispectral imager (see Sect. 2). To investigate whether a new approach, the proposed semiphysical log-linear model in Sect. 3.2, is a superior way to fit observations compared to the state-of-the-art approach, the sigmoidal fit described in Sect. 3.1, we took CERES Ed4SSF observations (Sect. 2) of liquid-phase clouds along the principal plane of an example 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 of two-stream cloud albedo (Eqs. 1012) that is used in this study as a framework to ingest MODIS-based cloud properties. Looking at the SD 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 by up to 1.5 Wm-2sr-1 (shown in Fig. 3). Only the central portion of sun-glint-affected geometries remained best explained by sigmoidal fits (and accompanied lookup-table approach as laid out in Sect. 3.1). In contrast, the Coakley–Chylek approximation using a perfectly absorbing lower boundary (Eq. 11) or the Eddington approximation (Eq. 10) performed only equally well or worse than sigmoidal fits.

Figure 3Applied to angular bins of the principal plane for θ0[20,22], we test a variety of two-stream solutions for cloud albedo (Eqs. 1012) as input to log-linear models as presented in Eq. (4). This plot shows the SD of resulting residuals and compares against the state-of-the-art sigmoidal fit (black line). As labeled, gray dashed lines mark the position of sun-glint and direct backscatter.


To ensure that the Coakley–Chylek approximation using a reflective lower boundary performed well for other sun-observer geometries, 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 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 θ020 and 1 % for a θ075) caused higher uncertainties in log-linear models, increasing with solar zenith angle and higher by up to 60 % compared to the sigmoidal approach. 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. A 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 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.

Figure 4Using all observed angular bins within θ0[6,82], we show how radiance residuals from proposed log-linear models compare against state-of-the-art sigmoidal fits. Results are presented by the CERES-defined cloud phase (vertically), by the newly defined phase (colors), and by whether the angular geometry is affected by sun glint (left) or free of sun glint (right). We show relative change in model uncertainty: δ=[σ(ΔILogLinear)-σ(ΔISigmoidal)]/σ(ΔISigmoidal)100%. Consequently, negative values relate to a better performance of the log-linear model, while positive values mark a better performance by the state-of-the-art methodology. Solid lines and dots mark the 50th percentile and shades show the interquartile range between the 25th and 75th percentiles. Point size relates to the average number of observations per angular bin. The dashed black line marks zero change.


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 the sigmoidal fit (Fig. 5b). Setting cloud fraction to zero (f1=f2=0) and using a range of 10 m wind speeds, log-linear and sigmoidal models produced again comparable radiance fields (Fig. 5c). This shows that the sensitivities of the state-of-the-art approach were captured by log-linear models. When varying cloud-top effective radius – a newly added sensitivity – we found radiances grew 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 (Sect. 3.2.2) made changing cloud bow and glory possible, and g(Re) 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.

Figure 5For angular bins along the principal plane for θ0[20,22] containing liquid-phase footprints, we present error metrics and sensitivities of proposed log-linear vs. state-of-the-art sigmoidal fits. Panel (a) shows the SD of residuals; colors mark the type of fit. Panel (b) displays the optimal g(Re) for three Re (by color). Panels (c), (d), and (e) demonstrate predicted radiances by both fits for varying cloud optical thickness (c), cloud-top effective radius (d), and above-cloud water vapor (e). Predictions from log-linear fits are colored, while predictions from sigmoidal fits are shown in black. Panel (f) presents the response of both fits to a variety of surface wind speeds. Properties held constant in panels (c), (d), (e), and (f) are listed in the each panel's top-right corner.


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 to sigmoidal fits (Fig. 4, bottom panels). Possible reasons will be discussed in Sect. 5. Similar to liquid-phase clouds, angular geometries affected by sun glint showed worse performance when using log-linear models, increasing residuals by up to 30 %. 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), which was unaccounted for. The response to above-cloud water vapor was consistent and covered much of the lower levels (0.03–0.17 kg 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 mixed phase. The proposed log-linear approach allows us to handle mixed-phase cases fundamentally differently. Instead of a footprint-effective optical depth (as used in Eq. 2), we can produce a footprint-effective albedo (Eq. 13) and account not only for cloud macrophysical (f1∕2, τ̃1/2) but also for microphysical (Re1/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 mixed-phase CERES footprint. Hence only A, B, and C from Eq. (4) needed to be estimated. Figure 6a illustrates the reduction in model uncertainty for many bins and of up to 2.5 Wm-2sr-1 when using the log-linear approach. Once again, the center of sun glint remained best captured by the sigmoidal approach, especially for solar zenith angles beyond 50 where semiphysical models produce up to 55 % higher residuals. Using a cloud-phase-specific albedo allowed us to account for radiance changes with varying amount of liquid vs. ice fraction within a footprint. Figure 6c shows radiance predictions for different 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 for the backscattering direction and 75 % liquid and 25 % ice cloud fractions for much of the forward scattering direction. indicating that sampled footprints shifted in liquid-to-ice proportions along the principal plane. Figure 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 and 70, we found model uncertainty of most bins reduced by as much as 35.8 %.

Figure 6For angular bins along the principal plane for θ0[20,22], we show details for mixed-phase footprints. Panel (a) presents the SD of residuals (colors mark the type of fit). Panel (b) shows optimal g(Re) from pure ice- and liquid-phase footprints employed for mixed-phase cases (colors mark liquid and ice particle effective radius). Panels (c) and (d) show predicted radiances for liquid cloud fraction (c) and cloud optical thickness (d). In panels (c) and (d), we show log-linear fits in color and sigmoidal fits in black. Quantities left constant are shown in the bottom-left corner.


Across all solar and viewing geometries, we calculated the median change in uncertainty when using the log-linear approach over the state-of-the-art approach to be 5.76 % for liquid-phase clouds, +0.34% for ice-phase clouds, and 10.81 % when both phases were present.

In summary, we showed that the proposed log-linear model had the ability to outperform the existing sigmoidal approach in capturing CERES radiance fluctuations per angular bin. For most geometries it produced lower uncertainties, added new radiance sensitivities, and allowed us to treat mixed-phase footprints in a fundamentally different manner. Drawbacks were typically found for geometries affected by sun glint.

5 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 incorporated additional parameters – namely cloud-top effective radius and above-cloud water vapor – via a semiphysical log-linear approach. We found this new approach to better explain radiance fluctuations for the majority of observed geometries and to produce plausible radiance fields. Weaker performance than the state-of-the-art approach was generally observed for solar zenith angles lower than 20 and for sun-glint-affected geometries that constitute between 1 % and 10 % of the hemispheric radiance field.

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 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 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 the potential reduction of radiance residuals and that liquid clouds samples, having near-asymmetric properties (few homogeneous samples, large loads of water vapor aloft, distinct 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(Re|θ0Δ,θvΔ,φΔ), captured observer-specific multiscattering. Could this byproduct contain information that allows inference on multiscattering properties? Monte Carlo radiative transfer simulations may help in answering this. Future work should simulate radiances, derive simulation-based g(Re|θ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 semiphysical approach added cloud particle size and above-cloud absorber mass to parameter combinations. A similarity condition explaining albedo through adjusted optical thickness, (1-g)τ̃, was found earlier using simulations (e.g., van de Hulst1996). 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 semiphysical approach can easily be applied to land surfaces. Imager-based bidirectional reflectance distribution function (BRDF) products, such as MCD43GF (MODIS BRDF and Albedo Global Gap-Filled, Snow-Free Version 6;, Schaaf2019), could provide land surface albedo and surface bidirectional 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 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.

Data availability

CERES Ed4SSF data were downloaded through (NASA Langley Research Center2020) and ERA20C was accessed via (ECMWF2020).

Author contributions

FT had the idea, designed the experiment, and conducted the research. CD, HWB, and RP had a major influence on the development of the methodology through discussion. CD and HWB further helped revise this paper. JF provided essential resources for data processing and evaluation.

Competing interests

The authors declare that they have no conflict of interest.


We thank Jason N. S. Cole, Almudena Velazquez Blazquez, Tobias Wehr, and all other members of the CLARA team as well as colleagues at the Institute for Space Sciences at Freie Universität Berlin for the 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 the Earth's Radiant Energy System Single Scanner Footprint TOA/Surface Fluxes and Clouds data product. We specifically thank Wenying Su for providing us with details on the CERES methodology. We further thank two anonymous referees very much for their feedback that helped to improve this paper substantially.

Financial support

This research has been supported by the European Space Agency (grant no. 4000112019/14/NL/CT).

We acknowledge support from the Open Access Publication
Initiative of Freie Universität Berlin.

Review statement

This paper was edited by Alexander Kokhanovsky and reviewed by two anonymous referees.


Barker, H. W. and Wehr, T.: Computation of Solar Radiative Fluxes by 1D and 3D Methods Using Cloudy Atmospheres Inferred from A-train Satellite Data, Surv. Geophys., 33, 657–676,, 2012. a

Barker, H. W., Wielicki, B. A., and Parker, L.: A Parameterization for Computing Grid-Averaged Solar Fluxes for Inhomogeneous Marine Boundary Layer Clouds. Part II: Validation Using Satellite Data, J. Atmos. Sci., 53, 2304–2316,<2304:APFCGA>2.0.CO;2, 1996. a

Barker, H. W., Jerg, M. P., Wehr, T., Kato, S., Donovan, D. P., and Hogan,  R. J.: A 3D cloud-construction algorithm for the EarthCARE satellite mission, Q. J. Roy. Meteor. Soc., 137, 1042–1058,, 2011. a

Baum, B. A., Heymsfield, A. J., Yang, P., and Bedka, S. T.: Bulk Scattering Properties for the Remote Sensing of Ice Clouds. Part I: Microphysical Data and Models, J. Appl. Meteorol., 44, 1885–1895,, 2005. a

Bender, F. A.-M., Rohde, H., Charloson, R. J., Ekman, A. M. L., and Loeb, N.: 22 views of the global albedo – comparison between 20 GCMs and two satellites, Tellus A, 58, 320–330,, 2006. a

Bohren, C. F. and Clothiaux, E. E.: Scattering: The Life of Photons, Chap. 3, John Wiley and Sons, Ltd, 125–184,, 2008. a

Bolton, D.: The Computation of Equivalent Potential Temperature, Mon. Weather Rev., 108, 1046–1053,<1046:TCOEPT>2.0.CO;2, 1980. a

Boucher, O., Randall, D., Artaxo, P., Bretherton, C., Feingold, G., Forster,  P., Kerminen, V.-M., Kondo, Y., Liao, H., Lohmann, U., Rasch, P., Satheesh,  S., Sherwood, S., Stevens, B., and Zhang, X.: Clouds and Aerosols, Sect. 7, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA,, 571–658, 2013. a

Calisto, M., Folini, D., Wild, M., and Bengtsson, L.: Cloud radiative forcing intercomparison between fully coupled CMIP5 models and CERES satellite data, Ann. Geophys., 32, 793–807,, 2014. a

Coakley, J. A. and Chylek, P.: The Two-Stream Approximation in Radiative Transfer: Including the Angle of the Incident Radiation, J. Atmos. Sci., 32, 409–418,<0409:TTSAIR>2.0.CO;2, 1975. a

Cox, C. and Munk, W.: Measurement of the Roughness of the Sea Surface from Photographs of the Sun's Glitter, J. Opt. Soc. Am., 44, 838–850,, 1954. a

European Centre for Medium-Range Weather Forecasts (ECMWF): ERA-20C Daily Dataset, available at:, last access: 17 July 2020. a

Gao, M., Huang, X., Yang, P., and Kattawar, G. W.: Angular distribution of diffuse reflectance from incoherent multiple scattering in tur bid media, Appl. Optics, 52, 5869–5879,, 2013. a

Illingworth, A. J., Barker, H. W., Beljaars, A., Chepfer, H., Delanoe, J., Domenech, C., Donovan, D. P., Fukuda, S., Hirakata, M., Hogan, R. J., Huenerbein, A., Kollias, P., Kubota, T., Nakajima, T., Nakajima, T. Y., Nishizawa, T., Ohno, Y., Okamoto, H., Oki, R., Sato, K., Satoh, M., Wandinger, U., Wehr, T., and van Zadelhoff, G.: The EarthCARE Satellite: the next step forward in global measurements of clouds, aerosols, precipitation and radiation, B. Am. Meteorol. Soc, 96, 1311–1332,, 2015. a

Kato, S., Rose, F. G., and Charlock, T. P.: Computation of Domain-Averaged Irradiance Using Satellite-Derived Cloud Properties, J. Atmos. Ocean. Tech., 22, 146–164,, 2005. a

King, M. D. and Harshvardhan: Comparative Accuracy of Selected Multiple Scattering Approximations, J. Atmos. Sci., 43, 784–801,<0784:CAOSMS>2.0.CO;2, 1986. a

Li, J., Yi, Y., Minnis, P., Huang, J., Yan, H., Ma, Y., Wang, W., and Ayers,  J. K.: Radiative effect differences between multi-layered and single-layer clouds derived from CERES, CALIPSO, and CloudSat data, J. Quant. Spectrosc. Ra., 112, 361–375,, international Symposium on Atmospheric Light Scattering and Remote Sensing (ISALSaRS'09), 2011. a

Loeb, N. G. and Manalo-Smith, N.: Top-of-Atmosphere Direct Radiative Effect of Aerosols over Global Oceans from Merged CERES and MODIS Observations, J. Climate, 18, 3506–3526,, 2005. a

Loeb, N. G., Manalo-Smith, N., Kato, S., Miller, W. F., Gupta, S. K., Minnis,  P., and Wielicki, B. A.: Angular Distribution Models for Top-of-Atmosphere Radiative Flux Estimation from the Clouds and the Earth's Radiant Energy System Instrument on the Tropical Rainfall Measuring Mission Satellite. Part I: Methodology, J. Appl. Meteorol., 42, 240–265,<0240:ADMFTO>2.0.CO;2, 2003. a, b

Loeb, N. G., Kato, S., Loukachine, K., and Manalo-Smith, N.: Angular Distribution Models for Top-of-Atmosphere Radiative Flux Estimation from the Clouds and the Earth's Radiant Energy System Instrument on the Terra Satellite. Part I: Methodology, J. Atmos. Ocean. Tech., 22, 338–351,, 2005. a, b

Loeb, N. G., Doelling, D. R., Wang, H., Su, W., Nguyen, C., Corbett, J. G., Liang, L., Mitrescu, C., Rose, F. G., and Kato, S.: Clouds and the Earth's Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) Top-of-Atmosphere (TOA) Edition-4.0 Data Product, J. Climate, 31, 895–918,, 2018. a

Meador, W. E. and Weaver, W. R.: Two-Stream Approximations to Radiative Transfer in Planetary Atmospheres: A Unified Description of Existing Methods and a New Improvement, J. Atmos. Sci., 37, 630–643,<0630:TSATRT>2.0.CO;2, 1980. a

Nam, C., Bony, S., Dufresne, J.-L., and Chepfer, H.: The “too few, too bright” tropical low-cloud problem in CMIP5 models, Geophys. Res. Lett., 39, L21801,, 2012. a

NASA Langley Research Center: CERES Data Products, available at:, last access: 17 July 2020. a

Poli, P., Hersbach, H., Dee, D. P., Berrisford, P., Simmons, A. J., Vitart, F., Laloyaux, P., Tan, D. G. H., Peubey, C., Thépaut, J.-N., Trémolet, Y., Hólm, E. V., Bonavita, M., Isaksen, L., and Fisher, M.: ERA-20C: An Atmospheric Reanalysis of the Twentieth Century, J. Climate, 29, 4083–4097,, 2016. a

Schaaf, C.: MODIS/Terra+Aqua BRDF/Albedo Gap-Filled Snow-Free Daily L3 Global 30ArcSec CMG V006, distributed by NASA EOSDIS Land Processes DAAC,, 2019. a

Scheck, L., Weissmann, M., and Mayer, B.: Efficient Methods to Account for Cloud-Top Inclination and Cloud Overlap in Synthetic Visible Satellite Images, J. Atmos. Ocean. Tech., 35, 665–685,, 2018. a

Shettle, E. P. and Weinman, J. A.: The Transfer of Solar Irradiance Through Inhomogeneous Turbid Atmospheres Evaluated by Eddington's Approximation,  J. Atmos. Sci., 27, 1048–1055,<1048:TTOSIT>2.0.CO;2, 1970. a

Smith, G. L., Green, R. N., Raschke, E., Avis, L. M., Suttles, J. T., Wielicki,  B. A., and Davies, R.: Inversion methods for satellite studies of the Earth's Radiation Budget: Development of algorithms for the ERBE Mission, Rev. Geophys., 24, 407–421,, 1986. a

Stephens, G., Li, J., and Wild, M.: An update on Earth's energy balance in light of the latest global observations, Nat. Geosci., 5, 691–696,, 2012. a

Su, W., Corbett, J., Eitzen, Z., and Liang, L.: Next-generation angular distribution models for top-of-atmosphere radiative flux calculation from CERES instruments: methodology, Atmos. Meas. Tech., 8, 611–632,, 2015. a, b, c

Sun-Mack, S., Minnis, P., Chen, Y., Doelling, D. R., Scarino, B. R., Haney,  C. O., and Smith, W. L. J.: Calibration Changes to Terra MODIS Collection-5 Radiances for CERES Edition 4 Cloud Retrievals, IEEE T. Geosci. Remote, 56, 6016–6032,, 2018. a

Suttles, J. T., Green, R. N., Minnis, P., Smith, G., Staylor, W. F., Wielicki,  B., J., W. I., Young, D., Taylor, V. R., and Stowe, L. L.: Angular radiation models for Earth-atmosphere systems, Volume 1 – Shortwave Radiation, Tech. Rep. 1184, NASA, Washington, D.C., USA, 1988. a

Thorsen, T. J., Kato, S., Loeb, N. G., and Rose, F. G.: Observation-Based Decomposition of Radiative Perturbations and Radiative Kernels, J. Climate, 31, 10039–10058,, 2018. a

Tornow, F., Preusker, R., Domenech, C., Carbajal Henken, C. K., Testorp, S., and Fischer, J.: Top-of-Atmosphere Shortwave Anisotropy over Liquid Clouds: Sensitivity to Clouds' Microphysical Structure and Cloud-Topped Moisture, Atmosphere, 9, 256,, 2018. a, b

Tornow, F., Domenech, C., and Fischer, J.: On the Use of Geophysical Parameters for the Top-of-Atmosphere Shortwave Clear-Sky Radiance-to-Flux Conversion in EarthCARE, J. Atmos. Ocean. Tech., 36, 717–732,, 2019. a

van de Hulst, H.: Scaling Laws in Multiple Light Scattering under very Small Angles. (Karl Schwarzschild Lecture 1995), Rev. Mod. Astron., 9, 1–16, 1996. a

Wald, L. and Monget, J. M.: Sea surface winds from sun glitter observations,  J. Geophys. Res.-Oceans, 88, 2547–2555,, 1983. a

Wielicki, B. A., Barkstrom, B. R., Harrison, E. F., III, R. B. L., Smith,  G. L., and Cooper, J. E.: Clouds and the Earth's Radiant Energy System (CERES): An Earth Observing System Experiment, B. Am. Meteorol. Soc., 77, 853–868,<0853:CATERE>2.0.CO;2, 1996.  a

Zhang, Y., Macke, A., and Albers, F.: Effect of crystal size spectrum and crystal shape on stratiform cirrus radiative forcing, Atmos. Res., 52, 59–75,, 1999. a

Short summary
Clouds reflect sunlight unevenly, which makes it difficult to quantify the portion reflected back to space via satellite observation. To improve quantification, we propose a new statistical model that incorporates more satellite-inferred cloud and atmospheric properties than state-of-the-art models. We use concepts from radiative transfer theory that we statistically optimize to fit observations. The new model often explains past satellite observations better and predicts reflection plausibly.