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 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-Smith, 2005; 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 Wehr, 2012). 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 $\stackrel{\mathrm{̃}}{\mathit{\tau }}\approx \mathrm{10}$ and identified additional parameters that influence ADMs: above-cloud water vapor (ACWV) and layer mean cloud-top effective radius $\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}$. They showed that ignoring these parameters could cause deviations in instantaneous flux estimates of about 10 Wm⁻². 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 $\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}$ 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 θ₀, 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^{*}\frac{S_{\mathrm{0}}\cos {\mathit{\theta }}_{\mathrm{0}}}{F^{\downarrow }}$, with solar constant S₀=1361.0 Wm⁻².

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 $\stackrel{\mathrm{̃}}{\mathit{\tau }}=e^{\stackrel{\mathrm{‾}}{\log \mathit{\tau }}}$, layer average $\overline{\mathit{\tau }}$, 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 $\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}$ (using band 3.7 µm), and cloud-top pressure p^ctop. 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 w_10 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 $e_{\mathrm{s}}=\mathrm{6.112}{e}^{\frac{\mathrm{17.67}T}{T+\mathrm{243.5}}}$ (using T in ^∘C) (Bolton, 1980), gravitational acceleration g, and molecular weights of water and dry air mol_h2o and mol_air respectively:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{ACWV}={\displaystyle \frac{\mathrm{1}}{g}}\underset{p^{\text{ctop}}}{\overset{\mathrm{0}}{\int }}\text{mr}(p,T,\text{RH})\mathrm{d}p\\ \text{(1)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{g}}\underset{p^{\text{ctop}}}{\overset{\mathrm{0}}{\int }}{\displaystyle \frac{e_{\mathrm{s}}\left(T\right)}{p}}{\displaystyle \frac{{\text{mol}}_{h\mathrm{20}}}{{\text{mol}}_{\text{air}}}}\text{RH}\left(p\right)\mathrm{d}p.\end{array}$$

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 1 Number of CERES footprints obtained after screening for marine clouds. Number are shown in millions (in total 1 711 937 663 footprints).

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 ${\mathit{\theta }}_{\mathrm{0}}^{\mathrm{\Delta }}$, ${\mathit{\theta }}_{\mathrm{v}}^{\mathrm{\Delta }}$, and φ^Δ, respectively. Combinations of ${\mathit{\theta }}_{\mathrm{0}}^{\mathrm{\Delta }}$, ${\mathit{\theta }}_{\mathrm{v}}^{\mathrm{\Delta }}$, 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 $\stackrel{\mathrm{̃}}{\mathit{\tau }}$.

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }I\left({\mathit{\theta }}_{\mathrm{0}}^{\mathrm{\Delta }},{\mathit{\theta }}_{\mathrm{v}}^{\mathrm{\Delta }},{\mathit{\phi }}^{\mathrm{\Delta }}\right)=I_{\mathrm{0}}+{\displaystyle \frac{a}{{\left[\mathrm{1}+e^{-\frac{(x-x_{\mathrm{0}})}{b}}\right]}^c}},\end{array}$$

where $x=\log f\stackrel{\mathrm{̃}}{\mathit{\tau }}$ for a single cloud layer or $x=\log \left[\right(f_{\mathrm{1}}+f_{\mathrm{2}}\left)e^{\frac{f_{\mathrm{1}}\log {\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathrm{1}}+f_{\mathrm{2}}\log {\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathrm{2}}}{f_{\mathrm{1}}+f_{\mathrm{2}}}}\right]$ for two layers, and I₀, a, b, c, and x₀ 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 1 For an example angular bin (${\mathit{\theta }}_{\mathrm{0}}\in [\mathrm{20}{}^{\circ },\mathrm{22}{}^{\circ }]$, ${\mathit{\theta }}_{\mathrm{v}}\in [\mathrm{6}{}^{\circ },\mathrm{8}{}^{\circ }]$, $\mathit{\phi }\in [\mathrm{12}{}^{\circ },\mathrm{14}{}^{\circ }]$), 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 f₁ and f₂ 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 ${\mathit{\varphi }}_{\text{eff}}=\frac{f_{\mathrm{1}}{\mathit{\varphi }}_{\mathrm{1}}+f_{\mathrm{2}}{\mathit{\varphi }}_{\mathrm{2}}}{f_{\mathrm{1}}+f_{\mathrm{2}}}$ for two layers, and the following thresholds: liquid for $\mathrm{1}<{\mathit{\varphi }}_{\text{eff}}<\mathrm{1.01}$, mixed for $\mathrm{1.01}\le {\mathit{\varphi }}_{\text{eff}}\le \mathrm{1.75}$, and ice for $\mathrm{1.75}<{\mathit{\varphi }}_{\text{eff}}\le \mathrm{2}$.

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 >10 m s⁻¹) 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 $\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}$ and ACWV into a radiance-predicting statistical model. One could divide each angular bin's samples into classes of $\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}$ 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:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }I\left({\mathit{\theta }}_{\mathrm{0}}^{\mathrm{\Delta }},{\mathit{\theta }}_{\mathrm{v}}^{\mathrm{\Delta }},{\mathit{\phi }}^{\mathrm{\Delta }}\right)\approx S_o\cos {\mathit{\theta }}_{\mathrm{0}}\mathit{\alpha }{e}^{-\mathrm{2}\text{ACWV}},\end{array}$$

with solar influx S_ocos θ₀ 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.

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }\log I\left({\mathit{\theta }}_{\mathrm{0}}^{\mathrm{\Delta }},{\mathit{\theta }}_{\mathrm{v}}^{\mathrm{\Delta }},{\mathit{\phi }}^{\mathrm{\Delta }}\right)\approx A+B\log \mathit{\alpha }+C\text{ACWV}\end{array}$$

Like the state-of-the-art methodology (Sect. 3.1), we applied this approach per angular bin (resolved by 2^∘ in θ₀, θ_v, and φ), allowing us to treat S_ocos θ₀ 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 ${\mathit{\varphi }}_{\mathrm{1}/\mathrm{2}}=\mathrm{1}$, as ice for ${\mathit{\varphi }}_{\mathrm{1}/\mathrm{2}}=\mathrm{2}$, and as mixed for ϕ₁=1 and ϕ₂=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, w_10 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 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:

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }{r}^{\text{SunGlint}}={\displaystyle \frac{\mathit{\pi }\mathit{\rho }\left(\mathit{\omega }\right)P({\mathit{\theta }}_n,W_{\mathrm{10}\mathrm{m}})}{\mathrm{4}\cos {\mathit{\theta }}_{\mathrm{0}}\cos {\mathit{\theta }}_{\mathrm{v}}{\cos }^{\mathrm{4}}{\mathit{\theta }}_n}},\end{array}$$

where

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }P({\mathit{\theta }}_n,W_{\mathrm{10}\mathrm{m}})={\displaystyle \frac{\mathrm{1}}{\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{{\tan }^{\mathrm{2}}{\mathit{\theta }}_n}{{\mathit{\sigma }}^{\mathrm{2}}}}\right),\text{(7)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}=\mathrm{0.003}+\mathrm{0.00512}{W}_{\mathrm{10}\mathrm{m}},\text{(8)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_n=\arccos \left({\displaystyle \frac{\cos {\mathit{\theta }}_{\mathrm{v}}+\cos {\mathit{\theta }}_{\mathrm{0}}}{\mathrm{2}\cos \mathit{\omega }}}\right),\text{(9)}& {\displaystyle }& {\displaystyle }\cos \mathrm{2}\mathit{\omega }=\cos {\mathit{\theta }}_{\mathrm{v}}\cos {\mathit{\theta }}_{\mathrm{0}}+\sin {\mathit{\theta }}_{\mathrm{v}}\sin {\mathit{\theta }}_{\mathrm{0}}\cos \mathit{\phi }.\end{array}$$

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 $\stackrel{\mathrm{̃}}{\mathit{\tau }}$ and cloud-top effective radius $\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}$ through asymmetry parameter $g\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ or backscattering fraction $\mathit{\beta }\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$, 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 Harshvardhan, 1986): the Eddington approximation (Shettle and Weinman, 1970) and the Coakley–Chylek approximation (using solution I of Coakley 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 μ₀=cos θ₀:

$$\begin{array}{}\text{(10)}& {\displaystyle }{\displaystyle }{\mathit{\alpha }}^{\text{TwoStream}}={\displaystyle \frac{\frac{\mathrm{3}}{\mathrm{4}}(\mathrm{1}-g)\stackrel{\mathrm{̃}}{\mathit{\tau }}-\frac{\mathrm{1}}{\mathrm{4}}(\mathrm{1}-\mathrm{3}{\mathit{\mu }}_{\mathrm{0}})\left(\mathrm{1}-e^{-\frac{\stackrel{\mathrm{̃}}{\mathit{\tau }}}{{\mathit{\mu }}_{\mathrm{0}}}}\right)}{\mathrm{1}+\frac{\mathrm{3}}{\mathrm{4}}(\mathrm{1}-g)\stackrel{\mathrm{̃}}{\mathit{\tau }}}}.\end{array}$$

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

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }{\mathit{\alpha }}^{\text{TwoStream}}={\displaystyle \frac{\frac{(\mathrm{1}-g)\stackrel{\mathrm{̃}}{\mathit{\tau }}}{\mathrm{2}}}{\mathrm{1}+\frac{(\mathrm{1}-g)\stackrel{\mathrm{̃}}{\mathit{\tau }}}{\mathrm{2}}}},\end{array}$$

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

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{\mathit{\alpha }}^{\text{TwoStream}}={\displaystyle \frac{{\mathit{\alpha }}^{\text{ocean}}+\frac{(\mathrm{1}-{\mathit{\alpha }}^{\text{ocean}})(\mathrm{1}-g)\stackrel{\mathrm{̃}}{\mathit{\tau }}}{\mathrm{2}}}{\mathrm{1}+\frac{(\mathrm{1}-{\mathit{\alpha }}^{\text{ocean}})(\mathrm{1}-g)\stackrel{\mathrm{̃}}{\mathit{\tau }}}{\mathrm{2}}}}.\end{array}$$

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 f₀ and cloud fractions of layer 1 and layer 2, f₁ and f₂, respectively:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\mathit{\alpha }& =f_{\mathrm{0}}\left({\mathit{\alpha }}^{\text{ocean}}+r^{\text{SunGlint}}\right)+f_{\mathrm{1}}{\mathit{\alpha }}_{\mathrm{1}}^{\text{TwoStream}}\\ & +f_{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{2}}^{\text{TwoStream}},\end{array}\end{array}$$

where $f_{\mathrm{0}}+f_{\mathrm{1}}+f_{\mathrm{2}}=\mathrm{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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ 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 2 For 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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ (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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ 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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ 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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ to choose. Initial tests that used a $g\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ 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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ for each angular bin and for each cloud phase (liquid and ice). Inspired by the shape of Mie-calculated $g\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$, we approximated $g\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ via a quadratic function,

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }g\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)=a+b\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}+c{\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}}^{\mathrm{2}},\end{array}$$

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 2 In search for optimal $g\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$, we list the range (minimum and maximum) and step size for each parameter in Eq. (14).

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 $\mathit{\nu }=\frac{{\overline{\mathit{\tau }}}^{\mathrm{2}}}{\mathit{\sigma }(\mathit{\tau }{)}^{\mathrm{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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ 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 3 Per solar geometry θ₀ 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.

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. 10–12) 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 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{sr}}^{-\mathrm{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 3 Applied to angular bins of the principal plane for ${\mathit{\theta }}_{\mathrm{0}}\in [\mathrm{20}{}^{\circ },\mathrm{22}{}^{\circ }]$, we test a variety of two-stream solutions for cloud albedo (Eqs. 10–12) 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 θ₀∼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 ${\mathit{\theta }}_{\mathrm{0}}>\mathrm{40}{}^{\circ }$, bins in sun-glint-affected geometries (constituting a portion of all bins in a hemisphere between 10 % for a ${\mathit{\theta }}_{\mathrm{0}}\sim \mathrm{20}{}^{\circ }$ and 1 % for a ${\mathit{\theta }}_{\mathrm{0}}\sim \mathrm{75}{}^{\circ }$) 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 ${\mathit{\theta }}_{\mathrm{0}}<\mathrm{20}{}^{\circ }$, 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 4 Using all observed angular bins within ${\mathit{\theta }}_{\mathrm{0}}\in [\mathrm{6}{}^{\circ },\mathrm{82}{}^{\circ }]$, 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: $\mathit{\delta }=\left[\mathit{\sigma }\right(\mathrm{\Delta }{I}^{\text{LogLinear}})-\mathit{\sigma }(\mathrm{\Delta }{I}^{\text{Sigmoidal}}\left)\right]/\mathit{\sigma }\left(\mathrm{\Delta }{I}^{\text{Sigmoidal}}\right)\cdot \mathrm{100}\mathrm{\%}$. 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 ($f_{\mathrm{1}}=f_{\mathrm{2}}=\mathrm{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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ 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 5 For angular bins along the principal plane for ${\mathit{\theta }}_{\mathrm{0}}\in [\mathrm{20}{}^{\circ },\mathrm{22}{}^{\circ }]$ 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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ for three R_e (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⁻²; 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 (f_1∕2, ${\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathrm{1}/\mathrm{2}}$) but also for microphysical (${\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}}_{\mathrm{1}/\mathrm{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 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{sr}}^{-\mathrm{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 θ₀∼20 and 70^∘, we found model uncertainty of most bins reduced by as much as 35.8 %.

Figure 6 For angular bins along the principal plane for ${\mathit{\theta }}_{\mathrm{0}}\in [\mathrm{20}{}^{\circ },\mathrm{22}{}^{\circ }]$, 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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right)$ 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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right|{\mathit{\theta }}_{\mathrm{0}}^{\mathrm{\Delta }},{\mathit{\theta }}_{\mathrm{v}}^{\mathrm{\Delta }},{\mathit{\phi }}^{\mathrm{\Delta }})$, 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\left(\stackrel{\mathrm{‾}}{R_{\mathrm{e}}}\right|{\mathit{\theta }}_{\mathrm{0}}^{\mathrm{\Delta }},{\mathit{\theta }}_{\mathrm{v}}^{\mathrm{\Delta }},{\mathit{\phi }}^{\mathrm{\Delta }})$, 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, $(\mathrm{1}-g)\stackrel{\mathrm{̃}}{\mathit{\tau }}$, 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 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; https://doi.org/10.5067/MODIS/MCD43GF.006, Schaaf, 2019), 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.