Influence of cloud retrieval errors due to three-dimensional radiative effects on calculations of broadband shortwave cloud radiative effect

We investigate how cloud retrieval errors due to the three-dimensional (3D) radiative effects affect broadband shortwave (SW) cloud radiative effects (CREs) in shallow cumulus clouds. A framework based on the combination of large eddy simulations (LESs) and radiative transfer (RT) models was developed to simulate both one-dimensional (1D) and 3D radiance, as well as SW broadband fluxes. Results show that the broadband SW fluxes reflected at top of the domain, transmitted at the surface, and absorbed in the atmosphere, computed from the cloud retrievals using 1D RT ($F_{\mathrm{1}\mathrm{D}}^{\ast }$), can provide reasonable broadband radiative energy estimates in comparison with those derived from the true cloud fields using 1D RT (F_1D). The difference between these 1D-RT-simulated fluxes ($F_{\mathrm{1}\mathrm{D}}^{\ast }$, F_1D) and the benchmark 3D RT simulations computed from the true cloud field (F_3D) depends primarily on the horizontal transport of photons in 3D RT, whose characteristics vary with the sun's geometry. When the solar zenith angle (SZA) is 5°, the domain-averaged $F_{\mathrm{1}\mathrm{D}}^{\ast }$ values are in excellent agreement with the F_3D, all within 7 % relative CRE bias. When the SZA is 60°, the CRE differences between calculations from $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and F_3D are determined by how the cloud side-brightening and darkening effects offset each other in the radiance, retrieval, and broadband fluxes. This study suggests that although the cloud property retrievals based on the 1D RT theory may be biased due to the 3D radiative effects, they still provide CRE estimates that are comparable to or better than CREs calculated from the true cloud properties using 1D RT.

1 Introduction

Covering about 60 %–70 % of the Earth's surface (Rossow and Schiffer, 1999; Vardavas and Taylor, 2011), clouds play a very important role in the Earth's climate system. Clouds can cool the Earth by reflecting shortwave (SW) solar radiative flux back to space and at the same time warm the Earth by retaining the outgoing longwave (LW) infrared radiative flux at the top of the atmosphere (TOA), known as cloud radiative effects (CREs). The annual global average TOA CRE is approximately −50 W m⁻² at SW and 30 W m⁻² at LW, resulting in a net CRE of about −20 W m⁻² (Stocker et al., 2013). These strong CREs show that clouds greatly affect the Earth's energy budget (Ramanathan et al., 1989; Kiehl and Trenberth, 1997; Trenberth et al., 2009). The CRE of clouds is largely determined by the optical and microphysical properties of clouds including the cloud optical thickness (τ), cloud droplet effective radius (r_e), and cloud liquid water path (LWP). Thus, continuous measurements of these cloud properties from regional to global scales are critical to better understand the role of clouds in the climate systems. Currently, satellite-based remote sensing is the only way to make such observations. Remotely “retrieved” cloud properties based on these satellite observations are often used to derive the radiative effects of clouds (e.g., Wielicki et al., 1996; Platnick et al., 2003; Loeb and Manalo-Smith, 2005; Oreopoulos et al., 2016) and evaluate the simulations of Earth system models (ESMs) (Kay et al., 2012; Nam et al., 2012; Song et al., 2018).

A commonly used retrieval technique in passive satellite remote sensing is the bispectral retrieval method first developed by Nakajima and King (1990). It retrieves τ and r_e simultaneously from a pair of total reflectance measurements, one from the non-absorbing visible or near-infrared (VNIR) band (e.g., 0.66 µm) and the other from the moderately absorbing shortwave infrared (SWIR) band (e.g., 2.13 µm). Since clouds in reality have three-dimensional (3D) structures, the simulation of radiative transfer (RT) in clouds should ideally consider the transport of radiation in both vertical and horizontal directions (referred to as “3D RT”). Unfortunately, the computational cost for 3D RT is extremely high. As a result, operational bispectral cloud retrievals are almost exclusively based on the one-dimensional (1D) RT theory that considers only the vertical and ignores the net horizontal transport of radiation. The radiative properties of clouds under 3D RT are substantially different from those under 1D RT. This is known as the 3D radiative effects and can lead to substantial biases in cloud property retrievals based on 1D RT (Várnai and Marshak, 2002; Marshak et al., 2006; Zhang et al., 2012, 2016). Although recent efforts have been made to employ machine learning techniques to retrieve cloud properties based on 3D RT theory (Okamura et al., 2017; Masuda et al., 2019; Nataraja et al., 2022), these machine-learning-based algorithms are still in their infancy and far from being used in operational algorithms.

Many previous studies have investigated the 3D radiative effects on satellite radiance observations and cloud property retrievals. For example, Welch and Wielicki (1984) used some “toy” cloud fields (e.g., cubic and cylindrical) to illustrate the impact of side-illuminating and mutual shadowing on cloud albedo. Várnai and Davis (1999) and Várnai (2000) elucidated several 3D RT mechanisms, e.g., upward–downward trapping–escaping, that can result in significant differences between 1D and 3D cloud albedo. Hogan et al. (2019) proposed a distinct mechanism, named “entrapment”, which plays a key role in the 3D radiative effect of clouds. Davis and Marshak (2001) pointed out that the channeling effect in 3D RT can smoothen out the small-scale cloud variations and lead to the reduction of cloud brightness at cloud edges. Marshak et al. (2006) explained how the radiance biases due to 3D radiative effects can lead to τ and r_e retrieval biases in MODIS (Moderate Resolution Imaging Spectroradiometer) cloud products. This study is built upon these classic papers but has a different objective.

Here, we investigate an important question: do cloud property retrievals based on 1D RT, which are potentially biased due to the 3D radiative effects, still provide an observational basis to estimate the broadband SW CREs? This is an important question because as mentioned above, operational cloud retrieval products from, for example, MODIS, are frequently used for CRE estimation and ESM evaluation. However, to our best knowledge, the impacts of retrieval bias due to the 3D radiative effects on such applications have never been examined systematically in previous studies. To better explain our objective and the difference of this study from many previous ones on the 3D radiative effects, we need to introduce the framework illustrated in Fig. 1. As conceptually illustrated in Fig. 1, the observed radiances are inherently 3D (i.e., from Box A to C) because the RT in nature is 3D. However, when 1D RT theory is used to interpret the observations, we get the “retrieved cloud properties” in Box D that can be significantly different from the “true” cloud properties in Box A. Although the retrieved cloud properties are often biased due to the 3D radiative effects, they are still widely used to compute the radiative fluxes by clouds (i.e., from Box D to E) using 1D RT and the results are often used for studying the climatic effects of clouds (e.g., Kato et al., 2011; Zelinka et al., 2012; Oreopoulos et al., 2016). In contrast, the true radiative fluxes in nature are also 3D (i.e., from Box A to F). A few recent studies have computed and compared the 1D and 3D radiative fluxes as well as heating rates by clouds. For example, Barker et al. (2011, 2012) and Okata et al. (2017) compared the 1D and 3D SW fluxes computed based on the constructed A-Train cloud scenes at the TOA and surface. A more recent study by Singer et al. (2021) utilized large eddy simulation (LES) cloud fields of different cloud regimes to assess the SW radiative flux and TOA albedo bias associated with the 3D effects. The main difference between their study and this current work is as follows: they compared the 3D (i.e., Box F in Fig. 1) with the 1D broadband fluxes (i.e., Box G in Fig. 1) both computed from the true clouds. In contrast, we argue that the true clouds are not known in practice and therefore we compare the 3D flux (i.e., Box F in Fig. 1) with the 1D flux computed from the retrieved cloud properties (i.e., Box E in Fig. 1); this approach enables us to measure the impact of cloud retrieval errors on the radiative flux and CRE.

Figure 1 Conceptual framework to understand the study. R_3D and R_1D are the reflectances from three-dimensional (3D) and one-dimensional (1D) radiative transfer (RT), respectively, while δR is their difference. X represents the true cloud field and X^∗(R_3D) represents the retrieved cloud properties from 3D RT reflectance, while δX is the cloud property retrieval bias. $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and F_1D are the radiative flux calculated using 1D RT on the retrieved cloud properties and true cloud properties, respectively. F_3D is the radiative flux derived from the true cloud field using 3D RT. δF₁ and δF₂ are the differences between the pairs ($F_{\mathrm{1}\mathrm{D}}^{\ast }$, F_3D) and (F_1D, F_3D), respectively.

To determine whether biased cloud retrievals of cloud properties can still provide an observational basis for CREs, we focus on three important scientific questions (SQs) as illustrated in Fig. 1.

• SQ 1: how does the reflectance simulated based on 3D RT (R_3D) compare with the reflectance simulated based on 1D RT (R_1D) for different types of clouds at different illuminating–viewing geometries (i.e., comparing Box C to B in Fig. 1)?

• SQ 2: how do the retrieved cloud properties, e.g., cloud optical thickness and cloud droplet effective radius derived from the 3D reflectance using 1D RT, compare to the true cloud properties (i.e., comparing Box D to A in Fig. 1)?

• SQ 3: comparing δF₁ to δF₂ in Fig. 1, i.e., how are the broadband SW radiative fluxes derived from the retrieved cloud properties using 1D RT ($F_{\mathrm{1}\mathrm{D}}^{\ast }$) (see Box E in Fig. 1) different from the true radiative fluxes computed from the true cloud fields using 3D RT (F_3D) (see Box F in Fig. 1)? And how does this result compare with the difference between F_3D and the broadband SW radiative fluxes computed from the true cloud properties using 1D RT (F_1D) (see Box G in Fig. 1)?

The paper's remaining structure is arranged as follows: Sect. 2 briefly describes the data and theory for the study. Section 3 presents and discusses results on how the 3D radiative effects influence the radiance fields, cloud property retrievals, and broadband radiative flux. The summary and conclusion are given in Sect. 4.

2 Data and theory

2.1 Cloud field data set

A great challenge facing 3D radiative effects studies is that the true clouds are always obscured by the 3D radiative effects, which are inevitable in real observations. To overcome this challenge, many previous studies (e.g., Zhang et al., 2012; Miller et al., 2018; Rajapakshe and Zhang, 2020) have used synthetic cloud fields and RT simulations to mimic the observation–retrieval process and study the 3D radiative effects. Building on these previous studies, we adopt the same state-of-the-art satellite retrieval simulator by Zhang et al. (2012) and add a broadband flux computation function to study the 3D radiative effect and its impact on the broadband SW radiative flux. As described in Zhang et al. (2012) and illustrated in Fig. 1, the framework consists of three major components: (1) synthetic cloud fields, (2) RT models (for radiance and broadband flux simulations), and (3) a cloud property (e.g., τ and r_e) retrieval simulator. LES cloud fields, which are commonly used in different cloud microphysical and 3D effect studies (e.g., Singer et al., 2021; Zhang et al., 2012), are based on computational models and mathematical equations to simulate atmospheric behavior and get the 3D cloud property. Certain studies (e.g., Levis et al., 2015; Loveridge et al., 2023) have developed atmospheric tomography techniques to reconstruct 3D cloud scenes from observational data but have yet to be widely used globally. Similar to Zhang et al. (2012), the synthetic cloud fields utilized in this study are based on LES cloud fields.

Since the 3D radiative effects on overcast clouds are minimal, two cloud fields of low and intermediate cloud fractions have been selected as a case study to illustrate the framework explained in Sect. 1. The selected cloud fields are from the LES Atmospheric Radiation Measurement (ARM) Symbiotic Simulation and Observation (LASSO) activity, conducted in the ARM Southern Great Plains (SGP) site located in Lamont, Oklahoma (Gustafson et al., 2020) (https://www.arm.gov/capabilities/modeling/lasso/, last access: 19 May 2023). LASSO enhances ARM's observations by using LES modeling to provide contextual and self-consistent representation of the atmosphere surrounding the ARM site. It also provides continuous observations from ground-based cloud and radiometric instruments, which is valuable for enhancing research on cloud–radiation interactions. For this study, the two snapshots of LASSO LES cloud field cases analyzed are 14:00 UTC on 27 June 2015 (simulation ID 108) and the other at 14:00 UTC on 18 August 2016 (simulation ID 113). For conciseness in this text, these snapshots will be referred to as “27 June” and “18 August”, respectively. We chose to use these specific LASSO LES cloud field data from the stated dates because they represent typical shallow cumulus clouds, do not contain ice (to avoid the complexities dealing with ice microphysics), and have better diagnostic statistics compared to other LES data streams. It is important to note that, because the impact of 3D radiative effects varies substantially for different cloud regimes and surface types, this study is constrained to shallow cumulus cloud types (over land surface) found at the LASSO SGP site.

Figure 2 Large eddy simulation (LES) of cloud liquid water path (LWP) for 14:00 UTC on 27 June 2015 (a) and 14:00 UTC on 18 August 2016 (b) at the ARM SGP atmospheric observatory. White areas are clear-sky regions where the cloud liquid water path (LWP) is 0.

The LASSO LES cloud fields for this study are characterized by broken cloud patterns spatially distributed across the domain as seen in the LWP maps in Fig. 2a and b for the 27 June and 18 August cases, respectively. The 3D distribution of cloud liquid water content (LWC) was obtained from the LASSO cloud field data, and a two-moment bulk microphysics scheme by Morrison and Gettelman (2008) (see their Eq. 5 in Sect. 2) was used to obtain the r_e associated with the corresponding LWC distribution. It is important to note that for this study, a cloudy column has been defined as a column with LWP > 0 (i.e., clear-sky regions have LWP = 0). The cloud fields have different domain sizes and microphysics distribution, and the cloud cover for the 18 August cloud field (47.08 %) is more than twice that of the 27 June cloud field (20.15 %). Information about the cloud properties and the LES domain is summarized in Table 1.

Table 1 Cloud property characteristics for the LES cloud field cases. The mean cloud effective radius (r_e), mean cloud optical thickness (τ), and in-cloud liquid water path are from the average of the cloudy regions only. The columns from left to right are the case name, cloud fraction, mean in-cloud liquid water path, mean cloud-base height (CBH), mean cloud-top height (CTH), mean r_e, mean τ, grid spacing, and domain size, respectively.

2.2 Radiative transfer setup

We use the spherical harmonics discrete ordinate method (SHDOM) RT model developed by Evans (1998) to handle both 1D and 3D radiance computations. We have benchmarked the SHDOM simulations against the results from our previous studies (Zhang et al., 2012; Miller et al., 2016). Broadband SW radiative flux computations, both 1D and 3D, were performed with the Intercomparison of 3D Radiation Codes (I3RC) Monte Carlo community model (Pincus and Evans, 2009), and atmospheric gaseous absorption was incorporated via the SW Rapid Radiative Transfer Model (RRTM) correlated k-distribution approach (Mlawer et al., 1997), which consists of 14 bands with a spectral range from 0.2 to 12 µm (this coupled broadband radiative flux solver is hereafter known as the I3RC+CKD model). Rayleigh scattering was included in the flux RT calculations, the background atmospheric profiles are taken to be horizontally homogeneous throughout the domain, and the profiles of atmospheric temperature, pressure, ozone, air density, and water vapor utilized for the RT flux calculations were obtained from the sounding data at the ARM SGP site on 27 June 2015. Several studies (e.g., Gristey et al., 2022) have shown that aerosol embedded in clouds with small aspect ratios (similar to our chosen LASSO LES cloud fields) has a significant influence on the 3D radiative effect. Thus, for simplicity in our study, ambient aerosols are neglected in the RT calculations. The 1D broadband RT flux calculations were performed with the same I3RC+CKD model by dividing the LES domain into individual columns, and RT was calculated on each LES column property separately and independently.

The spectral cloud optical properties were calculated using Mie scattering theory and were averaged over each of the RRTM spectral bands. The phase functions were represented using Legendre coefficients with 35 log-spaced effective radii spanning from 2 to 40 µm. The surface was assumed to be Lambertian with surface spectral albedos obtained from the ARM SGP site (see Fig. 4 in Coddington et al., 2013) applied for wavelength (λ) in the range 0.2 $\le \mathit{\lambda }\le$ 2.5 µm, while surface spectral albedo corresponding to a vegetation-covered surface (Zhuravleva et al., 2009) was utilized for λ> 2.5 µm (see Appendix B for the surface spectral albedo plot used in this study). In the Monte Carlo calculations, 10⁸ and 10⁴ photons were initiated for calculations of the 3D broadband SW flux and the column-independent 1D broadband SW flux, respectively. The radiative transfer calculations were implemented for two solar zenith angles (SZAs), a high-sun case with SZA of 5° and a low-sun case with SZA of 60°. In the broadband flux calculations, the downward flux at the top of the domain (TOD) corresponds to 1363 and 684.1 W m⁻² for SZA 5 and 60°, respectively. Throughout this study, we choose a constant 0° solar azimuth angle (SAA) and a constant 0° viewing zenith angle (VZA). Double periodic horizontal boundary conditions were applied for all the RT calculations, and all RT calculations have been conducted at the native LES resolution of 100 m. Current satellite remote sensing instruments have different footprints (e.g., 1 km footprint for MODIS), which can have different 3D effect signatures on the retrievals and impact the derived radiative flux. Therefore, future studies will investigate how 3D effect retrieval errors for different spatial resolutions (coarse and fine) affect the radiative flux estimates.

2.3 Bispectral retrieval method

The bispectral retrieval method introduced in Sect. 1 is solely based on the 1D RT theory to interpret the observed cloud reflectance. It is implemented using a precomputed lookup table (LUT), which consists of 1D reflectance functions for different τ and r_e combinations at the required solar viewing geometry (an example LUT is shown in Fig. 3). The observed cloud reflectance is then utilized as input to the LUT to simultaneously retrieve the τ and r_e via a two-dimensional (2D) interpolation between the observed cloud reflectance and the LUT grid. Notably, in the bispectral LUT regions with smaller τ, the retrieval uncertainty increases because the isolines of the LUT τ are less orthogonal and more tightly packed.

Figure 3 An example Nakajima and King (1990) bispectral lookup table (LUT) space. The solid lines are the reflectance function contours for fixed cloud effective radius (r_e), while the dashed lines are for fixed cloud optical thickness (τ). The surface is Lambertian with a surface albedo of 0.07. The solar zenith angle (SZA) is 60°, the view zenith angle (VZA) is 0°, and the solar azimuth angle (SAA) is 0°.

This nonlinearity in the LUT has high inhomogeneity consequences for cloud retrievals at the pixel level (Zhang et al., 2012, 2016). In this study, the VNIR reflectances were measured at 0.66 µm (identical to the central wavelength of the operational MODIS retrieval algorithm over a vegetated land surface), while the SWIR reflectances were measured at the 2.13 µm wavelength. The LUTs utilized for our bispectral retrievals have 19 effective radii spanning from 5 to 40 µm and 43 log-spaced τ values spanning from 0.05 to 158.48, while a constant effective variance (v_e) value of 0.1 is used for consistency with all other RT simulations in this study. The surface albedo at both 0.66 and 2.13 µm wavelengths for the LES radiance simulations and LUT RT calculations was 0.07. This value is consistent with the surface albedo of similar spectral bands in the broadband SW flux computations (see the spectral albedo plot in Appendix B).

2.4 Classification of failed and successful retrievals

One major challenge in cloud property retrievals from satellite remote sensing instruments like MODIS is a so-called “failed retrieval”. A retrieval can be considered failed if there is no r_e and τ LUT grid combination to interpret the reflectance observation or if there are no realistic cloud microphysics to explain the retrieved cloud property (e.g., a retrieved r_e> 40 µm). These can be due to several factors, such as the limits of the LUT, cloud overlapping effect, presence of partially cloudy pixels, extreme solar-satellite viewing geometries, strategy used in cloud mask implementation, and the optical characteristics of the underlying surface. Potential causes and rates of occurrence of failed MODIS retrievals for marine liquid-phase clouds have been studied extensively (Cho et al., 2015). In this study, we refer to the MODIS cloud property retrieval algorithm's classification of failed retrievals (Platnick et al., 2016) and the study by Cho et al. (2015) to classify a pixel as a successful or failed retrieval as explained below.

• 1.

  For observations with both VNIR and SWIR reflectance observations within the LUT solution space, the nearest interpolated τ and r_e values are retrieved (pink area bounded by the LUT lines in Fig. 3). If the observed VNIR reflectances exceed the upper limit of LUT τ but within the LUT r_e solution range (extended pink area in Fig. 3), the nearest LUT r_e is retrieved and the maximum LUT τ value (τ=158.48) is assigned to the retrieval. These explained categories are classified as “successful retrievals” for this study.

• 2.

  In other cases, for observations with VNIR reflectance within the LUT solution space but SWIR reflectance above the LUT solution space (purple area in Fig. 3), the nearest τ values are retrieved but the smallest LUT r_e value of 5 µm is assigned to the retrievals. This category of retrieval failure is called “r_e too small” failures. In cases where the VNIR reflectance observations are within the LUT τ solution space but the SWIR reflectances are below the LUT solution space (green area in Fig. 3), the nearest τ values are retrieved but the largest LUT r_e value of 40 µm is assigned to the retrieval. This category of retrieval failure is called the “r_e too large” failures. In cases where the observed VNIR reflectance is greater than the largest LUT τ value and the observed SWIR reflectance is smaller than the largest LUT r_e (i.e., the lower yellow region in Fig. 3), the retrievals are assigned the largest τ value (τ=158.48) and the largest r_e value (r_e=40 µm). For observations with VNIR reflectance greater than the largest LUT τ value and the SWIR reflectance greater than the smallest LUT r_e value (i.e., the upper yellow region in Fig. 3), the retrievals are assigned the largest τ value (τ=158.48) and smallest r_e value (r_e=5 µm). Lastly, for observations with VNIR reflectance below the minimum LUT τ (red area in Fig. 3), the r_e and τ retrievals are assigned fill values (which are represented by τ=0 in our flux calculations). These explained categories are called τ failures. The “r_e too small”, “r_e too large”, and τ failure categories are collectively classified as “failed retrievals” for this study.

2.5 Approach for radiative transfer simulation and result comparisons

To address the three SQs for our study (identified in Sect. 1), we performed a total of 14 experiments for each cloud field. The first four experiments were performed with the SHDOM model to study the 3D radiative effects on the observed reflectance and address SQ 1. It involves simulating and comparing R_3D with R_1D for the high- and low-sun cases. The next four experiments involve comparing cloud properties retrieved from R_3D (Box D in Fig. 1) and cloud properties retrieved from R_1D (Box B to A in Fig. 1) for both high and low sun to examine the influence of the 3D radiative effects on the retrieved cloud properties and address SQ 2. These experiments were conducted using the 3D- and 1D-RT-based reflectance as inputs to the precomputed LUT described in Sect. 2.3. The last six experiments were conducted with the I3RC+CKD to examine the impact of the 3D radiative effects on the broadband solar radiative flux for both high- and low-sun scenarios in the LES domains and address SQ 3. These experiments involve calculating $F_{\mathrm{1}\mathrm{D}}^{\ast }$ for each SZA from the retrieved cloud properties using 1D RT as well as computing F_3D and F_1D from the true cloud fields using 3D and 1D RT, respectively. It is important to note that in the $F_{\mathrm{1}\mathrm{D}}^{\ast }$ calculations, the retrieved cloud properties (τ^∗(R_3D) and $r_{\mathrm{e}}^{\ast }\left(R_{\mathrm{3}\mathrm{D}}\right)$) are utilized to calculate the retrieved LWP (using retrieved LWP $\cong \mathrm{2}{\mathit{\tau }}^{\ast }\mathit{\rho }{r}_{\mathrm{e}}^{\ast }/\mathrm{3}$, where ρ is the density of liquid water; Liou, 1992), which is then reconstructed into cloud effective radius and LWC distribution for each LES column while preserving the vertical structure of the original LES cloud field. 1D RT is then performed using the reconstructed retrieved clouds as inputs to obtain $F_{\mathrm{1}\mathrm{D}}^{\ast }$. Note that unless otherwise stated, for this study, the successful and failed retrievals (as described in Sect. 2.4) have been used to represent the total population of cloudy pixels in the cloud property inputs used to calculate $F_{\mathrm{1}\mathrm{D}}^{\ast }$. The calculation of F_1D is identical to that of F_3D except for the absence of the horizontal movement of photons between the LES grid columns. This enables us to determine the impact of neglecting the horizontal movement of photons on the broadband radiative fluxes. On the other hand, in reference to the F_3D, computing $F_{\mathrm{1}\mathrm{D}}^{\ast }$ will not only help us to better understand the implications of neglecting the horizontal transport of photons but will also enable us to measure how biases in the retrieved cloud properties (which are affected by the 3D radiative effects) impact the broadband radiative fluxes.

In order to describe the impact of the 3D radiative effects on the radiance fields, retrieved cloud properties, and broadband radiative flux, we first examine their effects across the LES domain and subsequently quantify their overall impact on the domain by computing the horizontally domain-averaged results to determine the absolute bias, hereafter referred to as “bias” for brevity, defined as $\stackrel{\mathrm{‾}}{y}-\stackrel{\mathrm{‾}}{x}$, where $\stackrel{\mathrm{‾}}{y}$ denotes the domain-averaged result from the 3D RT quantity (e.g., reflectance or flux), and $\stackrel{\mathrm{‾}}{x}$ denotes the domain-averaged result from the 1D RT quantity (e.g., reflectance or flux).

To quantify the difference between the CREs computed from the benchmark F_3D and the CREs computed from F_1D or $F_{\mathrm{1}\mathrm{D}}^{\ast }$, we define a domain-scale quantity known as the relative cloud radiative effects (rCRE) bias as

$$\begin{array}{}\text{(1)}& \text{rCRE bias}=\left(\mathrm{1}-{\displaystyle \frac{{\text{CRE}}_{\mathrm{1}\mathrm{D}}}{{\text{CRE}}_{\mathrm{3}\mathrm{D}}}}\right)\times \mathrm{100},\end{array}$$

where CRE_1D is the CRE calculated from either F_1D or $F_{\mathrm{1}\mathrm{D}}^{\ast }$ in units of W m⁻² and CRE_3D is the CRE calculated from F_3D in units of W m⁻². According to this definition, a rCRE bias of 0 % would indicate that there is no bias between the CRE computed from F_1D or $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and the CRE computed from F_3D. This implies that the CRE computed from F_1D or $F_{\mathrm{1}\mathrm{D}}^{\ast }$ is equivalent to the CRE computed from F_3D. A positive rCRE bias greater than 0 % would quantify the percentage by which the CRE computed from F_1D or $F_{\mathrm{1}\mathrm{D}}^{\ast }$ is lesser than the CRE computed from F_3D and thus indicate that the 1D calculations (F_1D, $F_{\mathrm{1}\mathrm{D}}^{\ast }$) underestimate the CRE relative to F_3D. Also, a negative rCRE bias less than 0 % would quantify the percentage by which the CRE computed from F_1D or $F_{\mathrm{1}\mathrm{D}}^{\ast }$ exceeds the CRE computed from F_3D and imply that the calculations from F_1D or $F_{\mathrm{1}\mathrm{D}}^{\ast }$ overestimate the CRE relative to F_3D.

3 Results and discussion

3.1 Investigating the 3D radiative effects on simulated reflectance

Focusing first on SQ 1, we compare R_1D and R_3D to assess the impact of the 3D radiative effects on the reflectance radiation field, i.e., Box B vs. Box C in the framework of Fig. 1. Specifically, we will investigate the reflectance bias, $\mathit{\delta }R\left(\mathit{\delta }R=R_{\mathrm{3}\mathrm{D}}-R_{\mathrm{1}\mathrm{D}}\right)$, at the two λ values (0.66 and 2.13 µm) required for our bispectral retrieval for both low-sun (SZA 60°) and high-sun (SZA 5°) cases. To describe the 3D radiative effects on the observed reflectance, classifications are made based on the increase in the brightness of a pixel in the LES domain. A pixel in the LES domain is considered “brightened” (“darkened”) if its 3D-RT-based reflectance is higher (lower) than its 1D counterpart.

Maps of δR at λ=0.66 µm (δR_{λ=0.66 µm}) for the two cloud fields when the sun is high and low are shown in Fig. 4. In the low-sun case, the deviation of the 1D-RT-based simulated reflectance from the 3D-RT-based simulated reflectance leads to δR with a distinct pattern of brightening and darkening observed in some pixels across the LES domain. A closer examination of δR_{λ=0.66 µm} within cloudy regions in the low-sun case for the two cloud fields (Fig. 4b, c) reveals a consistent pattern; the brightened pixels, where δR_{λ=0.66 µm} is positive, are predominantly observed in sunlit regions that directly face the sun located on the left (e.g., at X=3.5 km, Y=14 km in Fig. 4b). On the other hand, darkened pixels, where δR_{λ=0.66 µm} is negative, are observed on the opposite side of the cloud layer (e.g., at X=5 km, Y=14 km in Fig. 4b). These findings are consistent with previous 3D radiative effect studies for oblique solar geometry (e.g., Várnai and Davies, 1999; Várnai, 2000; Marshak et al., 2006). The observed opposing effects of brightening and darkening in the low-sun-angle case depend not only on the orientation of the cloud towards or away from the sun; other factors like cloud–cloud interactions, cloud geometry and aspect ratio, spatial distribution of the cloud in the domain, and the horizontal transport of photons also contribute to these behaviors (Várnai and Marshak, 2001, 2002; Marshak and Davis, 2005; Marshak et al., 2006; Zhang et al., 2012).

Figure 4 Maps of the reflectance bias ($\mathit{\delta }R=R_{\mathrm{3}\mathrm{D}}-R_{\mathrm{1}\mathrm{D}}$) for wavelength 0.66 µm at solar zenith angle (SZA) 5° (a, c) for the 27 June and 18 August cases, respectively, and SZA 60° (b, d) for the 27 June and 18 August cases, respectively. The direction of view is at nadir. For SZA 5°, the sun is almost perpendicular to the domain but slightly tilted to the left. For SZA 60° the sun is on the left of the domain.

In the case of the high sun, the sun is almost perpendicular (at SZA 5°), and its radiation interaction with clouds under 3D RT is different from that of the low-sun case. In 3D RT at high sun, the original direction of photons is downwards (due to the sun's small angle of inclination to the vertical), and on striking a cloud, some photons are scattered and some leak from optically thick to optically thin cloudy regions and even out of cloud sides (O'Hirok and Gauiter, 1998) down to the surface where they are absorbed. This is because for photons with a low number of scattering trajectories and high sun, photons leaking out of cloud sides are statistically more likely to continue moving downwards towards the surface where they are absorbed. This leaking of photons to surrounding clouds and the surface results in net photon loss in the thick cloud regions, which explains the darkening of the thick clouds and brightening of the surrounding thin clouds compared to 1D RT results. Hence, δR_{λ=0.66 µm} is mainly negative across the LES domain for the high sun (Fig. 4a, c). The darkening characteristics are more pronounced in the 18 August case because it consists of a larger distribution of thicker clouds compared to the 27 June cloud field; a large number of photons leaking from optically thicker clouds results in a more significant reduction in the reflectance values and more prominent darkening effect than photons leaking from optically thinner clouds. Similar reflectance characteristics are observed for the 2.13 µm band (not shown).

To examine the statistical characteristics of δR in the LES domain, the probability density function (PDF) of δR for “cloudy only” pixels is analyzed to investigate the 3D radiative effects on the observed cloud reflectance. Subsequently, we compared this PDF to δR for both “cloudy and clear-sky” pixels (i.e., the whole LES domain) to highlight the effects of cloud presence on the overall reflectance bias within the LES domain.

The PDFs of δR for cloudy only pixels in the low-sun case (broken black and gray lines in Fig. 5a, b) are characterized by a positive and negative distribution in both the VNIR and SWIR bands (corroborating the brightening and darkening effects in Fig. 4b, d). The overall positive δR observed in the VNIR and SWIR bands (domain mean δR of 0.0351 (0.0292) for the VNIR (SWIR) band in the 27 June case and 0.0379 for the VNIR band in the 18 August case) indicates that the brightening effects are predominant when only cloudy pixels are considered. Meanwhile, δR is −0.0233 for the SWIR band in the 18 August case. This negative δR is due to a high net loss of photons in 3D RT reflectance (more photons leak from clouds to the surface where they are absorbed than those reflected from clouds) compared to the 1D RT results. On the other hand, the PDFs of δR for the cloudy and clear-sky pixels (broken black and gray lines in Fig. 5c, d) are almost similar to that of the cloudy only but show a shift of the distribution leftwards, almost centered around zero. This is expected because clear-sky regions not in the vicinity of any clouds exhibit negligible 3D radiative effects, which causes the distribution to shift closer to zero, since the cloud fraction for both cloud cases is less than 50 %. The horizontal movement of photons from cloudy to surrounding clear-sky regions increase the 3D reflectance of clear-sky areas around the sunlit cloudy regions but the strong darkening effects on the clear-sky regions located opposite the sunlit direction dominate the clear-sky only areas and result in a negative mean bias when the reflectance of clear-sky only pixels is examined. Interestingly, the mean δR for the cloudy and clear-sky pixels is of the same sign as the cloudy only values, which indicates that the cloudy pixels have a significant effect on the domain-scale statistics.

Figure 5 PDF (probability density function) of reflectance bias (δR) for cloudy only pixels for the 27 June case (a) and 18 August case (b). PDF of reflectance bias for cloudy and clear-sky pixels for the 27 June case (c) and 18 August case (d). μ is the domain mean reflectance bias. A Gaussian distribution (solid black curve) with a standard deviation for the 0.66 µm band at SZA 5° and centered around zero is shown in all panels.

The PDFs of δR in the case of the high sun for cloudy only pixels show a larger distribution of pixels with positive δR in both the VNIR and SWIR band accompanied by longer tails to the left (solid red and blue lines in Fig. 5a, b). However, the δR values for both cloud cases are negative in the VNIR and SWIR bands. These observations suggest that large amounts of radiation and/or a large number of photons leak from a small number of thick cloud pixels to a larger number of thin clouds. This phenomenon therefore increases the number of thin clouds with positive reflectance bias, although this is of very small magnitude when compared to the negative biases.

Similar to the low-sun case, the PDF of δR when both cloudy and clear-sky pixels for the high-sun case are considered (solid red and blue lines in Fig. 5c, d) shows a significant distribution of values close to zero. Due to the leaking of photons from thick clouds to thin clouds and clear-sky regions surrounding the clouds, there is an increase in the 3D reflectance of clear-sky regions. Additionally, when the sun is high at SZA of 5°, there are very minimal shadows cast on the clear-sky regions. These two highlighted reasons result in a positive δR for the clear-sky-only region. Thus, the negative value of δR for cloudy and clear sky (same sign as the cloudy only) indicates that the domain-scale reflectance bias is dominated mainly by the cloudy only pixels and they play a significant role in the domain-scale statistics.

3.2 Investigating the 3D radiative effects on cloud retrievals

Focusing on SQ 2 in this section, we investigate how δR, as discussed in the previous section, affects r_e and τ retrievals (i.e., Box A vs. Box D in the framework of Fig. 1). We utilize R_1D as inputs for the LUT (explained in Sect. 2.3) to retrieve the 1D-RT-based cloud droplet effective radius $\left(r_{\mathrm{e}}^{\ast }\right(R_{\mathrm{1}\mathrm{D}}\left)\right)$ and cloud optical thickness (τ^∗(R_1D)) Additionally, we use R_3D as inputs for the LUT to retrieve the 3D-RT-based cloud droplet effective radius $\left(r_{\mathrm{e}}^{\ast }\right(R_{\mathrm{3}\mathrm{D}}\left)\right)$ and cloud optical thickness (τ^∗(R_3D)).

Before discussing analysis of the 3D- and 1D-RT-based retrieval comparison, we first check the accuracy of our retrievals by comparing the original LES cloud properties with our 1D-RT-based retrievals (i.e., comparing retrievals from 1D radiance in Box C with cloud properties in Box A in Fig. 1). For this purpose, the τ from the original LES (τ^true) is the vertical integration of the visible (0.66 µm) extinction coefficient of each column from cloud base to cloud top. For the LES r_e, we the follow Zhang et al. (2017) analytical vertical weighting function (see their Eq. 4) to get the vertically weighted cloud droplet effective radius ($r_{\mathrm{e}}^{\mathrm{VW}}$) where μ_o=0.5, μ=1, and the vertically weighting function parameter (b) associated with the 2.13 µm band was set to 2 to allow for a deeper penetration depth and for better correlation between the $r_{\mathrm{e}}^{\mathrm{VW}}(\mathrm{2.13}\mathrm{\mu }\mathrm{m}$) and bispectral retrievals.

Figure 6 shows the comparison between the $r_{\mathrm{e}}^{\mathrm{VW}}$ (2.13 µm) and the $r_{\mathrm{e}}^{\ast }(R_{\mathrm{1}\mathrm{D}}$) as well as τ^true with the τ^∗(R_1D) for the two cloud fields at SZA = 60° and VZA = 0°. For this comparison, the mean τ and r_e biases are ${\mathit{\mu }}_{\mathit{\tau }\mathrm{bias}}=\langle {\mathit{\tau }}^{\ast }\left(R_{\mathrm{1}\mathrm{D}}\right){\mathit{\tau }}^{\mathrm{true}}\rangle$ and ${\mathit{\mu }}_{r_{\mathrm{e}}\mathrm{bias}}=\langle r_{\mathrm{e}}^{\ast }\left(R_{\mathrm{1}\mathrm{D}}\right)-r_{\mathrm{e}}^{\mathrm{VW}}\left(\mathrm{2.13}\mathrm{\mu }\mathrm{m}\right)\rangle$.

Figure 6 Joint histogram of bispectral retrieved τ based on 1D-RT-simulated reflectance τ^∗(R_1D) vs. vertically integrated τ(τ^true) for the 27 June case and (a) 18 August case (b). Joint histogram of bispectral retrieved r_e based on 1D-RT-simulated reflectance ($r_{\mathrm{e}}^{\ast }\left(R_{\mathrm{1}\mathrm{D}}\right))$ vs. vertically weighted effective radius ($r_{\mathrm{e}}^{\mathrm{VW}}\left(\mathrm{2.13}\mathrm{\mu }\mathrm{m}\right))$ in (c) and (d). The μ_bias is calculated as 〈retrieved cloud properties−reference cloud property〉.

For the two cloud fields considered in this study, the τ^∗(R_1D) is highly correlated with the τ^true as seen in the joint histogram plots (Fig. 6a and b) with a correlation coefficient (r) of 0.9997 for the 27 June case and r of 0.9993 for the 18 August case, although both have a slight positive mean bias (μ_τ bias= 0.1107 and 0.3011 for the 27 June and 18 August cases, respectively). Also, the comparisons of the $r_{\mathrm{e}}^{\ast }(R_{\mathrm{1}\mathrm{D}}$) with the $r_{\mathrm{e}}^{\mathrm{VW}}(\mathrm{2.13}\mathrm{\mu }\mathrm{m}$) in Fig. 6c and d show good correlation (r>0.96) for both cloud cases and slightly positive mean biases $({\mathit{\mu }}_{r_{\mathrm{e}}\mathrm{bias}}=\mathrm{0.0456})$ for the 27 June case as well as a negative mean bias $({\mathit{\mu }}_{r_{\mathrm{e}}\mathrm{bias}}=-\mathrm{0.1343})$ for the 18 August case. Certain extreme outlier bias is observed in the r_e comparisons; these outliers are attributed to thin clouds and have been studied by Miller et al. (2018). Several studies (e.g., Miller et al., 2016, 2018; Zhang et al., 2012) have investigated the accuracy of 1D bispectral retrievals compared to vertically weighted retrievals as well as the impact of cloud vertical profile on bispectral retrievals. Since we have good agreement between retrievals from the 1D-RT-based reflectance and the original LES cloud field properties, this study will use the $r_{\mathrm{e}}^{\ast }(R_{\mathrm{1}\mathrm{D}}$) and τ^∗(R_1D) as the reference cloud properties and directly compare them with the $r_{\mathrm{e}}^{\ast }(R_{\mathrm{3}\mathrm{D}}$) and τ^∗(R_3D) to investigate the impacts of 3D radiative effects on the retrievals.

In the high-sun case retrievals, $r_{\mathrm{e}}^{\ast }(R_{\mathrm{3}\mathrm{D}}$) is overestimated and τ^∗(R_3D) is underestimated compared to their 1D counterpart. This is because photons leaking from optically thick regions to optically thin cloudy regions and out of cloud sides down to the surface where they are absorbed results in a net photon loss, which makes the 3D radiance field appear darker than its 1D counterpart (explained in Sect. 3.1). Consequently, for retrievals, darkening shifts the reflectance observation on the LUT space leftwards and downwards to regions where the LUT r_e grid isolines represent larger droplet sizes and the LUT τ isolines represents thinner clouds. For the low-sun case, $r_{\mathrm{e}}^{\ast }(R_{\mathrm{3}\mathrm{D}}$) is underestimated and τ^∗(R_3D) overestimated in brightened optically thick cloudy pixels (facing the sun), and $r_{\mathrm{e}}^{\ast }(R_{\mathrm{3}\mathrm{D}}$) is overestimated and τ^∗(R_3D) underestimated in darkened pixels on its opposite cloud side. Larger $r_{\mathrm{e}}^{\ast }(R_{\mathrm{3}\mathrm{D}}$) and smaller τ^∗(R_3D) compared to τ^∗(R_1D) and τ^∗(R_1D) in brightened pixels occur since brightening phenomena in the LUT space shift the observed reflectance upwards and rightwards where the LUT r_e grid isolines represent smaller droplet sizes and the LUT τ isolines represent thicker clouds. τ and r_e retrieval biases in satellite observations have been well documented in numerous studies (e.g., Várnai and Marshak, 2002; Zhang and Platnick, 2011; Zhang et al., 2012), and commonly, overestimation of τ^∗(R_3D) is coupled with the underestimation of $r_{\mathrm{e}}^{\ast }(R_{\mathrm{3}\mathrm{D}}$) and vice versa.

Table 2 shows the frequency of failed and successful retrievals from R_3D for the two cloud fields considered in this study. It is observed that the number of failed retrievals is small for the SZA 5° case (< 13 %), while the retrieval failures are larger for the SZA 60° case (> 40 %) for both cloud fields under consideration. The larger retrieval failures for the low-sun case are mostly attributed to multiple scattering in the 3D RT due to increased path length (since the original direction of travel of the photons from the sun is oblique), which increases radiation–cloud interaction and reflectance. Although this leads mostly to τ failures, the other r_e type failures can arise from very darkened pixels (from photon leaking or cloud shadow), which shifts observation outside the LUT lower range (for r_e too large), or brightened pixels from less absorbing clouds, which shifts the observations beyond the upper range of the LUT (for r_e too small) depending on the scenario.

Table 2 Statistics of successful and failed retrievals from the 3D-RT-based radiance for the 27 June and 18 August cloud fields at solar zenith angle (SZA) 5 and 60°. The columns from left to right are the case name (identified by date and time), solar zenith angle (SZA), number of pixels with successful retrievals only, pixels with failed retrievals, and total number of successful and failed retrievals.

Values in parentheses are percentage counts (percentage of counts equals the number of affected pixels divided by the total number of pixels).

3.3 Investigating the 3D radiative effects on the broadband radiative flux

Investigating the 3D radiative effects on the broadband radiative flux: using a combination of the successful and failed retrievals as the input cloud property

Focusing on SQ 3 in this section, we will compare F_3D and $F_{\mathrm{1}\mathrm{D}}^{\ast }$ to investigate the impact of cloud retrieval biases due to the 3D radiative effects on the broadband SW radiative flux. We will also compare F_3D and $F_{\mathrm{1}\mathrm{D}}^{\ast }$ to study the impact of neglecting horizontal photon transport on the broadband SW flux results. Additionally, we compare δF₁ (i.e., $F_{\mathrm{3}\mathrm{D}}-F_{\mathrm{1}\mathrm{D}}^{\ast }$) with δF₂ (i.e., F_3D−F_1D) to determine errors in radiative flux estimates and evaluate the CRE.

It is important to note here that both the successful and the failed retrievals as described in Sect. 2.4 are included in the RT simulations in the control simulations presented in this section. The motivation for including the failed retrievals is to preserve the impacts of this significant fraction of pixels on the domain-averaged fluxes and CRE simulations, even though the retrieval of τ and r_e based on the bispectral method fails for them. In addition to the controlled simulations, we have also conducted sensitivity studies, where we exclude the failed retrievals in the analysis. The results are shown and discussed in the Appendix A.

Maps of the simulated SW broadband radiative quantities (reflected flux at the TOD (F^↑), transmitted flux at the surface (F^↓), and column absorbed flux (F^abs)) for the 27 June case at the high and low sun angles are presented in Figs. 7 and 8, respectively. These figures reveal several interesting and important points. First, it is interesting to note that the reflected flux in Fig. 7d seems blurry in comparison with 1D results in Fig. 7a and g. The same is also seen comparing Fig. 8d with Fig. 8a and g. This is because in 1D RT, simulation of the upwelling hemispheric flux at a given point at the TOD is determined only by the cloud and surface properties in the column beneath such a point. In contrast, in 3D RT simulation, it depends on the cloud and surface properties of both the corresponding column and a large extent of the surrounding columns as a result of a simple parallax effect. Therefore, the contrast between two adjacent columns in the 1D simulation, for example a cloudy column and an adjacent clear-sky column next to it, is quite large, whereas the contrast for the same two columns in 3D simulation is much smaller because the two have a significant overlap in terms of the areas that have influences on their flux. Because of this fundamental difference between 1D and 3D simulations, a pixel-to-pixel comparison of the upwelling flux is not appropriate. Instead, we compare the domain-averaged statistics.

Figure 7 Simulated shortwave broadband reflected flux at the top of the domain (F^↑), transmitted flux at the surface (F^↓), and column absorbed flux (F^abs) derived from the retrieved cloud properties using 1D RT and $F_{\mathrm{1}\mathrm{D}}^{\ast }$ (a–c), derived from the true cloud properties using 3D RT and F_3D (d–f), derived from the true cloud properties using 1D RT and F_1D (g–i), and the difference between $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and F_1D (j–l) for a solar zenith angle of 5° and a view zenith angle of 0°. The sun is high and slightly on the left-hand side of the domain. The solar irradiance at the top of the domain scales with the cosine of the solar zenith angle.

Before we delve into that, we first aim to unravel how cloud property retrieval errors affect 1D RT flux solutions. For this purpose, we compare the F^↑ component of $F_{\mathrm{1}\mathrm{D}}^{\ast }$ (denoted by $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$) with the F^↑ component of F_1D (denoted by $F_{\mathrm{1}\mathrm{D}}^{\uparrow }$). The same comparison is done between the F^↓ component of $F_{\mathrm{1}\mathrm{D}}^{\ast }$ (denoted by $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$) and its counterpart from the F_1D (denoted by $F_{\mathrm{1}\mathrm{D}}^{\downarrow }$). The $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ has visible signatures of the input cloud property retrievals. For instance, in the high-sun case, smaller reflected flux values (recall the underestimated τ dominates retrievals from high sun radiance) dominate $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ (Fig. 7a) compared to $F_{\mathrm{1}\mathrm{D}}^{\uparrow }$ (Fig. 7g). The underestimation of $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ compared to $F_{\mathrm{1}\mathrm{D}}^{\uparrow }$ is evident in Fig. 7j. This difference is also well captured in the domain-averaged values, which will be discussed later in this section. In the low-sun case, comparison between $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ and corresponding $F_{\mathrm{1}\mathrm{D}}^{\uparrow }$ reveals that in $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$, the overestimated retrieved τ areas characterized by thicker clouds (i.e., retrieved from brightened pixels) provide larger reflected flux values and the underestimated retrieved τ areas characterized by thinner clouds (i.e., retrieved from darkened pixels) have smaller reflected flux values than their $F_{\mathrm{1}\mathrm{D}}^{\uparrow }$ counterpart (Fig. 8j). Their overall effect on the domain reflected flux values depends on how the opposite 3D radiative effects (cloud side brightening and darkening) mitigate each other.

Figure 8 Simulated shortwave broadband reflected flux at the top of the domain (F^↑), transmitted flux at the surface (F^↓), and column absorbed flux (F^abs) derived from the retrieved cloud properties using 1D RT and $F_{\mathrm{1}\mathrm{D}}^{\ast }$ (a–c), derived from the true cloud properties using 3D RT and F_3D (d–f), derived from the true cloud properties using 1D RT and F_1D (g–i), and the difference between $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and F_1D (j–l) for a solar zenith angle of 60° and view zenith angle of 0°. The sun is on the left-hand side of the domain. The solar irradiance at the top of the domain scales with the cosine of the solar zenith angle.

An examination of $F_{\mathrm{1}\mathrm{D}}^{\downarrow }$ and $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$ for the high-sun case reveals that $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$ beneath clouds is larger compared to $F_{\mathrm{1}\mathrm{D}}^{\downarrow }$, while they have same values in clear-sky regions (Fig. 7k) This is expected since $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ is less than $F_{\mathrm{1}\mathrm{D}}^{\uparrow }$. Thus, the amount of $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$ at the surface beneath the cloud increases. For the low-sun case, $F_{\mathrm{1}\mathrm{D}}^{\downarrow }$ and $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$ beneath the clouds have higher values where the TOD reflected flux is low and lower values where the TOD reflected flux is high (Fig. 8k).

The radiative quantities across the LES domain have different characteristics in F_3D stemming from the horizontal transport of photons across pixels. For the high-sun case, aside from the blurriness of the F^↑ component of F_3D (denoted by $F_{\mathrm{3}\mathrm{D}}^{\uparrow }$) which was previously explained, an examination of the F^↓ component of F_3D (denoted by $F_{\mathrm{3}\mathrm{D}}^{\downarrow }$) in Fig. 7e shows a slight tilt of the cloud shadows according to the angle of projection of the sun (located at SZA 5° to the left). It also reveals enhanced $F_{\mathrm{3}\mathrm{D}}^{\downarrow }$ values around cloud edges. Such is not the case in $F_{\mathrm{1}\mathrm{D}}^{\downarrow }$ (Fig. 7h) due to 1D RT setup where each cloudy column is considered independent. These observations are consistent with findings made by Gristey et al. (2020) for similar shallow cumulus cloud fields. Gristey et al. (2020) showed that enhanced $F_{\mathrm{3}\mathrm{D}}^{\downarrow }$ around cloud edges is primarily caused by the diffused component of the transmitted flux, which is scattered by clouds towards clear-sky regions beyond cloud shadows. In $F_{\mathrm{3}\mathrm{D}}^{\downarrow }$ for the low sun (Fig. 8e), due to 3D RT and the more oblique solar angle (SZA 60°), there is an increase in the total effective cloud cover (Di Giuseppe and Tompkins, 2003; Tompkins and Di Giuseppe, 2007), as well as an increase in the size of the cloud shadow, which reduces the transmitted flux at the surface (e.g., around X=7, Y=6 km) in Fig. 8e. Just as in the case of the high sun, these features are absent in the low-sun 1D RT runs ($F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$ and $F_{\mathrm{1}\mathrm{D}}^{\downarrow }$). An analysis of the domain-averaged statistics will help shed more light on the differences between the 3D RT and 1D RT radiative flux results on the domain scale.

Table 3 Statistics of successful and failed retrievals from the 3D-RT-based radiance for the 27 June and 18 August cloud fields at solar zenith angle (SZA) 5 and 60°. The columns from left to right are the case name (identified by date and time), solar zenith angle (SZA), number of pixels with successful retrievals only, pixels with failed retrievals, and total number of successful and failed retrievals.

Note: values before the parentheses are calculated from the combination of failed and successful retrievals representing the total cloudy population, while values in parentheses are calculated from successful retrievals only representing the total cloudy population. Clear-sky pixel values have been included in all calculations.

The domain-averaged broadband F^↑, F^↓, and F^abs components of $F_{\mathrm{1}\mathrm{D}}^{\ast }$, F_1D, and F_3D for the 27 June and 18 August cases at SZA 5 and SZA 60° are reported in Table 3. As previously explained, the predominant photon leaking associated with high-sun 3D RT and the ensuing underestimation of the retrieved τ, which dominate the cloud property retrievals from high-sun 3D-simulated reflectance, increase the number of retrieved optically thinner clouds (relative to the original LES τ used in F_1D calculations) utilized as inputs for the $F_{\mathrm{1}\mathrm{D}}^{\ast }$ calculations. This leads to the underestimation of the domain-averaged $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ compared to $F_{\mathrm{1}\mathrm{D}}^{\uparrow }$; In the 27 June case, the domain-averaged $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ (215.44 W m⁻²) is underestimated compared to the corresponding $F_{\mathrm{1}\mathrm{D}}^{\uparrow }$ value (225.37 W m⁻²) by about 9.93 W m⁻², while in the 18 August case, the domain-averaged $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ (315.16 W m⁻²) is underestimated compared to the corresponding $F_{\mathrm{1}\mathrm{D}}^{\uparrow }$ value (355.26 W m⁻²) by 40.1 W m⁻². The larger value of the underestimated domain-averaged $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ in the 18 August case stems from its larger cloud fraction and τ bias. The transmitted flux at the surface below clouds is dependent on the amount of flux reflected towards the TOD; lower reflected flux values indicate that less radiation is reflected from the clouds, which allows a greater amount of radiative flux to be transmitted to the surface beneath the clouds. This reason, coupled with the overestimation of the transmitted flux at the surface due to missed thin clouds in our bispectral retrievals (red regions in Fig. 3, retrieved τ=0 for VNIR reflectance less than the smallest LUT τ value), explains why for the high-sun case the domain-averaged $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$ values are higher compared to $F_{\mathrm{1}\mathrm{D}}^{\downarrow }$ values, resulting in differences of 8.21 and 35.13 W m⁻² for the 27 June and 18 August cases, respectively, although for the high sun angle the contribution of the missed thin clouds to the overestimation of $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$ beneath clouds in our case study is small (constituting about 0.23 % and 0.34 % of the domain-averaged surface transmitted flux for the 27 June and 18 August high cases, respectively).

Comparing results from the three sets of experiments in Table 3 reveals that for the high-sun case, the $F_{\mathrm{1}\mathrm{D}}^{\ast }$ results clearly agree better with the benchmark F_3D than the F_1D results. In the 27 June case, δF₁ values for the domain-averaged F^↑, F^↓, and F^abs are 0.49, −0.18, and −0.33 W m⁻², respectively, which are significantly smaller in magnitude than those for δF₂ (−9.44, 8.03, and 1.41 W m⁻², respectively). Similarly, for the 18 August case, δF₁ for the domain–averaged F^↑, F^↓, and F^abs are −6.48, 6.81, and −0.41 W m⁻², respectively, compared to corresponding biases of −46.58, 42.04, and 4.59 W m⁻² for δF₂. These results suggest that $F_{\mathrm{1}\mathrm{D}}^{\ast }$ gives an overall better radiative energy estimate than F_1D for the high-SZA case. In the low-sun case, the $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and F_1D are very close to each other and there is no clear winner when compared to the benchmark 3D RT results. In the 27 June case, the $F_{\mathrm{1}\mathrm{D}}^{\ast }$ agrees slightly better with 3D results than the F_1D, but the opposite is true in the 18 August case. This result seems to suggest that although in the low-sun case the brightening and darkening effects can lead to large retrieval biases, they tend to cancel each other out in the flux computations. Interestingly, both 1D results tend to underestimate F^↑ and overestimate F^↓. This is probably because the brightening effect is dominant in the 3D RT, leading to some extremely bright pixels. But they are not captured in the 1D RT computations, even in the $F_{\mathrm{1}\mathrm{D}}^{\ast }$ using the upper limit of τ=158.48 in the flux computation. Thus, the reflected flux quickly reaches the asymptotic value when τ is large, and therefore simply using larger τ value in 1D RT cannot simulate the extreme brightness of clouds due the brightening effect in 3D RT. Results for δF₂ computed from the transmitted flux at the surface for both cloud cases (27 June and 18 August) are positive when the sun is high (8.03 and 42.04 W m⁻²) and negative for the low sun angle (−6.61 and −15.39 W m⁻²), consistent with Gristey el al. (2020) for surface irradiance showing positive domain mean δF₂ in the afternoon (high sun) and negative domain mean δF₂ towards the end of the day (low sun).

Figure 9 Relative cloud radiative effect (rCRE) bias computed from the successful + failed retrievals for the (a) reflected flux at the top of the domain, (b) surface transmitted, and (c) absorbed for the two cloud fields.

Because both cloud cases have a cloud fraction lower than 50 %, the domain-averaged statistics include a large fraction of clear-sky pixels. Now we focus our scope on cloudy pixels and investigate the differences in CRE. The rCRE bias provides a quantitative estimate of how these biases affect the CRE. For the two cloud cases considered in this study, plots of rCRE bias computed from the F^↑, F^↓, and F^abs at SZA 5 and 60° for $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and F_1D relative to F_3D are presented in Fig. 9. In the 27 June case, the rCRE bias of 0.97 % computed from the high-sun $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ result indicates a negligible deviation (less than 1 %) from the benchmark CRE, while the rCRE bias of −19 % computed from the high-sun $F_{\mathrm{1}\mathrm{D}}^{\uparrow }$ result shows that the bias is quite substantial. Similarly, for the 18 August case, the rCRE bias computed from the high-sun $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ result is less than 5 %. On the other hand, the rCRE bias of −32.48 % computed from the high-sun $F_{\mathrm{1}\mathrm{D}}^{\uparrow }$ result shows that the bias is quite large. Similar results are obtained for rCRE bias computed from F^↓. In the 27 June case, the rCRE bias computed from the $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$ is 0.33 % (Fig. 9b, second bar on the left), which shows minimal bias less than 1 %, while the rCRE computed from the $F_{\mathrm{1}\mathrm{D}}^{\downarrow }$ is −14.5 % (Fig. 9b, first bar on the left). Similarly, for the 18 August case, the rCRE bias computed from the $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$ is −4.12 % (Fig. 9b, second bar on the right), while the rCRE bias computed from the $F_{\mathrm{1}\mathrm{D}}^{\downarrow }$ is −25.10 % (Fig. 9b, first bar on the right).

When the absorbed flux is taken into consideration, for the 27 June high-sun case, the rCRE bias computed from the absorbed component of $F_{\mathrm{1}\mathrm{D}}^{\ast }$ ($F_{\mathrm{1}\mathrm{D}}^{\ast \mathrm{abs}}$) is −6.05 % (Fig. 9c, second bar on the left), which is less bias compared to the 25.64 % rCRE bias computed from the absorbed component of F_1D ($F_{\mathrm{1}\mathrm{D}}^{\mathrm{abs}}$) (Fig. 9c, first bar on the left). Similarly, for the 18 August case, the rCRE bias computed from $F_{\mathrm{1}\mathrm{D}}^{\ast \mathrm{abs}}$ is −1.73 % (Fig. 9c, second bar on the right), while the rCRE bias computed from $F_{\mathrm{1}\mathrm{D}}^{\mathrm{abs}}$ is 19.09 % (Fig. 9c, first bar on the right). For the low-sun case, the rCRE biases computed from $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and F_1D are comparable, which is consistent with the domain-averaged statistics in Table 3. Evidently, both $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and F_1D overestimate the CRE at TOD and surface, which means an underestimation of cloud reflection and overestimation of transmission. This is consistent with the results in Table 3.

Overall, the above analysis indicates that the $F_{\mathrm{1}\mathrm{D}}^{\ast }$ provides results that are better than (in the high-sun case) or at least comparable to (in the low-sun case) F_1D for both domain-averaged flux statistics and CRE when compared to the benchmark F_3D results. With these results we can conclude that the CRE calculated with 1D RT using retrieved cloud properties, which are biased due to the 3D effects, is found to be comparable to or better than the CRE calculated with 1D RT using the true cloud properties.

4 Summary and conclusion

It is well known that bispectral cloud property retrievals based on the 1D RT have significant errors due to the 3D radiative effects. In this study, we investigate whether the biased retrievals can still be used to estimate the broadband flux and CRE. To address this question, we selected two cloud fields from the LASSO activity, one on 27 June 2015 and another on 18 August 2016, to serve as case studies for our research. The LES cloud fields have different microphysics with different CBH and CTH, and the value of the cloud fraction for the 18 August 2016 cloud field (47.08 %) is more than twice that of the 27 June 2015 (20.15 %) cloud field. Radiance simulations, bispectral retrievals, and broadband SW flux radiative transfer simulations were performed using these cloud fields at two SZAs, a high-sun case (SZA = 5°) and a low-sun case (SZA = 60°), and the results were analyzed. The flux computations were carried out in three sets. The reference broadband SW flux calculations were performed using the cloud properties from the original LES cloud field under 3D RT (F_3D). We also computed similar RT broadband SW flux calculations with the same cloud properties from the original LES cloud field except that the RT calculations were computed using 1D RT (F_1D). Additionally, we computed the last set of broadband SW flux calculations using 1D RT and bispectrally retrieved cloud properties as inputs ($F_{\mathrm{1}\mathrm{D}}^{\ast }$).

The high sun radiance results for the two cloud fields show that in the 3D RT high-sun case, the photons leaking from optically thick cloudy regions to optically thin cloudy regions and the surface dominate the LES reflectance field. This results in overestimated r_e and underestimated τ dominating the cloud property retrievals. Results from the low-sun case, for the two cloud fields considered, show that in comparison to the 1D RT radiance fields, brightening and darkening effects both occur in the 3D-RT-simulated radiance observation. Therefore, retrievals from the low-sun 3D radiance observations are characterized mainly by both overestimation of τ and underestimation of r_e in brightened pixels and underestimation of τ and overestimation of r_e in darkened pixels. The cumulative effects of these brightening and darkening or photon leaking effects and their impacts on the retrieved cloud properties dictate their impact on the broadband radiative flux.

The results from the broadband SW radiative flux computation showed that, although the bispectrally retrieved cloud properties are often biased due to the 3D radiative transfer effects, for high-sun cases, calculations of the CRE from these $F_{\mathrm{1}\mathrm{D}}^{\ast }$ values agree well with the benchmark values (which is the F_3D in our case) with agreement within 7 % for rCRE bias calculations from the reflected, transmitted, and absorbed fluxes in the high-sun cases. Conversely, the rCRE bias computed from the F_1D quantities could reach about 33 %. Thus, for high-sun situations, the $F_{\mathrm{1}\mathrm{D}}^{\ast }$ provides consistently better estimates of the CRE than the F_1D. For the low-sun case, the two 1D RT experiments provide comparable results, both underestimating cloud reflection and overestimating transmission, and there is not a clear winner when compared to the 3D RT benchmark.

The influence of the failed retrievals on the CRE was also investigated (see details in Appendix A), with results indicating that for the high-sun case, the impact of the failed retrievals on the radiative flux quantities is negligible, with less than 6 % changes observed in the rCRE bias computed from the domain-averaged TOD reflected, surface transmitted, and absorbed $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and F_1D results. Such is not the case for the low-sun case where the failed retrievals have a very huge impact on the radiative flux quantities. Excluding the failed retrievals from the domain-averaged reflected, transmitted, and absorbed $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and F_1D low-sun case analysis could increase the rCRE bias by as much as a factor of 6 compared to values which included the failed retrievals in the analysis. Whether or not to always use the failed retrievals in the radiative flux and CRE estimation is still an important question, especially how best to filter out the failed retrievals from cloud properties retrieved from instruments that rely on the bispectral method (e.g., in MODIS cloud products) for use in radiative flux estimation. We observed here that filtering out all failed retrievals, especially from the low sun angle, can greatly impact the radiative flux estimates. Thus, efforts should be made to study which category of failed retrievals is most relevant for use in CRE estimation.

In conclusion, despite the potential biases due to the 3D radiative effects, the retrieved cloud properties based on 1D RT from the bispectral method still provide CRE estimates that are comparable to or better than CRE calculated from the true cloud properties using 1D RT. Some future questions that warrant answers involve how the 3D radiative effects affect the broadband fluxes for different cloud arrangements and other types of clouds, such as deep convective clouds. Also, while we have considered only the nadir viewing angle in this work, previous studies (e.g., Várnai and Marshak, 2007) have shown that the biases of 1D cloud retrievals vary systematically with viewing direction; therefore, the impacts of off-nadir viewing directions on the broadband flux need to be investigated. Another important study will be to determine how changes in surface albedo and type affect our results. Additionally, while our case study mainly focused on the impact of the 3D radiative effects on SW fluxes, the impact of the 3D radiative effects on LW radiation is important and needs to be investigated.

Appendix A: Impacts of failed retrievals on the radiative flux

The calculations of $F_{\mathrm{1}\mathrm{D}}^{\ast }$ and domain radiative flux analysis in Sect. 3.3 utilize both the successful and failed retrievals (categorized in Sect. 2.4) to represent the total population of cloudy pixels. Henceforth, both successful and failed retrievals as a representative of the total population of cloudy pixels will be referred to as “all retrieved cloud pixels”. In this Appendix, our focus is to examine and compare the TOD reflected, surface transmitted, and column absorbed radiative fluxes when the failed retrievals are excluded from the radiative flux analysis. This will help to diagnose if using solely successful retrievals as a representative of the total population of cloudy pixels in the LES domain will produce the correct radiative energy estimates and thus provide information on the radiative properties of the excluded failed retrievals.

An examination of the high-sun domain-averaged F^↑, F^↓, and F^abs for both LES cloud cases, when only successful retrievals represent the total population of cloudy pixels in the $F_{\mathrm{1}\mathrm{D}}^{\ast }$ calculations, shows minimal changes (within the range ±1.9 W m⁻²) from previous values which utilized all retrieved cloud pixels in the radiative flux analysis (Table 3). This is due to the small number of failed retrievals in the high-sun scenario (< 14 % for both cloud cases; Table 2). But this is not the case for the low-sun case, where changes between the two aforementioned calculations are large, reaching up to ±35.96 W m⁻² (Table 3). These large changes are because of the large number of failed retrievals from strong 3D radiative effects (> 43 % for both cloud cases; Table 2) as well as different radiative behavior of the failed retrieval categories observed in the low-sun scenario. Figure A1 shows plots of successful and failed retrievals categories (classified as described in Sect. 2.4) from the high and low sun radiance for the 27 June and 18 August cases. From these plots, it is observed that when the SZA is 60°, the “r_e too small” failures are predominant around cloud edges in the sunlit areas. The τ failures are observed mostly in the illuminated sunlit cloudy regions and the “r_e too large” failures occur mostly on the opposite sides where the shadowing effect is dominant (Fig. A1b and d). For the high sun at SZA 5°, τ failures are almost negligible because the VNIR reflectance observations do not exceed the LUT τ upper limit of 158.48, while there is a small number of occurrences of the “r_e too large” and “r_e too small” failures (Fig. A1a and c).

It should be noted that when we exclude the failed retrievals from the broadband flux analysis, we keep the total cloud fraction constant. In other words, we scale the broadband flux based on the successful pixels by the ratio of total cloudy to successful pixels such that the effect of cloud fraction reduction is removed from the analysis. The impacts of excluding failed retrievals on the domain-averaged broadband flux can be assessed by comparing the values outside the parentheses with those inside in Table 3 and better understood in light of the failed retrieval statistics given in Table 2.

Results of $F_{\mathrm{1}\mathrm{D}}^{\ast }$ for the 27 June case at SZA 5° show that the domain-averaged $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ is underestimated by 1.50 W m⁻² (213.94 W m⁻² in comparison to 215.44 W m⁻²) when only successful pixels are used to represent the total population of cloudy pixels compared to results which utilize all retrieved cloud pixels in the radiative flux analysis. This is mainly because the dominant type of retrieval failure in this case is the “r_e too small” failure, accounting for about 71 % of the failed pixel retrieval statistics (see Table 2). Recall that “r_e too small” failure is mainly a result of the brightening effect, and therefore associated pixels appear brighter in 3D RT than 1D RT. As a result, excluding these pixels leads to an underestimate of the domain-averaged broadband reflected flux. For the same reason, excluding these pixels leads to an overestimation of transmitted flux at the domain bottom.

In contrast to the 27 June case, excluding the failed retrievals in the $F_{\mathrm{1}\mathrm{D}}^{\ast }$ for the 18 August case leads to an overestimation of domain-averaged $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ and underestimation of the $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$. This is probably because the dominant failed retrieval type is “r_e too large”, which is because of the darkening effect. These pixels appear darker from the perspective of TOD and more transmissive from the perspective of the bottom in 3D RT than 1D RT. For comparison purposes, we have also excluded the failed pixels from the F_1D calculations. Overall, the results are very similar to and consistent with those based on $F_{\mathrm{1}\mathrm{D}}^{\ast }$.

In comparison with the high-sun case, the impacts of failed retrievals on the broadband flux statistics are much larger in the low-sun SZA 60° case. In both LES cases, the exclusion of failed retrievals leads to a significant decrease in domain-averaged $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ and increase in $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$. For example, in the 27 June case, the $F_{\mathrm{1}\mathrm{D}}^{\ast \uparrow }$ decreases from 134.22 W m⁻² when failed pixels are included to 111.21 W m⁻² when they are excluded, which is accompanied by an increase in $F_{\mathrm{1}\mathrm{D}}^{\ast \downarrow }$ from 419.60 to 441.77 W m⁻². A close look at Table 2 reveals that in both LES cases, the combination of “r_e too small” and τ failures accounts for the majority of failed retrievals: 95 % in the case of 27 June and 68 % in the 18 August case. As mentioned above, both types of failures are because of the brightening effect. Excluding them is expected to cause underestimation of domain-averaged reflected flux and overestimation of the transmitted flux.

The impacts of excluding failed retrievals on the rCRE bias are shown in Fig. A2. A comparison to the results in Fig. 9 reveals two points. First, the biases in the low-sun cases become much larger, which is expected because there are many more failed retrievals in these cases. Second, it is evident that the flux estimates derived from the retrieved clouds using 1D RT still provide a better (in the case of high sun) or comparable (in the case of low sun) approximation to the flux estimates from the true cloud fields using 3D RT simulations in comparison with those derived from the true cloud fields using 1D RT. Therefore, our conclusion made based on the statistics of all retrievals still holds when failed retrievals are excluded from the analysis. On the other hand, it is also evident that to achieve a better comparison with the flux derived from the true clouds using 3D RT, it is better to include the failed retrievals to preserve the effects of 3D RT.

Figure A1 Plots of successful and failed retrieval categories for the 27 June 2015 and 18 August 2016 cases at solar zenith angle 5° (a, c) and solar zenith angle 60° (b, d).

Figure A2 Relative cloud radiative effect bias computed from only the successful retrievals: top of the domain reflected in panel (a), surface transmitted in panel (b), and column absorbed flux in panel (c) for the two cloud fields.

Appendix B: Surface spectral albedo plot

Figure B1 Surface spectral albedo plot utilized in the study. Bi (i ranges from 1 to 14) represents the bands in parentheses.