Spatial Aggregation of Satellite Observations Leads to an Overestimation of the Radiative Forcing Due To Aerosol‐Cloud Interactions

The estimation of cloud radiative forcing due to aerosol‐cloud interactions, RFaci (also known as the first indirect effect), relies on approximating the cloud albedo susceptibility to changes in droplet concentration, β. β depends on the cloud albedo and droplet concentration, both of which can be observed by satellites. Satellite observations are often spatially aggregated to coarser resolutions, typically 1 × 1° scenes. However, on such spatial scales, the cloud albedo tends to be heterogeneous, whereas the β approximation assumes homogeneity. Here, we demonstrate that the common practice of aggregating satellite data and neglecting cloud albedo heterogeneity results in an average overestimation of 10% in previous estimates of the RFaci. Additionally, we establish a relationship between the magnitude of the bias in β and Stratocumulus morphologies, providing a physical context for cloud heterogeneity and the associated bias. Lastly, we propose a correction method that can be applied to cloud albedo gridded data.


Introduction
The top of atmosphere albedo changes most effectively due to variations in cloud cover.The reason is that clouds are inherently brighter than the underlying surface, especially over oceans.However, clouds also vary in albedo, with some clouds showing differences in albedo that are as significant as the contrast between cloudy and cloudfree regions (I.L. McCoy et al., 2017).The albedo of clouds is a non-linear function of cloud optical thickness (τ c ) on liquid water path (LWP) and cloud droplet concentration, N d (King, 1987;Nakajima & King, 1990;Platnick & Twomey, 1994).Changes in N d were shown to affect cloud albedo through changes in droplets size, assuming that there is no change in LWP, also known as the "Twomey effect" (Twomey, 1974) or the "first indirect effect" (Lohmann & Feichter, 2001).A more recent definition is the instantaneous radiative forcing due to aerosol-cloud interactions, RF aci (Boucher et al., 2013).
The physical response of the cloud albedo to changes in N d is well described by established theory (Platnick & Twomey, 1994;Twomey, 1977).It follows that without a change in cloud water, a higher N d would lead to smaller droplets and therefore to a larger surface-to-volume ratio of the droplets.This results in a larger scattering cross section of all cloud droplets within a given cloud, and an increase in the cloud albedo (Twomey, 1977).The change in cloud albedo to a change in N d is termed cloud albedo susceptibility.Platnick and Twomey (1994) and Twomey (1991) further defined an approximation for the cloud albedo susceptibility, β: 10.1029/2023GL105282 2 of 10 with R being the observed cloud albedo.The susceptibility is expressed using the logarithmic change of N d because most cloud processes are sensitive to a relative change in N d (Carslaw et al., 2013;Kaufman & Koren, 2006).The above formulation assumes that the clouds are homogeneous in R and adiabatic, and that there is no change in the cloud water content (Twomey, 1991).The formulation also shows that the highest value of β occurs for clouds with an R of 0.5.From there, β decreases symmetrically as R deviates toward both higher and lower values.The relationship between β and R is illustrated in Figures 1d-1f.
Marine stratocumulus clouds (Sc) are low level, liquid clouds that cover extensive areas over the oceans.Due to their large spatial extent, changes in N d due to changes in aerosol concentrations can potentially exert a strong RF aci (Wood, 2012).Observations and model simulations of Sc show that their albedo is heterogenous on scales ranging from a few kilometers to hundreds of kilometers (I.L. McCoy et al., 2017;Rampal & Davies, 2020;Stevens et al., 2020;Wood & Hartmann, 2006;Zhou et al., 2021).As outlined in the previous paragraph, β strongly depends on the exact value of R of a given cloud.This means that β exhibits spatial heterogeneity on scales of kilometers to hundreds of kilometers, similar to R.
In studies that utilize satellite observations, it is a common practice to aggregate the data into coarser resolutions, typically 1 × 1° scenes (Christensen et al., 2020;Diamond et al., 2020;Gryspeerdt et al., 2017Gryspeerdt et al., , 2019;;Hasekamp et al., 2019;D. McCoy et al., 2017;D. T. McCoy & Hartmann, 2015).These scenes are represented by the mean values of the respective cloud properties.Among those is the scene mean albedo  ( R ) . Observations of  R are often retrieved from the Earth's Radiant Energy System (CERES) instrument with a footprint resolution of 20 km (Loeb et al., 2005).The 20 km scale is not able to capture smaller scale cloud heterogeneity, even more so when aggregated to 1 × 1° scenes.The process of aggregating the data into a 1 × 1° resolution smoothens out the spatial variance in R. As a result, spatial variations in cloud albedo within the scene are loss.
Previous studies have acknowledged the influence of cloud albedo heterogeneity on the parameter β (Barker, 2000;Feingold et al., 2022;Oreopoulos & Platnick, 2008).However, to the best of our knowledge, no studies have quantified the propagation of this bias in β to the bias in RF aci within the context of the common practice of spatial  , 2022), it is crucial to establish the dependence of β on cloud morphology, rather than on the generic term of cloud heterogeneity.

Methods
The analysis covers a full randomly chosen year (2014) of Moderate Resolution Imaging Spectroradiometer (MODIS) observations of Sc clouds in the south-east Pacific (SEP) ocean from the Aqua satellite.Observations were collected from the region east (west) of Longitude 70°W (115°W) and south (north) of Latitude 10°S (30°S).We chose to analyze Sc clouds because they contribute significantly to the uncertainties in cloud-radiation interactions, especially when comparing models to observations (Christensen et al., 2022;Gryspeerdt et al., 2020;Mülmenstädt et al., 2020;Neubauer et al., 2014).The SEP ocean was selected because it has the maximum annual amount of Sc clouds (Muhlbauer et al., 2014).As a proxy for the cloud albedo, we use the cloud reflectance at 0.6 μm.It is provided from MODIS level 2, Collection 6.1 cloud property data set with a spatial resolution of 1 × 1 km at nadir for cloudy pixels only (Platnick et al., 2003(Platnick et al., , 2016)).It should be noted that albedo observations in many studies are often retrieved from the Earth's Radiant Energy System (CERES) instrument with a footprint resolution of 20 km (Loeb et al., 2005).It therefore lacks the information to detect cloud heterogeneity in scales smaller than 20 × 20 km.Furthermore, the 20 × 20 km scale is often composed into 1 × 1° daily data (Wielicki et al., 1996), which is more commonly used.Our use of the 0.6 μm at a spatial resolution of 1 × 1 km, therefore allows us to detect smaller scale cloud heterogeneity.Three-dimensional radiative effects, partial pixel filling, or shadows (Cahalan et al., 1994;Cahalan & Snider, 1989) might lead to a non-physical cloud albedo heterogeneity.However, these effects are expected to be minimal in high cloud cover scenes of Sc and are more likely to occur when the sun is low in the sky, which is not the case for the subtropical regions during the Aqua or Terra overpasses.A thorough review of retrieval uncertainties can be found in D. P. Grosvenor et al. (2018), which highlighted that-while being important-uncertainties on cloud retrievals should be minimal in Sc regions.
Only data located within a satellite viewing angle of less than 45° were used to avoid biases due to the degraded spatial resolution that occurs at large viewing angles (D.P. Grosvenor et al., 2018).The MODIS data were regridded into 1 × 1° scenes, which were further filtered to include only those having less than 1% coverage of ice clouds, multi-layer clouds, and cloud top temperature <268°K.This was done to ensure that 99% of the observed clouds within a scene are liquid-containing low-level marine clouds.Only scenes with a cloud fraction greater than 0.3 were included.This was done to avoid scenes with small cloud elements, which are subject to biases in the reflectance retrievals due to partial pixel filling.Note that the low cloud fraction of such scenes means that they contribute relatively little to the RF aci .
The cloud albedo susceptibility β, is estimated from the observed 0.6 μm cloud reflectance using Equation 1.In every 1 × 1° scene, β is calculated twice.First, it is calculated directly from the 1 × 1° scene mean 0.6 μm reflectance,  R (denoted as    ).This approach is comparable to β calculated from the standard CERES 1 × 1° albedo (Wielicki et al., 1996).Second, it is calculated individually for each 1 × 1 km pixel within the scene, and then averaged over the 1 × 1° scene to obtain the scene mean value of β (denoted as   ).With respect to the constant LWP assumption that underlies the β approximation,   further assumes that the spatial LWP distribution remains the same.This is so because the 1 × 1 km pixel is treated as a single cloud entity in which LWP is assumed to remain constant, whereas the 1 × 1° scene is treated as a single unit of clouds in which the LWP can vary spatially, as long as the scene mean LWP remains the same (a spatial mean can be composed of different distributions, see, e.g., Figure 1).In this sense, we broaden the constant LWP assumption to include a constant spatial variability of the LWP.We consider   to be the true, unbiased β, which holds when neglecting potential subpixel albedo heterogeneity (Zhang et al., 2016).We then define β * , which is the relative bias in β due to cloud heterogeneity in the 1 × 1° scene: There are several measures that can be used to quantify the dispersion of a distribution.A common measure is the standard deviation.However, standard deviation does not provide information on the asymmetry of the distribution (i.e., its skewness).As can be seen in Figures 1d-1f, as well as implicitly in previous studies (Wood  & Hartmann, 2006;Zheng et al., 2018), the distribution of the cloud albedo is often skewed.The skewness, although providing information about the asymmetry of the distribution, lacks information about the range of the values that are skewed.The interpercentile range (IPR) is another measure of the dispersion of a distribution and is defined by the difference between the values of two chosen percentiles.IPR therefore accounts for the spread and asymmetry of the distribution with respect to the mean.Here, we define the IPR percentiles to be the 33rd (lowest tertile, T1) and the 66th (highest tertile, T3) of R within a given 1 × 1° scene, and further calculate the mean of  R for the upper and lower tertiles as follows: with n being the number of pixels of R within the tertile.The upper and lower tertiles include 66% of the data, while the remaining 33% are spread closer to the mean.The choice of tertiles was made to avoid the strong influence of outliers when using smaller range bins, but yet to represent the majority of the data that are distributed closely to the mean value of the whole scene.

Results
One of the underlying assumptions of the β approximation (Equation 1) is that clouds are homogeneous (Twomey, 1977).Studies have shown that even within a 1 × 1 km scale, which is the MODIS pixel resolution at the nadir, subpixel heterogeneity is present (Zhang et al., 2012(Zhang et al., , 2016)).It is even more strongly pronounced when considering a 1 × 1° scene, the typical scene resolution used to estimate RF aci from satellite observations (Christensen et al., 2020;Diamond et al., 2020;Gryspeerdt et al., 2017;Hasekamp et al., 2019;D. McCoy et al., 2017;D. T. McCoy & Hartmann, 2015).Figure 1 shows examples of 1 × 1° scenes that, despite having a similar  R = 0.35, have distinctly different visible cloud morphologies.The cloud morphology is in fact a manifestation of the cloud albedo heterogeneity, and is associated with the spatial heterogeneity of LWP (Wood & Hartmann, 2006;Zhou et al., 2021).
The dependence of β on R is demonstrated in Figures 1d and 1e.It can be seen that the β function (thin line) has a parabolic shape with a maximum at R = 0.5.Figures 1d and 1e also contain the histograms of R within the 1 × 1° scenes, from which it becomes obvious that each scene is composed of a different R distribution.Since the relationship between R and β is non-linear (thin lines in Figures 1d-1f), and because R varies within the scene (histograms in Figures 1d-1f),    cannot be equal to   .This is further explored statistically and conceptually in the following sections.

Exploring the Relationship Between Cloud Albedo Morphology and Bias in β
Figure 2a shows β * (Equation 2) as a function of  R .The color scale indicates the difference between the means of the upper and lower reflectance tertiles, ΔR 66th−33rd (see Equation 3).It can be seen that for any given  R , β * increases as ΔR 66th−33rd increases.This can be thought of as if a distribution with a given  R was to be stretched to have longer tails, while  R remains constant (see, e.g., Figure 1d with respect to Figure 1e).Consequently, the lower tail of the distribution would be linked to lower values of β (and would not be compensated by the higher values of β in the upper tail, see later Section 3.3), which would ultimately lead to a reduction in the average value of   with respect to    .When  R = 0.5, ΔR 66th−33rd is most effective in reducing   with respect to    (hence the largest β * ).This is because both tails of the distribution are associated with a lower β (see further discussion in Section 3.3).Consequently, β * is the largest when  R = 0.5 and when ΔR 66th−33rd is high, as can be seen in Figure 2a.
In Figure 2b β * is projected onto the parameter space of  R33rd and  R66th , providing an additional perspective on the influence of the dispersion of the R distribution.It can be seen that β * increases roughly perpendicular to the 1:1 line.This implies that β * depends strongly on ΔR 66th−33rd , that is, the dispersion of R.Although R is controlled by N d and LWP, the spatial variations in R are mainly driven by the variations in LWP (Wood & Hartmann, 2006;Zhou et al., 2021).In that sense, the dispersion of R can be considered as a property of the cloud morphology.) .This can be seen in Figure 1 by comparing the true-color images and the R distributions of areas b and c, in which area c has brighter thick clouds and darker thin clouds with respect to the clouds in area b.Following this argumentation, a larger ΔR 66th−33rd implies a larger geometrical difference between the deeper (high R) and the thinner (low R) clouds within a scene, which leads to a larger β * .Such inhomogeneous scenes are associated with the stage at which overcast Sc breakup to a cumulus regime (Wyant et al., 1997;Zheng et al., 2018).
The minimum β * is independent of R and occurs when ΔR 66th−33rd is close to 0 (Figure 2a).In such scenes  R ≈ R33rd ≈ R66th , which means that the scenes are homogeneous, in agreement with the assumption that is in the nature of Equation 1 (Twomey, 1977).An example of a rather homogeneous scene in Figures 1a  and 1d.An exception is when  R is within the range of 0.2 (and theoretically also for 0.8, given the symmetry of the β approximation function around  R = 0.5), where a bias close to 0 occurs also for ΔR 66th−33rd > 0 (Figure 2a).In that range the relationship between  R and β exhibits a close-to-linear behavior for relatively small ΔR 66th−33rd (see Figure 3b).The low β * , despite the positive ΔR 66th−33rd , is related to the fact that biases, introduced by both tails of the distribution of R, compensate for each other and effectively balancing out the overall effect (see Section 3.3).Negative β * can occur in scenes with low  R (Figure 2a) and low cloud cover (Figure 2d).Such scenes compose only 4.4% of the sampled scenes, which implies that the positive bias dominates (Figure 2d).

𝐴𝐴
R in both scenes, indicating that R follows a normal distribution (the blue example has a larger standard deviation than the red example).Despite the fact that the two scenes have different distribution of R,    is similar for both.This is simply because  R in both scenes is the same.  , on the other hand, would be comparably smaller in both scenes, leading to a bias in β (that would be larger for the blue example scene).The reason is the lower β values that are associated with both lower and higher R at the tails of the distribution, which reduce   with respect to    .These examples illustrate why the largest bias occurs at both  R = 0.5 and large ΔR 66th−33rd (Figure 2a), as discussed in the previous section.
Figures 3c and 3d, respectively, provide an illustration of the two examples that are shown by the red and blue arrows in Figure 3a.The figures show typical Sc closed cells with different morphologies that differ by the depth of their cores and their anvils.Closed cells are made of convective cores, in which air rises within the cores and diverges horizontally, filling the gaps between the cores (the anvils) (Goren et al., 2018;Wang & Feingold, 2009).The maximum depth of the clouds is located at their cores, and they are therefore brighter compared to the anvils, which are thinner and thus less reflective.The difference in R between the thick cores and the thin anvils 10.1029/2023GL105282 7 of 10 is ΔR 66th−33rd .Accordingly, ΔR 66th−33rd is larger in Figure 3d, which has deeper cores and thinner anvils with respect to the clouds illustrated in Figure 3c.In accordance, β * would be larger in the more heterogeneous clouds in Figure 3d.Figures 3c and 3d can be conceptually compared with the true-color images in Figures 1b and 1c, respectively.Zheng et al. (2018) have shown that the skewness of the LWP increases toward the Sc-Cu transition regions, where the marine boundary layer is deeper.The greater skewness reported by Zheng et al. (2018) corresponds to a morphology transition in which the thicker the cloudy cores are, the thinner the anvils (illustrated in Figures 3d and 3c, which complement Figure 2c), and implies increasing heterogeneity.This is aligned with the deepening-warming hypothesis that elucidates the morphology transition from Sc to cumulus clouds (Wyant et al., 1997).The increasing heterogeneity characterizing the morphology transition means that the cloud albedo is becoming prone to a larger β * .

A Conceptual View on the Bias in β Due To the Non-Linearity of the Albedo Susceptibility Function
Figure 3b illustrates a theoretical scene in which  R33rd < R < R66th < 0.5 (according to Figure 2d such scenes are the most common).In such a scene one would expect a smaller β * for a given ΔR 66th−33rd .The reason is the larger β associated with  R66th that can be thought to compensates to some extent the smaller β associated with  R33rd .However, Figure 2a contradicts this notion, indicating that in the majority of such scenes, β * for a given ΔR 66th−33rd always increases when  R deviates from 0.5 (for nearly 95% of the data β * > 0).The reason is the non-linearity of the β approximation with R (Equation 1): As  R33rd becomes smaller than  R (as in the example shown in Figure 3b), the corresponding β decreases progressively.This leads to a smaller β for  R33rd , which is larger in magnitude than the larger β associated with increasing  R66th .The same applies when  R > 0.5, but in an inverse manner.This ultimately leads to a smaller   with respect to    , thus to a positive β * .This non-linearity of β with R enhances the bias due to the dispersion effect (Section 3.2).Zhang et al. (2016) proposed a mathematical framework based on Taylor expansions (up to the second order) to quantify the impacts of subpixel reflectance variance on satellite retrievals of cloud microphysical properties.Here, we adapt this method to derive a correction for the impacts of the reflectance variability on

Correcting β for Albedo Heterogeneity
Figure 4a shows the relationship between the true β  ( and the biased β  (

𝛽𝛽 𝑅𝑅
) . It is evident that the bias generally increases as  R increases, its maximum at  R = 0.5, which is consistent with our findings in the previous sections.In Figure 4b, we present the impact of the correction method in reproducing   , yielding an R 2 value of 0.98.It should be noted that our correction may also address inhomogeneity originating from reflectance retrieval bias; however, these biases are assumed to be much smaller than the inhomogeneity caused by cloud morphology (Zhang et al., 2016).One should therefore exercise caution when using satellite observations in high latitudes, as biases may be exaggerated due to the low solar zenith angle (Grosvenor & Wood, 2014;Várnai & Marshak, 2007).
It should be noted that MODIS Level 3 daily 1 × 1° cloud product includes standard deviations for various cloud properties; however, it does not include the standard deviation for cloud reflectance.As a result, applying the aforementioned correction may not be straightforward from standard MODIS products.Attempts to use the variance of cloud optical thickness, which is correlated with cloud albedo and is provided in Level 3, did not yield an accurate correction due to retrieval uncertainties and missing pixels in the retrievals.Therefore, our 10.1029/2023GL105282 8 of 10 study emphasizes the necessity to incorporate higher-order statistical measures in Level 3 data also for the cloud reflectance.Those can be applied to improve the analysis of CERES albedo observations by accounting for cloud heterogeneity.

Conclusions
Our analysis shows that the larger the cloud albedo heterogeneity, the larger the bias in β, in accordance with previous studies (Barker, 2000;Oreopoulos & Platnick, 2008).The contribution of our study lies in quantifying the bias in β in terms of RF aci and with respect to the widespread practice of spatial data aggregation.We demonstrate that using the common grid format of 1 × 1° results in an average positive bias of 10.8% (±7.3%) in β.The bias can reach up to 50% for the most heterogeneous scenes.Considering the range of global estimates for RF aci , which vary between −0.33 and −1.1 W/m 2 (Bellouin et al., 2020;Diamond et al., 2020;Quaas et al., 2020), it can be inferred that the RF aci is overestimated by approximately 0.03-0.11W/m 2 .This approximation assumes that the bias in β exhibited by the Sc clouds in the SEP ocean represents the global bias.To provide a sense of scale, this overestimation of the RF aci is larger by 2 orders of magnitude than the positive radiative forcing attributed to contrails (Kärcher, 2018).
Our study also attributes the bias in β to cloud morphology, thus providing a more physically grounded perspective rather than the more generic cloud heterogeneity term.I. L. McCoy et al. (2022) have shown that environmentally driven transitions between Sc cloud morphologies coincide with variations in thin cloud features, which as shown here, control the magnitude of the bias in β.Future changes in the frequency of occurrence of Sc cloud morphologies could therefore influence the magnitude of RF aci .Such changes could be attributed to various factors, including cloud feedbacks (I.L. McCoy et al., 2022), changes in cloud thickness resulting from rising CO 2 concentrations (Schneider et al., 2019), and changes in aerosol levels that can influence Sc transitions (Goren et al., 2019(Goren et al., , 2022)).Our results also show that for a given cloud albedo, β tends to be larger under more heterogeneous conditions.This indicates that cloud albedo enhancement aiming at cooling the Earth by injecting aerosols to increase N d (Wood, 2021), would be most effective in homogeneous cloud conditions.Marine cloud brightening should therefore favor homogeneous cloud fields to achieve the highest level of effectiveness in brightening efforts.
Finally, a correction developed by Zhang et al. (2016) was applied to our problem.The correction method relies on subgrid reflectance variance derived from MODIS and proves to be efficient in significantly reducing the bias in β, and can be applied to scene mean albedo observations from any instrument (e.g., CERES).However, the standard MODIS gridded product does not provide subgrid reflectance variance (whereas it is provided for 10.1029/2023GL105282 9 of 10 other cloud microphysical properties).Our results therefore also emphasize the importance of including statistical properties for the cloud reflectance, in particular the standard deviation, which would allow accounting for biases arising from cloud heterogeneity.

Figure 1 .
Figure 1.Dependence of β on cloud morphology and cloud albedo distribution.(a-c) Moderate Resolution Imaging Spectroradiometer true-color images displaying 1 × 1° degree scenes containing Sc clouds with different morphologies but the same  R of 0.35.The shift from homogeneous to heterogeneous scenes (from area a to area c, respectively) can be seen by the different spatial organization of the cloud fields.Each pixel within the 1 × 1° scene covers approximately 1 × 1 km.(d-f) Histograms displaying the 0.6 μm cloud reflectance of the scenes that are shown above (a-c).The thin parabolic line is the cloud albedo susceptibility to N d , β, as a function of R (calculated using Equation 1).

Figure
Figure 2c shows that a negative relationship between  R33rd and  R66th exits for a given  R .The negative relationship is expected since a larger  R66th requires a lower  R33rd to maintain the same  R .It again emphasizes that a given  R

Figure 2 .
Figure 2. Dependence of cloud albedo susceptibility bias, β * , on cloud inhomogeneity.(a) The dependence of ΔR 66th−33rd (in color) on β * and  R .(b) The dependence of β * (in color) on  R33 and  R66 .(c) The dependence of  R (in color) on  R33 and  R66 .(d) The dependence of the number of observed scenes (in color) and cloud cover (contours) on β * and  R .

Figure 3 .
Figure 3. Illustration of the effect of cloud morphology on β.(a) β as a function of R for 2 theoretical scenes having  R = 0.5 (  R is represented by the black arrow), but different standard deviations (red and blue arrows).The arrow to the left (right) of  R represents  R33rd  ( R66th ) .(b) Same as (a) but for a scene having  R = 0.3.The black solid line is to show that for small ΔR 66th−33rd , β as function of R can be assumed to be linear.(c) A cross-section illustration of the cloud morphology corresponding to the scene indicated by the red arrow in (a).(d) A cross-section illustration of the cloud morphology corresponding to the scene indicated by the blue arrows (a) and (b).The black curved arrows indicate the flow of the air as it rises within the core of the clouds and diverges at the top of the clouds.The clouds illustrated in (d) exhibit characteristics that are more typical of mature Sc clouds along the Sc-cumulus transition.Such clouds are more inhomogeneous and therefore subject to larger β * .