Impacts of the three-dimensional radiative effects on cloud droplet number concentration retrieval and regression-based albedo susceptibility estimates

Cloud droplet number concentration (N_d) in warm liquid clouds plays a crucial role in understanding cloud microphysical processes and the influence of aerosol–cloud interactions (ACI) on Earth's climate. N_d from satellite-retrieved cloud properties such as the cloud optical thickness (τ) and cloud droplet effective radius (r_e) can be biased due to the three-dimensional (3D) radiative transfer (RT) effects. Using Large-Eddy Simulation (LES) cloud fields and RT simulations, this study investigates how biases in cloud property retrievals caused by 3D-RT effects impact the derived N_d and subsequent regression-based albedo-susceptibility estimates. Our sensitivity studies confirm that the bi-spectral retrievals using the 3.7 µm channel – whose r_e retrieval is closest to cloud top – shows better agreement with N_d from our LES models, compared to results based on the 1.6 and 2.1 µm retrievals. At native LES resolution, N_d across all absorbing channels is strongly impacted by the 3D-effects, with the magnitude depending on the solar zenith angles (SZAs); on average, for high/low sun conditions N_d under 1D-RT overestimates/underestimate its 3D-RT counterpart, which indicates dominant darkening/brightening effects. At coarser satellite-like resolutions, average statistics between 1D and 3D retrievals agree better, indicating compensation between 3D and plane-parallel averaging assumptions. Furthermore, our regression-based albedo-susceptibility results evaluated at the top of cloud domain show that 3D RT greatly modifies the local α–N_d relationship at native LES resolution, especially under more oblique solar geometry, but the differences between 1D- and 3D-regression-based susceptibility estimates are substantially reduced at coarser resolution and low-to-moderate LWP/τ regimes. These results indicate that although 3D RT can strongly affect domain-top regression-based susceptibility at the LES resolution, its impact is reduced at satellite-like spatial resolution for low to moderate LWP/τ regimes.

1 Introduction

Aerosols influence Earth's radiative budget both directly by scattering and absorbing radiation and indirectly by modifying cloud properties. An increase in aerosol particles that serve as cloud condensation nuclei (CCN) in warm clouds typically leads to a decrease in cloud droplet effective radius (r_e), and an increase in the cloud droplet number concentration (N_d). This results in a more reflective cloud for a given liquid water path (LWP), known as the first aerosol indirect effect or Twomey effect (Twomey, 1974, 1977), leading to the radiative forcing from aerosol–cloud interactions (ACI), (RF_aci) (Bellouin et al., 2020; Forster et al., 2021). In addition to the Twomey effect, increased aerosol concentrations can cause additional cloud responses, such as changes in LWP and suppression of precipitation (Albrecht, 1989).

In most remote sensing–based studies of ACI, aerosol effects on clouds are estimated through changes in N_d. Therefore, continuous and accurate measurements of N_d at regional and global scales are essential to improve understanding of aerosol impacts on cloud microphysical and optical properties, as well as to evaluate RF_aci. Satellite-based remote sensing is commonly employed for such observations because alternative methods for measuring N_d are often limited or unavailable.

Operational passive satellite instruments such as the Moderate-resolution Imaging Spectroradiometer (MODIS) do not directly retrieve N_d. Instead, N_d is derived from other retrieved parameters – namely, cloud optical thickness (τ) and r_e – under assumptions about cloud adiabatic growth and the constancy of N_d throughout the cloud depth (Boers et al., 2006; Grosvenor et al., 2018; Quaas et al., 2006). Such estimates of N_d derived from passive imager observations rely on the τ and r_e retrieved via the so-called bi-spectral retrieval techniques (e.g., Nakajima and King, 1990; Twomey and Seton, 1980). The bi-spectral retrieval method simultaneously retrieves τ and r_e using cloud reflectance measurements from two spectral bands. Typically, one band is selected from a non-absorbing visible or near-infrared (VNIR) spectral region (e.g., wavelength centered near 0.86 µm), which is primarily sensitive to τ, while the other is chosen from a moderately absorbing shortwave infrared (SWIR) (e.g., wavelength centered near 2.1 µm) or mid-wave infrared (MWIR) (e.g., wavelength centered near 3.7 µm) spectral region, which is more sensitive to r_e.

While bi-spectral retrievals from passive instruments like MODIS have significantly advanced our understanding of clouds, N_d derived from these retrievals may contain large errors and uncertainties. These arise from assumptions made both in the retrieval process and in the computation of N_d, which affect their applicability for aerosol–cloud interaction and other process studies (Gryspeerdt et al., 2016; McCoy et al., 2017; Quaas et al., 2020). All operational bi-spectral retrievals rely on one-dimensional (1D) radiative transfer (RT) theory, which assumes that the atmosphere within each pixel is horizontally homogeneous (plane-parallel assumption) and that each pixel is independent of its neighbors (independent pixel assumption), primarily for computational efficiency. However, real clouds have complex 3D structures, requiring RT simulations that account for both vertical and horizontal radiation transport, known as “3D RT”. This is challenging because detailed 3D cloud structures are generally unknown a priori, and even when available, 3D RT simulations are computationally expensive and not operationally feasible. The deviation of the 3D RT in real clouds from the 1D RT theory as the backbone of operational cloud retrieval algorithm is often referred to as the 3D radiative effects. These effects can introduce substantial biases in cloud property retrievals and radiative quantities that are based on 1D RT (Marshak et al., 2006; Várnai and Marshak, 2002; Zhang et al., 2012, 2016; Zinner et al., 2010; Cornet and Davies, 2008). This paper focuses on investigating how errors due to 3D radiative effects associated with retrieved cloud properties (τ and r_e) impact derived N_d statistics, and how these errors affect regression-based albedo susceptibility estimates and interpretations of the first aerosol indirect effect. Evaluating 3D effect errors on derived N_d from LES model fields can improve constraints on N_d, thereby provide better understanding of ACI and RF_aci.

From the perspective of results, 3D RT can cause brightening or darkening phenomena in observed cloud reflectance compared to 1D RT. Brightening occurs when cloud reflectance in 3D RT is higher than in 1D RT, while darkening occurs when clouds appear darker under 3D RT relative to 1D RT. These effects can significantly impact the retrievals of τ and r_e (Kato and Marshak, 2009; Marshak et al., 2006; Várnai and Marshak, 2002). For example, brightening effects typically lead to overestimated τ and underestimated r_e retrievals, while the darkening effects typically result in underestimated τ and overestimated r_e when compared to 1D RT retrieved results (Marshak et al., 2006). Several factors can contribute to 3D radiative effects in broken cloud fields, including variability in cloud-top height, cloud horizontal and vertical heterogeneity, solar and viewing geometry, and light scattering from optically thick to optically thin regions. Previous studies have shown that solar and observation geometry strongly modulate how these 3D effects appear in bi-spectral cloud-property retrievals, influencing both the sign and magnitude of retrieval errors in τ and r_e (Kato and Marshak, 2009; Marshak et al., 2006; Zhang et al., 2012), with these errors being generally more pronounced in broken cloud fields.

Over the past decade, several studies related to aerosol–cloud interactions have utilized N_d estimates or the N_d–LWP relationship obtained from bi-spectral retrievals to study the impact of aerosols on clouds, evaluate RF_aci, and quantify related uncertainties. For example, Grosvenor et al. (2018) estimated that N_d from passive remote sensing instruments at the pixel scale have a relative uncertainty of 78 % (dominated by errors in r_e retrievals) for optically thick single layer stratiform clouds. Quaas et al. (2020) reviewed the challenges and uncertainties in constraining the Twomey effect from satellite observations, suggesting that past studies have likely underestimated the sensitivity of N_d to aerosol perturbations, leading to an underestimation of the Twomey effect's radiative forcing. Arola et al. (2022) used satellite and simulated data to show that natural spatial variability and errors in bi-spectral retrievals of τ and r_e propagate into N_d and LWP estimates and can cause positive LWP adjustments to be misinterpreted as negative, often leading to an underestimation of the cooling effect of ACI. Recent studies by Loveridge and Di Girolamo (2024) used a combination of synthetic cloud fields, RT simulations, and satellite observations to show that commonly used sampling strategies for N_d derived from bi-spectral retrievals, do not eliminate systematic biases caused by cloud heterogeneity, potentially leading to an overestimation of aerosol indirect effects in climate studies.

While previous studies have quantified uncertainties in N_d derivations, investigated cloud adjustments under varying environmental conditions, constrained the Twomey effect using satellite observations, and tested filtering strategies on retrieved N_d, few have specifically examined how errors arising from 3D radiative effects propagate into N_d calculations and influence the estimation of the Twomey effect from satellite data (e.g., Arola et al., 2022). Understanding this linkage is crucial, as 3D radiative effects can introduce biases in cloud property retrievals that may significantly affect albedo susceptibility estimates and our interpretation of ACI.

The goal of this work is to build on previous studies by investigating how 3D radiative effects, which influence retrievals of τ and r_e, propagate into N_d derivations calculated from the retrievals. We further aim to examine how these effects impact regression-based albedo susceptibility estimates derived from such retrievals. Specifically, using our model LES fields, we seek to answer the following questions: How do the brightening and darkening effects that impact the retrieved r_e and τ affect calculations of N_d? Does the derived N_d change significantly when different wavelength pairs are utilized in the bi-spectral retrievals used for the N_d calculations? And how do the 3D effects impact these changes? Finally, how do 3D effects impact albedo susceptibility estimates derived from bi-spectral retrievals?

The paper's remaining structure is arranged as follows: Sect. 2 briefly describes the data and theory for the study. Results and discussions on how the 3D radiative effects influence calculations of N_d and regression-based albedo susceptibility estimates are presented in Sect. 3. The summary and conclusion are given in Sect. 4.

2 Data and Theory

2.1 Cloud Field Dataset

Since clouds in the real atmosphere are always influenced to some extent by 3D radiative effects, relying solely on observational data to understand these effects and their impacts on derived N_d is challenging. To address this challenge, many studies (e.g., Ademakinwa et al., 2024; Miller et al., 2018; Rajapakshe and Zhang, 2020; Zhang et al., 2012) have utilized synthetic cloud fields and RT simulations to mimic the satellite observation-retrieval process and study the 3D radiative effects on retrieved τ and r_e. A major advantage of this approach is that the LES cloud field provides the model “truth” which is difficult to obtain in real world observations. Also, the flexibility of this approach allows investigation of various cloud mechanisms at different levels of complexity. For example, we can analyze the 3D RT impacts by comparing retrievals from 3D RT with those from 1D RT. Disadvantages of this approach are that the usefulness of such a study depends on the realism of the synthetic microphysical and water content fields; which, like the satellite retrievals themselves are difficult to validate, as well as the limited variability in scene properties that can be run and processed compared to global satellite sampling.

This study utilizes a similar state-of-the-art satellite retrieval simulator as in Zhang et al. (2012) and Ademakinwa et al. (2024) that applies 1D and 3D RT simulations to synthetic cloud fields. In this work, our cloud fields are Large Eddy Simulation (LES) models (Distributed Hydrodynamic-Aerosol-Radiation Modeling Application (DHARMA)) with bin microphysics (Ackerman et al., 2004; Miller et al., 2016; Zhang et al., 2012). They consist of different initial aerosol loadings derived from an idealized case study by Stevens et al. (2001) conducted during the Atlantic Trade Wind Experiment (ATEX). Three LES cases of different initial aerosol loadings are examined throughout this study. The first case has a CCN loading of 40 cm⁻³ (hereafter referred to as “ATEX clean”), the second case has a CCN loading of 75 cm⁻³ (hereafter referred to as “ATEX control”), while the third case has a CCN loading of 600 cm⁻³ (hereafter referred to as “ATEX polluted”). The LES provides a model “truth” for the 3D cloud microphysical properties, which are used as a baseline for comparisons with numerically simulated retrievals. The droplet sizes in the LES consist of 25 bins (0.874 µm ≤ radius ≤ 281.76 µm) used to represent the distribution of droplet sizes (Ackerman et al., 1995) and Mie scattering properties are bulk averaged over a highly resolved flat sub-bin to obtain the optical properties for each size bin. This serves as input for RT simulations based on the size distributions of the LES cloud fields. The LES cases have a domain size of $\mathrm{9.6}\times \mathrm{9.6}\times \mathrm{3}$ km ($x\times y\times z$), with a horizontal spatial resolution of $\mathrm{\Delta }x=\mathrm{\Delta }y=\mathrm{100}$ m and a constant vertical grid spacing of Δz=40 m. These grid resolutions were selected such that individual cumulus clouds are appropriately resolved; however, it is expected to exhibit vertical artifacts near the inversion due to unresolved entrainment processes at these resolutions (Stevens et al., 2002). Additional information about the model setup for these LES cases can be found in Fridlind and Ackerman (2011) and Zhang et al. (2012). For each LES case, snapshots of cloud microphysical properties are taken every half hour after the first 4 h of each simulation, resulting in 9 cloud scenes per LES case and a total of 24 cloud fields.

The LES cloud fields comprise of spatially inhomogeneous microphysical properties, and these cloud properties change as the LES cloud field evolves at each time step. Figure 1 provides a map of the τ derived from the LES droplet size distributions for the ATEX clean, control and polluted cases at 4.0 h of simulation time. The map shows that the scenes are typically characterized by broken clouds with cloud fraction (CF; defined as a fraction of columns with τ>0.3) greater than 70 % in each LES case (CF of 73.7 %, 77. 2 % and 73.3 % for the ATEX clean, control and polluted cases respectively).

Figure 1 Map of the Large Eddy Simulation (LES) cloud optical thickness (τ) for the (a) ATEX clean, (b) ATEX control, and (c) ATEX polluted case at the 4 h simulation time.

2.2 Radiative Transfer Setup

The spherical harmonics discrete ordinate method (SHDOM) RT model developed by Evans (1998) was utilized to model reflectance for both 3D and 1D RT at the native LES resolution of 100 m. The RT simulations were performed for 6 solar zenith angles (SZAs) spanning from 10 to 60° in increments of 10°. A fixed nadir viewing zenith angle (VZA), and a constant relative azimuthal angle (ΔΦ=30°) were used throughout the study. The simulations for cloud retrievals are performed based on MODIS spectral response function for MODIS Band 2 (0.86 µm), Band 6 (1.6 µm), Band 7 (2.1 µm), and Band 20 (3.7 µm). We hereafter refer to these bands by their respective nominal wavelengths. It is important to note that although the 3.7 µm channel normally comprise of both solar and thermal portion and recent study by Loveridge and Di Girolamo (2024) suggest that contribution of the thermal emission to the error budget in the 3.7 µm retrieved r_e can sometimes be significant – partly due to their analysis utilizing partly cloudy pixels. This work assumes that the thermal emission in the 3.7 µm band can be removed with no error and thus have used only the solar portion in our RT calculations. For all RT simulations in this study, the surface was assumed to be Lambertian with an albedo of 5 % (an assumption approximating ocean surfaces under diffuse illumination), and double periodic horizontal boundary conditions were applied.

Observational studies usually estimate cloud susceptibility from derived N_d and broadband flux observations from Clouds and the Earth's Radiant Energy System (CERES) (e.g Painemal and Minnis, 2012). The CERES TOA flux is obtained from directional radiance measurements using scene-dependent anisotropy factors (often described via angular distribution models, ADMs), and therefore the CERES radiance to flux procedure explicitly accounts for the angular redistribution of radiation. In this study, rather than approximating fluxes from a single viewing direction or adopting an ADM assumption at LES scales, we compute the shortwave flux and albedo directly from radiative transfer by hemispherically integrating the upwelling radiance (I) at the top of the RT domain over a discrete angular grid. The integration is performed using radiance simulations at 37 viewing zenith angles (θ_v) spanning 0–90° in 2.5° increments and 16 viewing azimuth angles (ϕ_v) spanning 0–360° in 22.5° increments. The upwelling flux (F^↑) at solar zenith angle θ_o is defined as:

$$\begin{array}{}\text{(1)}& F^{\uparrow }\left({\mathit{\theta }}_{\mathrm{o}}\right)=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{\mathrm{0}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}I\left({\mathit{\theta }}_{\mathrm{v}},\mathit{\varphi };{\mathit{\theta }}_{\mathrm{o}}\right)\cos {\mathit{\theta }}_{\mathrm{v}}\sin {\mathit{\theta }}_{\mathrm{v}}\mathrm{d}{\mathit{\theta }}_{\mathrm{v}}\mathrm{d}\mathit{\varphi },\end{array}$$

where ϕ is the azimuth integration variable, $\mathit{\varphi }={\mathit{\varphi }}_{\mathrm{v}}-{\mathit{\varphi }}_{\mathrm{o}}$ is the relative azimuth angle, with ϕ_o and ϕ_v denoting the viewing and solar azimuth angles, respectively.

To ensure accurate flux and albedo estimates from the discrete angular sampling, the hemispheric integration was evaluated using finite-volume angular quadrature, in which each viewing-zenith bin was weighted by the exact integral of cos θ_vsin θ_v over its bin edges and each azimuth bin was weighted by its corresponding Δϕ. The resulting hemispherically integrated upwelling flux was then normalized by the incident solar flux at the top of the domain to obtain the broadband albedo (α), defined as:

$$\begin{array}{}\text{(2)}& \mathit{\alpha }\left({\mathit{\theta }}_{\mathrm{o}}\right)={\displaystyle \frac{F^{\uparrow }\left({\mathit{\theta }}_{\mathrm{o}}\right)}{F_{\mathrm{o}}{\mathit{\mu }}_{\mathrm{o}}}}\end{array}$$

where F_o is the incident solar flux at domain top and μ_o=cos θ_o. The broadband shortwave (0.3–5 µm) radiances used in the hemispheric integration are calculated for 13 out of the 14 Rapid Radiative Transfer Model (RRTM) spectral bands in the shortwave. Notably, we have excluded RRTM band 28 (0.2–0.26 µm) from these calculations because our Mie property computations failed for large droplet size bins (large size parameter). This exclusion is unlikely to impact the results significantly because this band contributes very little to the solar shortwave radiative energy budget. Unless otherwise stated, all albedo and susceptibility results presented in this study are based on these flux-derived broadband quantities computed consistently for both 1D and 3D forward RT simulations. Within a 1D RT framework, the upwelling flux and corresponding albedo can be interpreted locally because each column is radiatively independent from surrounding columns. However, within a 3D RT framework the horizontally resolved albedo field is not a uniquely local cloud-column quantity. The value of α calculated for a particular pixel depends on the height or reference at which the upwelling irradiance is registered (for our study the reference is at the top of the RT domain), because horizontal photon transport and angular integration redistribute radiative contributions from neighboring columns. Thus, different registration levels can produce different spatial variances of α and different covariances between α and local cloud properties. For example, if the irradiance field is evaluated at a level far above the cloud relative to the horizontal domain size, the locally registered albedo field can become increasingly smooth and approach a spatially uniform value, even though the area-averaged reflected flux remains physically well defined. Therefore, the local α used in the regression-based susceptibility calculations in this study should be interpreted as a domain-top, locally registered albedo diagnostic from the forward RT calculation, not as a unique column-resolved TOA albedo or as the exact 3D RT albedo susceptibility. The physically invariant 3D susceptibility is instead the response of area-averaged albedo to perturbations in cloud optical or microphysical properties, which we discuss in Sect. 3.4.

2.3 Bi-Spectral Retrieval Method

The bi-spectral retrieval method (Nakajima and King, 1990) described in Sect. 1 was applied to the simulated reflectance (see Sect. 2.3 in Ademakinwa et al., 2024). This method relies exclusively on homogeneous 1D RT assumptions to interpret the observed cloud reflectance. Its implementation is done using a precomputed lookup table (LUT), which includes computed 1D spectral reflectance for different τ and r_e combinations and solar-view geometries. The LUT then is used to find the τ and r_e whose computed reflectance best matches the observed (here, modeled LES) cloud reflectance. Notably, for small τ, retrieval uncertainty increases because the isolines of the LUT are less orthogonal and more tightly packed. This LUT non-orthogonality has significant consequences for observations with high spatial inhomogeneity below the pixel level resolution (Zhang et al., 2012, 2016) leading to the plane parallel homogeneous approximation (PPHA) bias. The consequences of both unresolved inhomogeneity and 3D radiative effects changes with spatial resolution and are examined in Sect. 3.3.

The VNIR reflectances used for the bi-spectral retrievals were calculated for 0.86 µm, while the SWIR and MWIR reflectances were calculate for the 1.6, 2.1 and 3.7 µm channel. We use a simplified cloud masking method that makes use of the 0.86 µm reflectance (i.e., R (0.86 µm) >0.07) to identify cloudy pixels after radiative transfer simulations. The LUTs for this study have 73 effective radii spanning from 4 to 40 µm and 105 log-spaced τ values spanning from 0.1 to 150. A constant effective variance (v_e) value of 0.1 in the gamma representation of the droplet size distribution is applied, consistent with operational MODIS retrievals. For consistency, the same surface albedo of 5 % used in all our LES RT reflectance simulations is applied in the LUT reflectance simulations.

We implement retrievals using reflectance simulations at the native LES resolution (100 m) and spatially average to a coarser MODIS-like resolution. For the MODIS-like resolution retrievals, we averaged reflectance from the LES resolution of 100 to 800 m resolution and utilized this area-averaged reflectance as input into the LUT for retrievals. Note that we have used 800 m resolution instead of the MODIS nadir 1 km because it divides evenly into our LES domain. Also, sensitivity analyses (not presented) suggest minimal differences between τ and r_e retrievals at 800 m and 1 km resolutions.

When averaging from high resolution to coarse resolution, we distinguish between overcast and partially cloudy pixels, depending on whether or not clear-sky pixels are included in the average reflectance. Firstly, if one or more pixels used for the area average is a clear-sky pixel, we classify the resulting pixel as a “partially cloudy pixel”. Secondly, if all pixels utilized for the area average are cloudy, we classify the averaged resulting pixel as an “overcast cloudy pixel”. With this, we define a category “all cloudy pixels” to consist of both partially cloudy and overcast cloudy pixels.

Cloud property retrievals that utilize area-averaged reflectances are impacted by unresolved (i.e., sub-pixel) spatial cloud inhomogeneity that can bias retrieval results, especially if the average is done over a highly inhomogeneous τ region. For example, Zhang et al. (2012) demonstrated using RT simulations how the nonlinearity in the LUT space (non-orthogonality of the τ and r_e LUT grid) can lead to plane-parallel albedo bias (Cahalan et al., 1994; the retrieved τ from the average reflectance of inhomogeneous cloud pixels tends to be smaller than the average of the sub-pixel τ) and plane-parallel r_e bias (r_e retrieved from area-averaged reflectances over inhomogeneous τ region is overestimated compared the original r_e).

2.4 N_d and LWP Estimates

2.4.1 N_d and LWP from LES

Cloud vertical structures in realistic cloud fields, such as the LES cases considered in this study, have microphysical properties (e.g., N_d and LWC) that vary vertically within each profile. Therefore, describing a single representative microphysical property value as a reference in the LES, requires accounting for the vertical distribution of the profile and on the application of interest. For example, most ACI and in situ studies typically take the reference N_d as an average value of N_d across the cloud vertical extent from cloud base to cloud top, or cloud-base N_d (e.g., Gryspeerdt et al., 2022; Painemal and Zuidema, 2011). In contrast, we adopt a different definition of N_d due to our application of interest: this study evaluates 3D effects biases in passive, radiance-based bi-spectral retrievals and examines how these biases propagate into N_d and albedo susceptibility estimates. Because passive retrievals and reflected SW radiative flux are most sensitive to the optically active portions of the cloud, a simple vertical mean does not necessarily provide the most radiatively relevant LES reference. We therefore define a radiation-relevant reference droplet number concentration in the LES as an extinction-weighted vertical mean, denoted N_{d_LES}. This is achieved by weighting N_d over the extinction coefficient (β_ext) at each layer (z) as given by:

$$\begin{array}{}\text{(3)}& N_{\mathrm{d}\mathrm{\_}\mathrm{LES}}={\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\underset{\mathrm{CTH}}{\overset{\mathrm{CBH}}{\int }}{N}_{\mathrm{d}}\left(z\right){\mathit{\beta }}_{\mathrm{ext}}\left(z\right)\mathrm{d}z,\end{array}$$

where $\mathit{\tau }={\int }_{\mathrm{CTH}}^{\mathrm{CBH}}\left({\mathit{\beta }}_{\mathrm{ext}}\right)\mathrm{d}z$. With this definition, layers with higher β_ext contribute more to the weighted value because they play a larger role in how the cloud interacts with scattered radiation, while reducing the influence of optically tenuous layers near cloud boundaries. We emphasize that N_{d_LES} is not intended to replace the vertically averaged or cloud-base N_d definitions commonly used in in situ and field-validation studies. Rather, it is used here as a radiatively weighted LES reference that is more directly aligned with the passive retrievals and broadband albedo analyzed in this study. Therefore, agreement between retrieval-derived N_d and N_{d_LES} should be interpreted as agreement with a radiatively effective cloud-column quantity, not necessarily as agreement with the vertically averaged or cloud-base N_d. This choice limits the theoretical interpretation of the retrieved N_d as a purely microphysical quantity, but it provides a more appropriate reference for evaluating how retrieval errors propagate into albedo susceptibility estimates. Also, the column LWP from the LES (LWP_LES) is obtained by integrating vertically the cloud liquid water content (LWC) from CBH to CTH in each column, expressed as:

$$\begin{array}{}\text{(4)}& {\mathrm{LWP}}_{\mathrm{LES}}=\underset{\mathrm{CTH}}{\overset{\mathrm{CBH}}{\int }}\mathrm{LWC}\left(z\right)\mathrm{d}z.\end{array}$$

2.4.2 Satellite-Based N_d and LWP Retrievals

Following previous studies, we calculate N_d from bi-spectral retrievals of τ and r_e following a pseudo-adiabatic model which assumes that (i) the LWC increases linearly (as a fixed fraction of its adiabatic value) with the cloud geometrical height and that (ii) N_d is constant vertically. Here, N_d is calculated from τ and r_e (denoted as N_{d_cal}) according to Grosvenor et al. (2018) which consolidates on prior effort by several previous studies (e.g., Brenguier et al., 2000; Quaas et al., 2006; Boers et al., 2006). Thus, we follow Grosvenor et al. (2018) simplified equation to calculate N_d from satellite-based retrievals given by:

$$\begin{array}{}\text{(5)}& N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}={\displaystyle \frac{\sqrt{\mathrm{5}}}{\mathrm{2}k\mathit{\pi }}}{\left({\displaystyle \frac{f_{\mathrm{ad}}{c}_{\mathrm{w}}\mathit{\tau }}{Q_{\mathrm{ext}}{\mathit{\rho }}_{\mathrm{w}}{r}_{\mathrm{e}}^{\mathrm{5}}}}\right)}^{\mathrm{1}/\mathrm{2}}\end{array}$$

where f_ad is the pseudo-adiabatic factor (constrained between 0 and 1), c_w is the effective condensation rate (g m⁻⁴), Q_ext is the bulk extinction efficiency derived from Mie theory (approximated as 2 in the geometric scattering limit), and ρ_w is the density of liquid water (taken as 1 g cm⁻³) and the parameter k is defined as:

$$\begin{array}{}\text{(6)}& k={\left({\displaystyle \frac{r_{\mathrm{v}}}{r_{\mathrm{e}}}}\right)}^{\mathrm{3}}\end{array}$$

where r_v is the droplet mean volume radius. Several values of k have been applied in N_d related studies and in operational remote-sensing that derive N_d from τ and r_e retrieved by the bi-spectral method. For example, the MODIS N_d calculation assume a constant value of k=0.72 (Grosvenor et al., 2018). Brenguier et al. (2011) utilized in situ measurements collected during five distinct field experiments to show that k values vary between 0.7 and 0.9, with uncertainties ranging from 10 % to 14 % across different cloud types and atmospheric conditions. All these suggests different values of k depending on the cloud regime. Studies based on in situ measurements have shown that k can vary significantly (Martin et al., 1994) and may also change under particularly extreme aerosol conditions (Noone et al., 2000). Thus, if k is not properly represented, it can introduce significant uncertainty into the calculation of N_d. To avoid these unconstrained constants driving our understanding of retrieval behavior, throughout this work we have calculated both k, and the quasi-adiabatic lapse rate (f_adc_w) in each cloudy column directly from the LES cloud field.

Several studies that derive N_d from bi-spectral retrieved cloud properties commonly utilize filtering techniques based on several criteria to minimize uncertainties in N_d estimates, which arise from from errors in retrieving r_e and τ from the bi-spectral method, especially in highly variable cloud fields (Grosvenor et al., 2018; Gryspeerdt et al., 2022; Loveridge and Di Girolamo, 2024; Quaas et al., 2006). However, this filtering methodology can unintentionally exclude certain types of clouds, limiting the analysis to specific cloud regimes. To avoid this bias, and to effectively quantify all possible sources of uncertainty, we incorporate all available N_d data without applying any pre-filtering method beyond the cloud mask. Furthermore, prior filtering schemes (e.g., Gryspeerdt et al., 2022; Zhu et al., 2018) implemented for coarser-resolution satellite observations may not be appropriate for the native LES resolution (100 m). For analysis involving LWP, we obtain the derived LWP (LWP_cal) from the retrieved τ and r_e by the adiabatic relationship given as (Wood and Hartmann, 2006):

$$\begin{array}{}\text{(7)}& {\mathrm{LWP}}_{\mathrm{cal}}={\displaystyle \frac{\mathrm{5}}{\mathrm{9}}}{\mathit{\rho }}_{\mathrm{w}}{r}_{\mathrm{e}}\mathit{\tau }\end{array}$$

3 Results and discussion

3.1 Comparison of LES N_d and N_d Derived from Cloud Properties Retrieved from 1D RT Simulated Reflectance

We first perform a RT closure comparison between the reference N_d obtained directly from the LES (N_{d_LES}) and the N_d derived from Eq. (5) using τ and r_e retrieved from homogeneous plane parallel RT-based simulated reflectance (i.e., 1D RT) from the 0.86 µm channel paired with an absorbing wavelength channel λ_abs, $\left(N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\right({\mathit{\lambda }}_{\mathrm{abs}}\left)\right)$. This closure comparison allows us to (i) check the methodology of our retrieval algorithm and (ii) investigate N_d errors associated with adiabatic cloud model assumptions. Figure 2 shows these comparisons as a regression plot of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ vs. N_{d_LES} for the ATEX clean, control and polluted cases at 4.0 h of simulation time, when λ_abs is 1.6, 2.1 and 3.7 µm. The comparison shows a general good agreement between $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ and N_{d_LES} for the ATEX clean, ATEX control and to some extent the ATEX polluted case with correlation coefficient (R) of 0.847, 0.852 and 0.748, respectively, when λ_abs is 1.6 µm; 0.885, 0.867 and 0.748 when λ_abs is 2.1 µm; and 0.859, 0.856 and 0.790 when λ_abs is 3.7 µm. The disparities observed between $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ and N_{d_LES} reflect the combined influence of vertical microphysical inhomogeneity, extinction weighting, and the assumptions embedded in the adiabatic N_d retrieval equation (Eq. 5). Therefore, agreement with N_{d_LES} indicates consistency with the optically dominant portions of the cloud column sampled by the passive retrieval.

Figure 2 Comparison of N_d derived from bi-spectral retrievals based on 1D RT reflectance ($N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$) versus N_d obtained from the LES (N_{d_LES}) for (a) ATEX clean, (b) ATEX control and (c) ATEX polluted for retrievals at absorbing channel, (λ_abs) of 1.6 µm, for (d) ATEX clean, (e) ATEX control and (f) ATEX polluted for retrievals at λ_abs of 2.1 µm and for (g) ATEX clean, (h) ATEX control and (i) ATEX polluted for retrievals at λ_abs of 3.7 µm at 4 h of simulation time. Note: 1D RT and bi-spectral retrievals utilized for this comparison were carried out at SZA 50°. Dotted lines indicate 1:1 relationship between the derived and LES N_d.

Establishing the underlying physical processes influencing the N_d estimates obtained from the LES and retrievals are challenging when relying solely on direct one-to-one comparison. Instead, we evaluate the probability density functions (PDFs) of N_{d_LES} against those obtained from the derived N_d. For the PDF comparison, data across all time steps within each LES case are combined to ensure robust statistical analysis. Apart from evaluating the PDFs of N_{d_LES} and $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$, we also investigate errors associated with the adiabatic cloud model assumptions in Eq. (5). This is done by introducing constraints on the retrieved r_e and τ used to compute N_d, which increasingly force conformity with the level of adiabaticity in the LES itself. First, we utilize the vertically weighted (VW) r_e (denoted as r_e(VW)), computed using the two-way-transmittance-weighted effective radius relationship derived from the optically weighted droplet size distribution defined by Miller et al. (2016) (check their Eq. 14), and optical thickness (τ_tot) from the LES into Eq. (5) to derive N_d (denoted as $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{r_{\mathrm{e}}\left(\mathrm{VW}\right),{\mathit{\tau }}_{\mathrm{tot}}}$). Second, we use droplet effective radius of the layer where the LWC is maximum as a representative of the adiabatic cloud top r_e (denoted as $r_{\mathrm{e}}^{\ast }$) and τ_tot as inputs in Eq. (5) to derive N_d (denoted as $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{{r_{\mathrm{e}}}^{\ast },{\mathit{\tau }}_{\mathrm{tot}}}$). The idea of testing with $r_{\mathrm{e}}^{\ast }$ is because it is representative of the droplet radius profile that follows the adiabatic assumption of the retrieval – absent significant entrainment modification at cloud top or model resolution artifacts. PDFs of N_d obtained from the LES, those inferred from the retrievals obtained from 1D RT simulated reflectance (at SZA 50°) when λ_abs is 1.6, 2.1 and 3.7 µm, and those derived from retrievals when the optical thickness and droplet effective radius are constrained as earlier discussed, for the ATEX clean, ATEX control, and ATEX polluted cases are presented in Fig. 3. In the ATEX clean and control cases (Fig. 3a and b), the N_{d_LES} distributions peak around their respective aerosol loading values (i.e., peak around 40 and 75 cm⁻³ for the ATEX clean and control cases respectively). This indicates that most of the available CCN are activated into cloud droplets. On the other hand, the lower values of the N_{d_LES} observed in the distribution are attributed to increased collision-coalescence as well as N_d removal processes driven by entrainment and precipitation.

Figure 3 Probability Density Functions (PDFs) of N_d obtained from the LES (N_{d_LES}) and N_d calculated from the adiabatic equation using 1D RT derived cloud properties (for absorbing channels 1.6, 2.1 and 3.7 µm), and under different levels of controls across all time steps for the ATEX clean, control and polluted cases in (a), (b) and (c) respectively.

Interestingly, the PDF distribution of N_{d_LES} for the ATEX polluted case (black line in Fig. 3c) is somewhat different from those of the clean and control cases. Its PDF has a peak around 300 to 350 cm⁻³, which is significantly less than the initial aerosol loading value (CCN 600 cm⁻³). This is probably because not all the CCNs are activated in the ATEX polluted case. When we consider the PDF of the $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ (for the three absorbing channels considered in this study) for the ATEX clean and control cases, similar cloud processes to the N_{d_LES} are observed, but with PDFs having a broader range and higher N_d peak values, and notably long tails towards larger values, indicating overestimation. The N_d overestimation statistics is more pronounced in $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}$(1.6 µm), reduces for $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}$(2.1 µm) and is smallest for $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}$(3.7 µm). This suggests that $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}$(3.7 µm) agrees the most with N_{d_LES} compared to results from the other two absorbing channels. When we consider the $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ PDFs from the ATEX polluted case (blue, green and purple lines in Fig. 3c), the peaks for all three absorbing channels are at smaller values than that of the N_{d_LES}, although all have long tails towards larger values indicating significant overestimation remains. For the ATEX polluted case, the $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}$(3.7 µm) PDF still has the smallest overestimation compared to the $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}$(1.6 µm) and $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}$(2.1 µm) results.

In general, the disparity in $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ obtained for all λ_abs (1.6, 2.1 and 3.7 µm) is due primarily to differences in the retrieved r_e – in vertically inhomogeneous clouds, satellite-based r_e retrievals via the bi-spectral method can vary for different channels due to spectrally-varying liquid water absorption (Meyer et al., 2025) leading to different penetration depths within a cloud (Platnick, 2000). The 3.7 µm being the most absorptive and 1.6 µm being the least absorptive of the three channels. For ideal clouds, the 3.7 µm channel is most absorptive and yield r_e retrievals closer to the cloud top compared to the less absorbing 2.1 µm channel, and the 1.6 µm channel is the least absorptive and yields retrievals deepest into the cloud. Since the N_d equation (Eq. 5) is most sensitive to r_e (r_e has the largest exponent in Eq. (5); $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}\mathit{\alpha }{r}_{\mathrm{e}}^{-\mathrm{5}/\mathrm{2}}$), small changes in r_e will greatly impact the derived N_d. Thus, the reasonable agreement observed between the PDF of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}$(3.7 µm) and N_{d_LES} is not surprising since Eq. (5) assumes r_e is at cloud top. This is corroborated by previous studies (e.g., Zhang and Platnick, 2011) that indicate that in adiabatic clouds, the 3.7 µm r_e retrieval is expected to be larger than retrievals from shorter wavelengths (and 2.1 µm r_e retrieval >1.6 µm r_e retrieval), although this relationship can be affected by other factors including retrieval biases.

In general, the $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ results in Fig. 3 from the three absorbing channels indicate that, even when the RT makes the same plane-parallel assumptions as the retrievals, there are notable variabilities in the derived $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ when compared to N_{d_LES}. These variabilities could be due to various reasons, which include inconsistencies between the vertical distribution of cloud microphysics in the LES truth and the assumptions that form the basis of the N_d equation utilized for the retrievals. Recall, the N_d derived from Eq. (5) assumes that the cloud is adiabatic i.e., LWC increases monotonically with height from cloud base to cloud top and N_d is constant vertically. But not all columns in the LES are adiabatic. Various processes such as entrainment, coalescence, etc. which are prevalent in natural clouds (and evident in this study LES cloud fields) introduce sub- or super-adiabatic behaviors in cloud vertical profiles, which impacts derived N_d.

When we consider the PDF of calculated N_d when τ and r_e are constrained to various degrees using values derived directly from the LES cloud field, the $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{r_{\mathrm{e}}\left(\mathrm{VW}\right),{\mathit{\tau }}_{\mathrm{tot}}}$ PDF shows substantial deviations from the PDF of N_{d_LES}. A reason for this is that r_e(VW) represents a weighted vertical average which is highly sensitive to the cloud microphysics and vertical structure of r_e, especially if the optical extinction in the cloud entrainment region is large enough to impact the r_e(VW). The agreement between the PDF's of N_d derived under the τ and r_e constraints (i.e., $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{{r_{\mathrm{e}}}^{\ast },{\mathit{\tau }}_{\mathrm{tot}}}$ and $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{r_{\mathrm{e}}\left(\mathrm{VW}\right),{\mathit{\tau }}_{\mathrm{tot}}}$) and N_{d_LES} reduces as the adiabatic realism of the input properties decreases; $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{{r_{\mathrm{e}}}^{\ast },{\mathit{\tau }}_{\mathrm{tot}}}$ showing better agreement with N_{d_LES} than $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{r_{\mathrm{e}}\left(VW\right),{\mathit{\tau }}_{\mathrm{tot}}}$ and N_{d_LES} comparison because $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{{r_{\mathrm{e}}}^{\ast },{\mathit{\tau }}_{\mathrm{tot}}}$ utilizes $r_{\mathrm{e}}^{\ast }$ which is most representative of an adiabatic cloud top r_e value which is a key requirement for formulation of Eq. (5). A sensitivity test was also performed where N_d is calculated using $r_{\mathrm{e}}^{\ast }$ and cloud optical depth computed from the LES that excludes cloud optical depth contributions from drizzle drop sizes (r) by isolating r>30 µm from the droplet size distributions. The resulting N_{d_cal} PDFs were mostly unchanged from the PDFs of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{{r_{\mathrm{e}}}^{\ast },{\mathit{\tau }}_{\mathrm{tot}}}$ (not shown). This was expected because drizzle-size drops typically occur in small numbers and contribute weakly to N_d, although their radiative impact could become important in cases with substantial drizzle populations.

3.2 Impact of 3D Radiative Effects on Retrieval-Derived N_d at the Native LES Resolution

To investigate the 3D radiative effects impacts on the retrieval-derived N_d, we compare $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ with N_d derived from Eq. (5) using τ and r_e retrieved from 3D RT-based simulated reflectance at λ_abs ($N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$). We carry out this comparison at different SZAs, ranging from high to low sun positions. Figure 4 shows the PDF of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ and $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ derived at the LES resolution (100 m) for a high sun (SZA 10°), moderate sun (SZA 30°) and a more oblique sun position (SZA 50°), as well as the reference N_d (N_{d_LES}) across all time steps (combined) in each of the ATEX clean, control and polluted LES cases. Here, N_{d_LES} PDF provides the LES-based reference distribution, so shifts or broadening in the $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ and $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ PDFs relative to N_{d_LES} indicate how retrieval assumptions and 3D radiative effects alter the inferred droplet number concentration.

When the sun is high, the PDFs of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ are generally similar to their $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ counterparts (where λ_abs=1.6, 2.1, and 3.7 µm), especially for the λ_abs=2.1 µm retrieval (Fig. 4a, d, and g). A reason for this similarity is the minimal 3D effects (lack of extreme darkening and brightening in the simulated 3D RT reflectance) when the sun is high, leading to the reflectance utilized for the cloud property retrieval under 3D RT at high sun comparable to its 1D RT counterpart. This in turn yields comparable derived N_d. Results for SZA 30°, show that the PDF of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ broadens slightly at both ends (Fig. 4b, e and h) compared to the high sun case (Fig. 4a, d, and g). This slight broadening of the $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ PDFs occurs because 3D-induced darkening and brightening effects occurs simultaneously. The darkening effects typically lead to retrieval of smaller τ and larger r_e, which together produce smaller $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ values compared to $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$.

Figure 4 Probability Density Functions (PDFs) of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ and $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ at native LES 100 m resolution at SZA 10, 30 and 50° and the reference LES N_d (N_{d_LES}) for the ATEX clean, control and polluted case at all time steps. Where λ_abs is the absorbing channel paired with 0.86 µm in the bi-spectral retrieval.

Also, the PDF of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ extends to larger N_d values, compared to $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$. These larger $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ values correspond to where τ is overestimated and r_e is underestimated (a result of the brightening effects). Overall, though, it appears that the shadowing effect is predominant in $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ for both high and moderate sun positions (discussed further in the domain-averaged plot in Fig. 5). When the sun is low, the PDF of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ become very broad, with larger peaks at smaller N_d, longer tails at large N_d, and a flattening of the peaks at moderate N_d that are observed in the high and moderate sun cases. This observed pattern is due to both brightening and darkening effects occurring simultaneously in the cloud property retrievals and contribute to the overestimate of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ (darkening effect) and underestimate of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ (brightening effect) compared to $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$. To examine the overall 3D radiative effect impact on the derived N_d at the LES domain scale, we calculate the in-cloud domain average of the retrieved N_d for each LES case at the different SZA's listed in Sect. 2 (SZA = [10°:60°] in steps of 10°). These domain-averaged $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ and $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ comparisons will provide insight on which 3D radiative effect (brightening or darkening) error is dominant (on the domain scale) in the cloud property retrievals (r_e and τ) that are utilized for the N_d calculations. Figure 5 shows the domain average of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ and $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ at the six SZA's for λ_abs of 1.6, 2.1, and 3.7 µm as well as domain average of the N_d reference for the ATEX clean, control and polluted cases averaged over all time steps.

Figure 5 Domain average of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ and $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ as a function of solar zenith angle as well as domain average of the N_{d_LES} reference for the ATEX clean, ATEX controlled and ATEX polluted case averaged across all time steps for λ_abs=1.6, 2.1, and 3.7 µm. The shaded region is the associated standard error computed from the population variance. Note that the y scale is different for each LES case due to the range of N_d associated with each LES case.

In Fig. 5, we generally observe that the domain-averaged $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ is smaller compare to the domain-averaged $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ under high to moderate sun positions. However, this reverses towards more oblique SZAs, with the domain-averaged $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ becoming larger than $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ when the SZA increases towards larger values. It is important to note that domain-averaged N_{d_LES}, computed directly from the LES microphysical fields and independent of SZA, absorbing channel, and RT retrieval assumptions is therefore used as the reference for assessing whether the retrieval-derived $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ and $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ are biased high or low at the domain scale. The separation between the red/blue lines and the green reference line shows that the magnitude and sign of the retrieval-derived N_d bias depend on SZA, absorbing channel, and LES case.

The physical impact of the 3D effects changes with the sun position: when the sun is high, absorption dominates horizontal transport differences between bands, but when the sun is low, horizontal transport is dominated by enhanced forward scattering across cloud boundaries. Thus our N_d result indicates that when the sun is high, horizontal transport is limited by absorption leading to darkening effects in the SWIR that dominate in the $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ domain-averaged statistics, whereas brightening effects in the VNIR becomes dominant when the sun is oblique. The SZA at which the domain-averaged $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ becomes larger than its 1D counterpart varies with the different scene configuration and absorbing channel utilized for the bi-spectral retrieval.

Interestingly, due to the strong absorption in the 3.7 µm band, in the ATEX polluted case (Fig. 5i), the $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}$(3.7 µm) is always larger than the $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}$(3.7 µm) for all SZAs considered in this study.

3.3 N_d Derived at at Coarse Spatial Resolution

Because the relative amount of horizontal energy transport between pixels is reduced as the pixel size becomes larger, the behavior of 3D radiative effects is spatial resolution dependent. So far, this study has focused on retrievals at the native LES resolution (100 m). At these high spatial resolutions, retrievals are more sensitive to 3D radiative effects, including cloud vertical inhomogeneity, and microphysical assumptions. However, at the coarser spatial resolutions common in global satellite imager retrievals (e.g., MODIS, Visible Infrared Imaging Radiometer Suite (VIIRS)), cloud-property retrievals are also affected by sub-pixel horizontal heterogeneity. This can introduce plane-parallel homogeneous approximation (PPHA) biases because the retrieval assumes each pixel is horizontally homogeneous, whereas the actual coarse pixel may contain unresolved variability in cloud optical thickness, cloud fraction, and cloud structure. Thus, at coarse resolution, the retrieved N_d can be influenced by both 3D radiative effects between neighboring cloudy regions and PPHA biases associated with unresolved sub-pixel heterogeneity. Therefore, extending the analysis of derived N_d to coarser satellite-like resolutions is necessary to determine whether the 3D-RT-induced retrieval errors observed at the native LES resolution persist after spatial aggregation and under resolutions more comparable to operational passive imager observations.

We examine the behavior of N_d derived from a MODIS or VIIRS-like resolution (800 m is convenient for our LES models) using retrievals from area-averaged 3D RT simulated reflectances, and then compare those N_d PDFs with derived N_d obtained from 3D RT simulations at the native LES resolution (100 m). The PDF comparison of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}$(2.1 µm) obtained at both LES and coarse 800 m resolution when all cloudy pixels (i.e., partially cloudy and overcast cloudy pixels; refer to Sect. 2.3 for definition) are considered, is shown in Fig. 6 and when only overcast cloudy pixels are considered is presented in Fig. 7. In Figs. 6 and 7, the green N_{d_LES} PDFs provides the LES reference distribution for each case and allows the 100 and 800 m retrieval-derived N_d PDFs to be interpreted relative to them. From Fig. 6, we clearly see that when all cloudy pixels are considered, the derived N_d PDF at the coarse 800 m resolution are distinct from those at the native LES resolution: the PDF of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}$(2.1 µm) at the coarse resolution is mostly skewed towards smaller N_d values, for the three solar geometries shown, with smaller in-cloud mean values for the coarse resolution compared to the corresponding LES resolution. The reason for this is due to the contribution of partially cloudy pixels in the coarse resolution retrievals, further enhancing the impact of the plane-parallel r_e and τ bias contribution and result in smaller N_d values.

Figure 6 Probability Density Functions (PDFs) comparison of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}$(2.1 µm) for LES (100 m) and coarse satellite-like (800 m) resolution for all cloudy pixels, at SZA 10, 30 and 50° for the ATEX clean, control and polluted case across all time steps. Vertical lines represent mean values. The green line represents the LES reference (N_{d_LES})

Figure 7 Probability Density Functions (PDFs) comparison of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}$(2.1 µm) for LES (100 m) and coarse satellite-like (800 m) resolution for overcast cloudy only pixels, at SZA 10, 30 and 50° for the ATEX clean, control and polluted case across all time steps. Vertical lines represent mean values. The green line represents the LES reference (N_{d_LES})

When only overcast cloudy pixels are considered (Fig. 7), the contribution of the partially cloudy pixels is removed and the smaller N_d values observed in Fig. 6 are mostly absent. As seen in Fig. 7, the PDFs of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}$(2.1 µm) at the coarse (800 m) resolution, are mostly centered around moderate N_d values, and perhaps surprisingly, happen to have mean values that closely match corresponding values at the LES (100 m) resolution compared to results which utilize all cloudy pixels in the retrievals. These results can provide some additional confidence in using N_d calculations at coarse resolution for our LES cloud fields, as the effects of cloud heterogeneity and opposing brightening and darkening 3D effects appear to cancel out (to a large extent) and still provide comparable mean values (although this depends on how well the satellite can identify partially cloudy pixels).

We also examined the sensitivity of $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ to spatial aggregation at the other MODIS absorbing bands utilized for r_e retrieval by comparing native LES-resolution (100 m) results with coarse satellite-like resolution (800 m) results for λ_abs=1.6 and 3.7 µm (not shown). In general, the coarse-resolution results based on overcast 800 m pixels are more consistent with the native-resolution mean $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}$ than results based on all cloudy 800 m pixels. This behavior is found for both absorbing channels, although for λ_abs=3.7 µm at SZA = 10°, the all-cloudy and overcast-only estimates are both comparable to the native-resolution result. These tests indicate that the treatment of partially cloudy coarse pixels can influence aggregated N_d, but the main resolution-dependent behavior remains robust across the two absorbing channels considered.

3.4 Implications for Regression-Based Albedo Susceptibility Estimates for Evaluating the Twomey Effect

Having investigated the impacts of 3D effects on the retrievals of τ and r_e, and the subsequently derived N_d, in this section we examine the consequent implications for estimating Twomey effect from these retrievals. In many previous studies, Twomey effect is usually estimated from observations using the so-called absolute cloud albedo susceptibility (S_α) that connects the change of cloud albedo (α) with the change of N_d at a given LWP (Ackerman et al., 2000; Platnick and Twomey, 1994),

$$\begin{array}{}\text{(8)}& S_{\mathit{\alpha }}\equiv {\left.{\displaystyle \frac{\partial \mathit{\alpha }}{\partial N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}}}\right|}_{\mathrm{LWP}}\approx {\displaystyle \frac{\partial \mathit{\alpha }}{\partial \mathit{\tau }}}\cdot {\left.{\displaystyle \frac{\partial \mathit{\tau }}{\partial N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}}}\right|}_{\mathrm{LWP}}.\end{array}$$

While S_α is often derived empirically through numerical regression of satellite retrievals, many researchers try to develop theoretical expectations for these relationships. This involves using the chain rule of differentiation (Eq. 8), where the first term – how cloud albedo depends on cloud properties (e.g., $\frac{\partial \mathit{\alpha }}{\partial \mathit{\tau }}$), is usually estimated from simplified radiative transfer models such as the two-stream approximation. The second term – how τ changes with respect to N_d at constant LWP ($\frac{\partial \mathit{\tau }}{\partial N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}}{\mathrm{|}}_{\mathrm{LWP}}$) is often derived from idealized cloud models assuming constant LWC or sub-adiabatic conditions. Of course, the change of N_d can also potentially lead to the change of LWP known as cloud adjustment, but it is beyond the scope of this study. Here, we first examine the impacts of 3D effect on the estimation of $\frac{\partial \mathit{\tau }}{\partial N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}}{\mathrm{|}}_{\mathrm{LWP}}$, which we define as the retrieval sensitivity s(τ). Because a constant N_d is assumed, regardless of the cloud vertical profile assumptions, the relationship $\mathit{\tau }\propto N_{\mathrm{d}}^{\mathrm{1}/\mathrm{3}}{\mathrm{LWP}}^{\mathrm{5}/\mathrm{6}}$ holds (Ackerman et al., 2000; Platnick and Twomey, 1994), and implies that,

$$\begin{array}{}\text{(9)}& s\left(\mathit{\tau }\right)\equiv {\left.{\displaystyle \frac{\partial \ln \mathit{\tau }}{\partial \ln N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}}}\right|}_{\mathrm{LWP}}\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}.\end{array}$$

Constraining fixed LWP regimes is required to evaluate s(τ) and S_α. Thus, we group LWP into eight bins ranges (B1 to B8) as given in Table 1. The bin edges were derived from percentile boundaries of valid retrieved LWP distribution computed at 10 % intervals; however, the first two percentile ranges characterized by small LWP were excluded, leaving the eight retained bins now labeled B1 to B8. We will point it out explicitly in the analysis whether the LWP bin categorization is based on the 1D or 3D retrievals or from the LES LWP (LWP_LES) cloud field.

Table 1 Liquid Water Path (LWP) bins used in this study. B1 to B8 are the bin number representation. Max(LWP) = 3333.33, is maximum LWP value obtained from the retrievals (computed by utilizing Max allowable r_e (40 µm) and Max allowable τ (150) in Eq. (7). The bin edges were derived from percentile boundaries of the valid retrieved LWP distribution computed at 10 % intervals. The first two percentile ranges were excluded, leaving the eight retained bins labeled B1–B8.

Describing retrieval sensitivity requires numerical regression of ∂ln τ with ∂ln N_d. As will be demonstrated in this section, 3D effects can influence the estimation of s(τ) though two main mechanisms:

• The errors in τ and N_{d_cal} retrievals caused by the 3D effects as shown in Sect. 3.2 and 3.3 can in turn lead to errors in the ∂ln τ–∂ln N_{d_cal} regression. We refer to this as the “retrieval error effect”.

• To obtain the theoretical $s\left(\mathit{\tau }\right)\approx \frac{\mathrm{1}}{\mathrm{3}}$, the ∂ln τ–∂ln N_{d_cal} regression must be performed at constant LWP. Since LWP is also derived from retrievals in reality, it is subject to errors caused by 3D effects. These errors can cause mis-categorization of LWP bins in the regression analysis. We refer to this as the “LWP categorization error”.

Figure 8 Regression plots of retrieved cloud optical thickness (τ) versus derived cloud droplet number concentration (N_d) for Liquid Water Path (LWP) bin 7 (B7; range 142.61–121.19 g m⁻²), based on bi-spectral retrievals using the 0.86 and 2.1 µm channels at a large-eddy simulation (LES) resolution of 100 m and a solar zenith angle (SZA) of 50°. Panel (a) shows results under 1D radiative transfer (RT) using LWP bin categorization based on 1D retrievals (i.e., LWP Self-Categorization); panel (b) shows 3D RT results using LWP binning from 3D retrievals (i.e., LWP Self-Categorization); and panel (c) presents 3D RT results using LWP bin categorization based on 1D retrievals. Values in bold print indicate sensitivity s(τ) computed according to Eq. (9) and R is the correlation coefficient.

We first demonstrate the impact of retrieval and LWP categorization errors by calculating s(τ) from regressions of a single LWP bin (B7; 142.61–212.19 g m⁻²). This analysis combines all three ATEX cases (clean CCN = 40 cm⁻³, control CCN = 75 cm⁻³, and polluted CCN = 600 cm⁻³) utilizing the 0.86 and 2.1 µm channels at SZA 50°. The regressions shown in Fig. 8 depict τ vs. N_{d_cal} and corresponding s(τ) obtained from regression. For both 1D and 3D retrievals based on LWP self-categorization (i.e., where LWP is determined consistently from cloud properties retrieved from the respective RT reflectance fields), the regressions in Fig. 8a and b clearly demonstrate a consistent power law scaling with different slopes; broadened variability in 3D retrievals, caused by 3D effects, which subsequently leads to a more consistent s(τ) value (s(τ)=0.334), even better than 1D RT result (s(τ)=0.349). However, in Fig. 8c, we observe that clustering 3D retrievals into LWP bins using a reference dataset that is insensitive to the 3D radiative effects (in this case the reference LWP is calculated from 1D retrievals) leads to noticeably different distributions, where regression becomes less well correlated (correlation coefficient (R)=0.857, compared to R>0.9 for both LWP Self-categorization cases) and alters the calculated s(τ) value (s(τ)=0.363). This mismatch between retrievals and LWP classification highlights the importance of using a self-consistent LWP binning approach when calculating sensitivity. Maintaining consistency in LWP categorization is also important for accurately analyzing susceptibility estimates; thus, the use of consistent LWP categorization will be explicitly stated wherever it is applied in the remainder of this study.

Having demonstrated the impact of 3D effects on s(τ) using a single representative LWP bin, we next use our simulations to investigate how these errors influence the cloud susceptibility. To evaluate this, it is useful to examine how albedo responds to variations in the relative droplet number concentration. Thus, we follow the approach of Painemal and Minnis (2012) and define relative albedo susceptibility (S_α,rel), which accounts for the dependence of droplet number concentration on spatial variability, as:

$$\begin{array}{}\text{(10)}& S_{\mathit{\alpha },\mathrm{rel}}=N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}{\left.{\displaystyle \frac{\partial \mathit{\alpha }}{\partial N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}}}\right|}_{\mathrm{LWP}}={\left.{\displaystyle \frac{\partial \mathit{\alpha }}{\partial \ln N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}}}\right|}_{\mathrm{LWP}}.\end{array}$$

Analyzing S_α,rel, which is based on fractional (logarithmic) changes in number concentration, will help to minimize the impact of absolute error biases in N_{d_cal} retrievals.

Before interpreting the susceptibility results, it is important to distinguish the exact 3D RT susceptibility from the regression-based albedo susceptibility estimated in this study. The physically invariant 3D albedo susceptibility is the functional derivative of an area-averaged albedo with respect to perturbations in the cloud optical or microphysical properties, which describes how the area-mean albedo responds to a perturbation in a particular cloud column or group of columns. Unlike a regression slope based on locally registered albedo, this derivative is not determined by the arbitrary height chosen to register the spatially varying flux/albedo field. Computing this exact quantity in 3D requires a linearized or adjoint 3D RT calculation, or an equivalent finite-difference perturbation framework. Such a calculation requires additional modeling and computational framework and is out of scope of this present study.

As stated in Sect. 2.2 and earlier in this section, the quantity evaluated here is an approximate, regression-based susceptibility diagnostic rather than an exact 3D RT susceptibility. Specifically, within each LWP bin, we estimate the slope of the locally registered, domain-top broadband albedo, α, with respect to ln N_d (Eq. 10). This approach is similar to regression-based observational albedo-susceptibility analyses, where TOA albedo from CERES and N_d inferred from passive cloud retrievals such as MODIS are used to estimate susceptibility (e.g., Painemal and Minnis, 2012). However, under 3D RT, the locally registered α field is sensitive to the chosen reference level because horizontal photon transport and angular integration redistribute reflected radiation across neighboring columns. As a result, the spatial variance of α and its covariance with N_d can depend on the albedo registration convention. Therefore, the susceptibility values reported in this study should not be interpreted as exact column-resolved 3D RT derivatives of TOA albedo. Rather, they quantify how 3D RT and retrieval errors affect regression-based albedo–N_d susceptibility estimates under the specific domain-top registration convention used in this study.

Our interpretation does not rely on the absolute value of the regression slope as a physically invariant susceptibility. Instead, we focus on how the slope changes when the same regression framework and albedo-registration convention are applied to albedo and N_d fields affected by different RT assumptions, spatial resolutions, and absorbing-channel choices. In this sense, the diagnostic quantifies how 3D RT effects would modify the albedo–N_d covariance and associated regression-inferred susceptibility under an observational-style regression framework, rather than the exact differential response of area-mean TOA albedo to a controlled perturbation in N_d.

Figure 9 Scatter plots for similar liquid water path (LWP) bins showing albedo (α) vs. N_d (logarithmic scale). α under 1D RT vs. derived N_d obtained from retrievals based on 1D RT reflectance when LWP bin is categorized by 1D retrievals (LWP Self Categorization) in (a), α under 1D RT vs. N_{d_LES} when LWP bin is categorized by the LES LWP in (b), α under 3D RT vs. derived N_d obtained from retrievals based on 3D RT reflectance when LWP bin is categorized by 3D retrievals (LWP Self Categorization) in (c), α under 3D RT vs. N_{d_LES} when LWP bin is categorized by the LES LWP in (d). Map of α for the ATEX clean, control and polluted LES scenes under 1D RT at 4 h of simulation time in (e), (f) and (g) respectively and under 3D RT in (h), (i) and (j) respectively. All RT simulations are carried out at SZA 20° and bi-spectral retrievals utilize the 0.86 and 2.1 µm channel at LES resolution of 100 m.

Figures 9 and 10a–d present the α–N_d relationships at the native LES resolution of 100 m across all LES cases and time steps, stratified by the eight LWP bins defined in Table 1. Figure 9 corresponds to SZA = 20°, while Fig. 10 shows the same set of diagnostics for a more oblique illumination condition at SZA = 60°. In both figures, panels (a) and (c) show α versus N_{d_cal}, where both N_{d_cal} and LWP are inferred from bi-spectral retrievals applied to 1D- and 3D-RT simulated reflectances, respectively, using the 0.86 and 2.1 µm wavelength pair. Panels (b) and (d) in both figures show α versus the LES reference droplet concentration (N_{d_LES}), with LWP binning based directly on the LWP_LES, and with α obtained from 1D and 3D RT, respectively. The corresponding albedo maps for the ATEX clean, control, and polluted cases at 4.0 h are shown in panels (e)–(g) for 1D RT and panels (h)–(j) for 3D RT. Together, Figs. 9a–d and 10a–d illustrate a key feature of the first aerosol indirect effect (Twomey, 1977): under fixed LWP conditions, increases in droplet number concentration are generally associated with increases in cloud optical thickness and albedo.

Figure 10 Scatter plots for similar liquid water path (LWP) bins showing albedo (α) vs. N_d (logarithmic scale). α under 1D RT vs. derived N_d obtained from retrievals based on 1D RT reflectance when LWP bin is categorized by 1D retrievals (LWP Self Categorization) in (a), α under 1D RT vs. N_{d_LES} when LWP bin is categorized by the LES LWP in (b), α under 3D RT vs. derived N_d obtained from retrievals based on 3D RT reflectance when LWP bin is categorized by 3D retrievals (LWP Self Categorization) in (c), α under 3D RT vs. N_{d_LES} when LWP bin is categorized by the LES LWP in (d). Map of α for the ATEX clean, control and polluted LES scenes under 1D RT at 4 h of simulation time in (e), (f) and (g) respectively and under 3D RT in (h), (i) and (j) respectively. All RT simulations are carried out at SZA 60° and bi-spectral retrievals utilize the 0.86 and 2.1 µm channel at LES resolution of 100 m.

For the retrieval-based comparisons, the α–N_{d_cal} relationships are well correlated across all LWP bins, with correlation coefficients greater than 0.77 in Fig. 9a and c, as well as in Fig. 10a and c. However, the rate at which α increases with N_{d_cal} differs among the cases and depends on whether the RT used to compute α is treated as 1D or 3D and also the RT method used to generate the reflectances from which N_{d_cal} is obtained. For 1D RT based results, the α(1D RT) vs. $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}$ relationships in Figs. 9a and 10a are very similar, showing increasing α with N_{d_cal} and mainly distinct across each LWP bins. When we consider the corresponding 3D RT cases in Figs. 9c and 10c, the α(3D RT) vs. $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}$ relationship still shows an overall increase in α(3D RT) with increasing $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}$. However, the relationship becomes more dispersed under the more oblique illumination condition. Specifically, Fig. 10c, corresponding to SZA = 60°, shows a larger spread in α and stronger overlap among LWP bins than Fig. 9c, corresponding to SZA = 20°. This indicates that 3D radiative effects become more pronounced under the more oblique SZA 60° case, as seen by the larger spread and stronger overlap among LWP bins in Fig. 10c (SZA = 60°) compared to Fig. 9c (SZA = 20°). This increased variability reflects the combined influence of 3D radiative effects on the retrieved τ and r_e, which subsequently affect $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}$ and derived LWP, together with the fact that α from 3D RT is less local than α from 1D RT. This nonlocality is illustrated by the albedo maps in Figs. 9e–j and 10e–j. For the ATEX clean, control, and polluted cases at 4.0 h, the 1D RT albedo fields retain sharper small-scale cloud structure, whereas the corresponding 3D RT albedo fields are smoother and more spatially blurred. This blurring is consistent with radiative smoothing caused by horizontal photon transport and multiple scattering, such that the top-of-domain albedo assigned to a given pixel contains contributions from surrounding cloud pixels. Consequently, the α–N_d relationship under 3D RT becomes less tightly controlled by the local cloud column properties alone, especially at larger SZA.

When we consider the plots of α(1D RT) versus N_{d_LES} in Figs. 9b and 10b, they show a pattern that is broadly similar to the 1D RT based retrieval results in Figs. 9a and 10a, respectively. In both SZA cases, α increases with increasing N_d, and the LWP bins remain relatively well separated. This similarity is expected because both panels (a) and (b) use α simulated under 1D RT, where the radiative response is largely controlled by the cloud and surface properties in the cloud column. Therefore, whether N_d, is taken from the retrievals in panels (a) or from the LES reference in panels (b), the α–N_d, relationship still retains a clear Twomey-like behavior. Panels (d) in Figs. 9 and 10 provide the clearest view of how 3D RT affects the α–N_d relationship when retrieval-related effects are removed, because they use α from 3D RT, the LES reference N_{d_LES}, and LWP_LES binning. Compared with panels (a) and (b), the relationship in panels (d) is less organized, with larger spread in α and stronger overlap among LWP bins. This indicates that even when N_{d_LES} and LWP are used, α from 3D RT is still influenced by surrounding cloud structure through horizontal photon transport and radiative smoothing, rather than by radiative properties of the local cloud column alone. The effect is more pronounced in Fig. 10d than in Fig. 9d, showing that the nonlocal influence of 3D RT increases under the more oblique illumination condition at SZA = 60°. Thus, panels (d) show that 3D RT can weaken the local Twomey-like α–N_d relationship, with direct implications for susceptibility estimates under 3D RT. Figure 11 shows regression-based S_α,rel from diagnostics in Fig. 9 (at SZA 20°) and Fig. 10 (at SZA 60°) computed within the same LWP bins B1 to B8 and plotted as a function of the mean LWP and mean τ of each bin. The LES-based benchmark is computed from the regression of domain-top α(3D RT) versus. N_{d_LES} relationship using the LWP_LES binning (the green curve). This is useful for evaluating the retrieval-based regression estimates because it uses the same forward 3D RT albedo field and LES microphysical information, while avoiding retrieval errors in N_d and LWP.

Figure 11 Mean relative cloud albedo susceptibility (S_α,rel) across all LES cases and time steps at solar zenith angles (SZAs) 20 and 60°, shown as a function of liquid water path (LWP) in panel (a), and as a function of cloud optical thickness (τ) in panels (b). Where, N_d and corresponding LWP or τ are retrieved by bi-spectral method using the reflectance pairs of 0.86 combined with 2.1 µm at the LES-native resolution (100 m). Green curves show the LES-based regression benchmark computed from domain-top 3D RT albedo, N_{d_LES}, and LWP_LES binning; these curves are used as diagnostic benchmarks and are not exact linearized 3D RT susceptibility derivatives.

For SZA = 20°, the 1D RT result follows the LES-based regression benchmark closely at low to moderate LWP/τ values, where S_α,rel increases and then begins to decrease. The 3D RT result is lower than the LES-based benchmark at the smallest LWP/τ, but becomes comparable at intermediate values and remains larger than the LES-based benchmark at higher LWP/τ, suggesting that 3D radiative effects modify the LWP/τ dependence of the susceptibility. For the more oblique case at SZA = 60°, all susceptibility values are generally smaller than those at SZA 20°, and the 3D RT retrieval produces a flatter response with increasing LWP compared with the LES-based benchmark and the 1D RT estimate. This indicates that, under more oblique illumination, 3D radiative smoothing and retrieval effects weaken the local α–N_d sensitivity within the LWP bins. The α–N_{d_cal} relationship obtained when N_{d_cal} is derived from cloud properties retrieved by pairing the 0.86 µm with either the 1.6 or 3.7 µm band exhibits a similar pattern to that observed with the 0.86 and 2.1 µm combination (not shown).

Figure 12 Mean relative cloud albedo susceptibility (S_α,rel) across all LES cases, time steps and six solar zenith angles (SZAs), shown as a function of liquid water path (LWP) in panels (a), (c), and (e), and as a function of cloud optical thickness (τ) in panels (b), (d) and (f). In each case, N_d and corresponding LWP or τ are retrieved by bi-spectral method using the reflectance pairs of 0.86 combined with 1.6 in (a) and (b), 2.1 in (c) and (d), and 3.7 µm in (e) and (f). Triangles represent S_α,rel derived from 3D RT results, while circles denote values from 1D RT. broken lines indicate result of 800 m coarse resolution, and solid lines represent LES-native resolution (100 m). Green curves show the LES-based regression benchmark computed from domain-top 3D RT albedo, N_{d_LES}, and LWP_LES binning; these curves are used as diagnostic benchmarks and are not exact linearized 3D RT susceptibility derivatives.

After examining the low- and high-SZA cases, we next evaluate whether the same susceptibility behavior persists across the full set of illumination conditions used in this study i.e., across all SZA's (SZA = [10°:60°] in 10° intervals). To achieve this, we calculate the mean relative cloud albedo susceptibility, S_α,rel, averaged across all LES cases, time steps, and all six SZAs. Figure 12 presents these mean S_α,rel as a function of mean LWP within each bin in panels (a), (c) and (e) and as a function of mean τ in panels (b), (d) and (f) for retrievals using the 0.86 µm band paired with 1.6, 2.1, and 3.7 µm, respectively. The LES-based regression benchmark curves are also included at both 100 and 800 m scales so that the retrieval-based susceptibility estimates can be compared against a consistently defined regression benchmark defined at the same spatial resolution. From Fig. 12, we observe that, overall, S_α,rel generally increases from the smallest LWP/τ bins to low-to-moderate LWP/τ values, and then decreases toward larger LWP/τ regimes. This behavior is consistent with the reduction in albedo sensitivity as clouds become optically thicker and albedo begins to saturate. The LES-based regression benchmark (green line) shows the same general pattern, indicating that the regression-based susceptibility estimated from retrievals, broadly capture the expected dependence of albedo sensitivity on clouds.

At the native LES resolution of 100 m, the 1D RT retrievals generally follow the LES-based regression benchmark more closely than the 3D RT retrievals, especially at low-to-moderate LWP/τ. The 3D RT retrievals exhibit smaller S_α,rel in the lower LWP/τ regime and show a more gradual decrease at larger LWP/τ. This difference reflects the combined influence of 3D radiative smoothing on α and 3D effects induced biases in the retrieved τ and r_e, which propagate into N_d and subsequent susceptibility estimate. At the coarser 800 m resolution, the differences between the 1D and 3D RT susceptibility estimates are reduced at the low LWP/τ regimes although notable differences occur at the larger LWP/τ values. The 800 m retrieval curves are also more comparable to the 800 m LES-based regression benchmark than would be expected from comparison with the 100 m reference alone, confirming that the susceptibility benchmark is resolution dependent. This indicates that spatial aggregation partly reduces fine-scale 3D variability and promotes compensation between 3D radiative effects, plane-parallel retrieval biases, and subpixel averaging.

The same broad behavior of S_α,rel is observed for all three absorbing-channel combinations at the 100 m resolution, with S_α,rel generally increasing at low LWP/τ before gradually decreasing towards larger LWP/τ values. At the coarse 800 m resolution, the S_α,rel vs. LWP/τ relationship is similar to those at 100 m, but the 1D vs. 3D differences is band dependent. The 1.6 and 2.1 µm retrieval results show somewhat better agreement between the 1D and 3D curves at low to moderate LWP/τ values. They all have significant differences between 1D and 3D at larger LWP/τ values. This indicates that 3D RT effects can greatly modify S_α,rel at native LES resolution, but their impact is substantially reduced at coarse satellite-like resolution in low LWP/τ regimes.

4 Summary and Conclusion

Errors associated with 3D radiative effects and cloud inhomogeneity can influence bi-spectral retrievals of cloud properties (τ and r_e), and subsequently impact estimates of N_d derived from these properties. Because N_d is so important in the representation of clouds in models and the interpretation of aerosol–cloud effects, it is important to disentangle and explain these sources of retrieval bias. Therefore, this study focuses on investigating the impact of the 3D radiative effects on N_d calculated from τ and r_e retrieved by the bi-spectral method (N_{d_cal}) and further probes how an understanding of regression-based albedo susceptibility which quantifies the first aerosol indirect effect is impacted when interpreted from such retrievals. We address this by a satellite observation-retrieval simulator framework, which consists of synthetic cloud fields, a radiative transfer solver and retrieval algorithm. The cloudy scenes examined here spanned numerous temporal snapshots from three LES cases, each of which was initialized based on observations made during the 1969 ATEX field campaign (NE Atlantic trade wind region ∼ 12° N, 35° W). These LES cases are each characterized by different initial aerosol loadings: a clean (CCN = 40 cm⁻³), control (CCN = 75 cm⁻³), and polluted (CCN = 600 cm⁻³) case. From these LES cloud fields, simulated reflectances (based on 1D and 3D RT) for MODIS Band 2 (0.86 µm), Band 6 (1.6 µm), Band 7 (2.1 µm), and Band 20 (3.7 µm) were carried out with the SHDOM RT model (Evans, 1998). Thereafter, bi-spectral retrievals were implemented on the simulated reflectance by pairing the 0.86 µm VNIR band with each of the absorbing wavelength bands (λ_abs=1.6, 2.1 and 3.7 µm). Subsequently, N_{d_cal} were obtained from the bi-spectral microphysical retrievals at six solar geometries (SZA = [10°:60°] in steps of 10°) and the results were analyzed.

We utilize N_d obtained from the LES as a reference and compare with N_{d_cal} derived from 1D RT simulations, and with N_{d_cal} inferred from the adiabatic assumption (Eq. 5) when τ and r_e retrievals are constrained toward scenes that more closely approximate that assumption. Results demonstrates that, N_{d_cal} from τ and r_e retrieved from homogeneous plane parallel RT-based simulated reflectance can still be biased when compared to the reference N_d. This bias arises from the complexity of how cloud microphysics is vertically distributed within realistic clouds, which may sometimes deviate from the adiabatic assumptions in which the N_{d_cal} equation is based. Also, the choice of the absorbing channel in the bi-spectral retrieval is a factor that contributes to this derived bias. Due to different liquid water absorption and therefore vertical sensitivities at different channels, the 3.7 µm band which has the strongest absorption, retrieves r_e closer to the cloud top, and agrees the most with the reference LES N_d, compared to results derived from pairing the 0.86 µm band with the 1.6 and 2.1 µm band.

To investigate how 3D radiative effects impact N_d estimates, we calculate N_{d_cal} using cloud properties derived from 1D RT-based reflectances and compare directly with N_{d_cal} calculated using cloud properties derived from 3D RT-based reflectances. These comparisons are performed for multiple bi-spectral band pairs (pairing VNIR = 0.86 µm with either of λ_abs=1.6, 2.1, and 3.7). Results from this comparison carried out at different SZAs show that, for high sun (SZA = 10°), the statistical distributions of N_{d_cal}(λ_abs) from 3D and 1D RT agree well for all three LES cases. Domain-averaged $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ and $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ was also compared to determine the overall impact of the 3D (brightening and darkening) effects on the domain scale. For the high sun case, the domain-averaged $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ is smaller compared to $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$. This indicates that the darkening effect which dominates cloud property retrievals when the sun is high also dominates the domain-averaged $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ result (i.e., overestimated r_e and underestimated τ dominates the cloud property retrievals and yields lower $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ values compared to $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ results). A similar pattern is also observed for the domain-averaged $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ derived from not too oblique SZAs (e.g., SZA 30°).

When the sun is low at SZA 50°, 3D radiative effects become stronger and produce larger variability in the retrieval-derived N_d. This occurs because both brightening and darkening effects occur in significant amount at the low sun angle. Brightening effects tend to produce larger τ and smaller r_e, leading to larger $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ compared to corresponding $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ values, while darkening effects tend to produce smaller τ and larger r_e, leading to smaller $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ compared to $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{1}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ values. At the domain scale, the relative magnitude of these effects changes with SZA and absorbing channel used for retrievals. For larger SZAs, the brightening effect becomes more dominant in several cases, causing the domain-averaged $N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}^{\mathrm{3}\mathrm{D}}\left({\mathit{\lambda }}_{\mathrm{abs}}\right)$ to exceed its 1D RT counterpart, although the SZA at which this transition occurs depends on the absorbing channel and cloud scene.

The statistical distributions of N_{d_cal} at satellite-like (MODIS, VIIRS) pixel footprints were also examined. Compared to retrievals at the native LES resolution of 100 m, the 800 m results exhibit some subtle differences depending on how one defines the cloud mask for the coarse resolution observations. We divide the retrieval population into two groups: all cloud retrievals (including partly cloudy) or overcast cloudy pixels (i.e., removing partly cloudy pixels). When all valid cloudy pixels are included, including partially cloudy pixels, the N_{d_cal} distributions tend to shift toward smaller values, indicating that partially cloudy coarse pixels can influence the retrieved N_d. However, when only confident overcast cloudy pixels are retained, the mean N_{d_cal} values at the native LES and coarse satellite-like resolutions become more comparable. This behavior is observed consistently for the 1.6, 2.1, and 3.7 µm absorbing-channel retrievals, which suggests that it is not just restricted to a specific band. Therefore, satellite-like coarse-resolution retrievals can provide comparable mean N_d estimates to the native-resolution results when partly cloudy pixels are excluded, and only confidently cloudy pixels are utilized.

We show that the sensitivity s(τ), defined as the change in τ with respect to N_d at constant LWP $\left(\frac{\partial \mathit{\tau }}{\partial N_{\mathrm{d}\mathrm{\_}\mathrm{cal}}}{\mathrm{|}}_{\mathrm{LWP}}\right)$, is influenced by both retrieval errors due to 3D radiative effects and LWP categorization errors. The results reveal that 3D effects retrieval error cancel out in s(τ) calculated under 3D RT using LWP self-categorization (i.e., LWP derived from r_e and τ retrieved from 3D reflectance) to give s(τ) values that aligns closely to theoretical expectation $(\mathrm{1}/\mathrm{3})$, compared to corresponding s(τ) values obtained when non-self-consistent LWP categorization is applied.

At the native LES resolution, the average regression-based relative susceptibility (S_α,rel) computed from retrievals which utilize 3D RT reflectance is generally smaller than the corresponding 1D RT estimates at low LWP/τ for all three absorbing-channel retrievals. However, this behavior reverses toward larger LWP/τ, where the S_α,rel from 3D becomes comparable to, or larger than, the 1D RT result. This indicates that 3D RT modifies not only the magnitude of S_α,rel, but also its dependence on cloud optical state. At the coarser 800 m resolution, the differences between the 1D- and 3D-based susceptibility estimates are substantially reduced over the low-to-moderate LWP/τ regimes, although the degree of agreement remains somewhat dependent on the absorbing channel. This closer agreement at coarse resolution low- to moderate LWP regimes suggest compensation between plane-parallel retrieval biases and 3D radiative effects, yielding more comparable susceptibility estimates than those obtained at the native LES resolution.

In conclusion, we demonstrate the impact of 3D effects on bi-spectral retrieval of τ and r_e and their subsequent impact on the ability to infer N_d from such retrievals. Retrieval errors are known to have far-reaching implications, biasing aerosol–cloud interaction studies that attempt to constrain the impact of the first aerosol indirect effect (Quaas et al., 2020). Our study shows that 3D effects have significant impact on regression-based relative albedo-susceptibility at the native LES resolution than at coarser satellite-like footprints and low-to-moderate LWP and optical depth regimes. We emphasize that these susceptibility estimates are not exact 3D RT derivatives of area-averaged TOA albedo. They are regression-based diagnostics computed from locally registered domain-top albedo and are therefore sensitive to the chosen flux-registration convention. At the coarse 800 m resolution, differences between 1D- and 3D-based relative susceptibility estimates are substantially reduced across the low-to-moderate LWP/τ regimes, suggesting partial compensation between plane-parallel retrieval biases, and 3D radiative effects. Improvements in satellite retrievals of cloud properties, from which N_d is inferred, remain important – especially the development of retrieval methods that are less sensitive to 3D radiative effects. Recently, there has been increased development of polarimetric retrievals of cloud properties, offering the potential to constrain biases associated with r_e retrieval to the extent that the uppermost cloud layer microphysics captures the desired physics. However, the method of retrieving τ still utilizes the spectral technique, so any N_d retrieval derived from a polarimeter may not benefit from some of the compensating error impacts observed in this study for low to moderate LWP and τ regimes at the coarse-resolution observations. Other retrieval techniques such as cloud property retrievals using a combination of active and passive instruments can also provide a better means to retrieve N_d and should be greatly encouraged.