Reconstructing albedo from mean cloud properties

Liquid marine clouds exert a substantial control on the Earth-atmosphere energy system through their large global coverage and high reflectivity of shortwave radiation, resulting in overall negative radiative impact. Previous studies showed that the two dominant factors determining their albedo are cloud fraction (CF) and liquid water path (LWP), but this relationship varies in regions of high aerosol loading. In this work, a simplified kernel was built to assess how well the top of atmosphere (TOA) all-sky albedo (α) can be estimated from the given properties of marine liquid clouds: CF, LWP and cloud droplet number concentration (N_d), and to what extent this approach applies globally. The study uses data retrieved from MODIS and CERES instruments for a near-global ocean domain (60° S–60° N) covering the period 2003–2021. The results showed that the albedo is only reconstructed to within 10 % in less than 40 % of cases. Several modifications of investigated method were tested for the improvement in albedo reconstructions. It was found that the number of biases decreases when the maximum solar zenith angle is considered, as well as if the CF–LWP–N_d–α kernel is calculated on a higher spatial resolution grid. The findings show that the relationship between the TOA albedo of a scene of clouds and the retrieved mean cloud properties is not universal and while accounting for regional variation is one way to address this, a better understanding of this effect is still needed to reduce uncertainty in aerosol-cloud interactions.

1 Introduction

The top of atmosphere (TOA) all-sky albedo (α, or, albedo), the fraction of incoming shortwave solar radiation reflected back into space, is a key factor in the planetary energy balance. It governs the difference between absorbed and reflected energy, with the processes that define global albedo playing an essential role in Earth's climate system (Loeb et al., 2007; Trenberth et al., 2009). The present-day estimate of the Earth's mean albedo is 0.29 (Stephens et al., 2015), which aligns with the first satellite-based measurement from the 1970s, where it was determined to be 0.30 (Vonder Haar and Suomi, 1971). Small changes in the albedo can have significant impacts on global mean temperature (Cess, 1976; North et al., 1981), both forcing climate change (Twomey, 1974) and acting as feedbacks to dampen or enhance the climate response to human activity (Budyko, 1969; Hansen et al., 1984). Despite the apparent stability in global mean albedo since the 1970s, satellite records reveal pronounced regional and temporal fluctuations, underscoring the need to understand the cloud processes that shape this balance (Loeb et al., 2024).

The global albedo can be influenced by a wide range of factors (Wielicki et al., 2005), including changes in surface characteristics (Hao et al., 2018, 2019; Miao et al., 2022; Nkemdirim, 1972; Sailor, 1995), as well as atmospheric factors, such as aerosol loading and properties (Herman and Browning, 1975), and change to the coverage and properties of clouds (Bender et al., 2011; Engström et al., 2014). With clouds being one of the most important (responsible for approximately half of the Earth's total albedo; Mueller et al., 2011) and variable (Hartmann and Short, 1980) factors, it is vital to understand how changes in cloud properties can modify planetary albedo and hence the overall energy budget.

While the albedo of a field of clouds is known to depend strongly on the cloud coverage and water path, processes that modify these, such as aerosol-cloud interactions (Bellouin et al., 2020) and cloud feedbacks (Stephens, 2005) are expected to have significant impacts on planetary albedo. This has led to several techniques to reconstruct the “scene” albedo (or changes in it) from cloud properties and their variations. Radiative kernels (Zelinka et al., 2012) link discretised cloud properties to TOA radiative fluxes and have been extensively used to calculate the strength of cloud feedbacks (Ceppi et al., 2016; Zhang et al., 2021) and the impact of aerosols (Wall et al., 2022). With their discretised and almost-linear nature, they are typically applied to monthly-mean cloud and radiative properties. The cloud radiative kernel from Zelinka et al. (2012) uses changes in the mean cloud fraction, optical depth and cloud top pressure to calculate the gridbox TOA shortwave change, although the cloud top pressure has a minimal impact. Kernel methods have been shown to accurately reconstruct cloud feedbacks in model output, comparing well to more complex methods (Zelinka et al., 2012).

Single/multi-variable regressions have also been used to characterise the relationship between cloud properties and albedo (Quaas et al., 2008), forming a critical part of observation-based calculations of the radiative forcing (Bellouin et al., 2020). Some studies have used a single, global relationship (Quaas et al., 2008), whilst others use local regressions, but often only one predictive variable at a time (creating uncertainty in the final result; e.g. Feingold et al., 2017). CF and LWP have been shown as insufficient for constraining albedo in studies with varying aerosol concentrations (Engström et al., 2015), such that aerosol-cloud-interaction studies typically constrain albedo using the mean CF, LWP and N_d. With relationships between albedo and cloud properties calculated at a climate-model gridbox scale, this method has been shown sufficient for calculating the aerosol forcing in model output (Gryspeerdt et al., 2020).

However, these kernels/relationships assume some linearity between cloud properties and albedo. While this can be partially addressed by using a larger number of bins in the radiative kernel (e.g. Gryspeerdt et al., 2019), this approach does not account for the sub-gridbox distribution of cloud properties. The impact of cloud heterogeneity on radiative transfer has been extensively studied, with cloud inhomogeneity shown to introduce biases in grid-averaged albedo calculations (Cahalan et al., 1994; Oreopoulos and Davies, 1998) and significant effects on the Earth's radiation budget through mesoscale cloud variations (Rossow et al., 2002). Climate models have developed parameterizations to account for subgrid-scale cloud variability (Barker and Räisänen, 2005), though uncertainties remain in how well mean properties capture scene albedo. For example, with a non-linear relationship between LWP and cloud albedo, the mean LWP of a field of clouds does not uniquely determine the mean cloud albedo (Zhang and Feingold, 2023), as cloud morphology and sub-grid cloud heterogeneity play key roles in controlling the relationship between albedo and CF, LWP, and/or N_d (Choudhury and Goren, 2024; Goren et al., 2023; McCoy et al., 2023).

Relationships between cloud properties and albedo are commonly used to quantify changes in cloud albedo associated with variations in cloud macrophysical and microphysical properties, for example in the context of aerosol forcing and cloud feedback studies (Wall et al., 2022; Zhang et al., 2021). It remains unclear to what extent scene albedo is accurately captured by a three-parameter decomposition (CF, LWP and N_d). While climate model results suggest this is sufficient, they do not represent the details of the sub-gridbox cloud distribution needed to accurately test this assumption.

This study aims to assess whether kernel-based decompositions of cloud albedo are justified and how well mean cloud properties (CF, LWP and N_d) can be used to reconstruct the albedo of liquid cloud scenes across the globe. Using a joint-histogram/kernel approach, this work reconstructs albedo from average cloud properties at 100 km scales, characterising regional variations in the error in the reconstructed albedo compared to observations. Different methods for accounting for these biases were assessed, providing recommendations for future observation-based calculations of aerosol forcing and cloud feedbacks.

2 Materials and methods

2.1 Data

Cloud property retrievals were obtained from the Collection 6.1 MODIS Level-3 (L3) Atmosphere Daily Global Product (MOD08_D3 and MYD08_D3 for Terra and Aqua, respectively; Platnick et al., 2017a, b). The Level-3 MODIS Atmosphere Daily Global Product is a daily global spatial aggregation of the parameters generated from the Level-2 MODIS Atmosphere Products: Aerosol (MOD04_L2, MYD04_L2), Water Vapor (MOD05_L2, MYD05_L2), Cloud (MOD06_L2, MYD06_L2), and Atmosphere Profile (MOD07_L2, MYD07_L2). The L3 statistics are summarized over a 1° equal-angle latitude-longitude grid. The MODIS L3 parameters primarily considered in this study are CF, defined here as the MODIS data field: Cloud Retrieval Fraction Liquid – Liquid Water Clouds from Successful Optical Properties Retrievals, and LWP. Ice cloud fraction (ICF) was additionally used for filtering scenes with overlying ice cloud.

Daily gridded N_d estimates from MODIS were obtained from Gryspeerdt et al. (2022a). N_d retrievals in this dataset are based on the Level-2 Collection-6 MODIS Cloud Product (MOD06_L2, MYD06_L2) and follow several sampling strategies. This work considers the strategy named G18; proposed by Grosvenor et al. (2018), it balances data quantity with accuracy, accounting for several known biases in the retrieval (Gryspeerdt et al., 2022a).

Albedo retrievals were derived from the CERES Daily Time-Interpolated TOA/Surface Fluxes, Clouds, and Aerosols (SSF1deg-Day) product (NASA/LARC/SD/ASDC, 2015b, a) and calculated as the ratio of reflected shortwave radiation to incident solar radiation:

$$\begin{array}{}\text{(1)}& {\mathit{\alpha }}_{\mathrm{CERES}}={\displaystyle \frac{F_{\mathrm{SW}}^{\mathrm{TOA}}}{F_{\mathrm{solar}}^{\mathrm{TOA}}}}\end{array}$$

where α_CERES represents the observed CERES albedo, $F_{\mathrm{SW}}^{\mathrm{TOA}}$ is the observed TOA shortwave flux, and $F_{\mathrm{solar}}^{\mathrm{TOA}}$ is the observed TOA solar insolation flux. CERES and MODIS retrievals were paired by satellite platform (Aqua–Aqua and Terra–Terra) and analysed jointly.

Finally, for investigating albedo–cloud sensitivity in stratocumulus region, the ECMWF ERA5 reanalysis (temperature at 700 hPa) was used in order to calculate the estimated inversion strength (EIS), with the formula:

$$\begin{array}{}\text{(2)}& \mathrm{EIS}=\mathrm{LTS}-{\mathrm{\Gamma }}_{\mathrm{m}}^{\mathrm{850}}\left(z_{\mathrm{700}}-\mathrm{LCL}\right)\end{array}$$

where LTS is the lower-tropospheric stability, ${\mathrm{\Gamma }}_{\mathrm{m}}^{\mathrm{850}}$ is the moist adiabat at 850 hPa, z₇₀₀ is the height of p = 700 hPa surface, and LCL is the lifting condensation level. EIS is known to be a predictor of stratus cloud amount (Wood and Bretherton, 2006), and this study assesses the relationship between EIS and differences in albedo estimates to help explain stratocumulus-to-shallow-cumulus transition.

2.2 Method

In order to reduce the overall impact of surface albedo variations, the study area was geographically limited to ocean and latitudes between 60° S and 60° N. Only cases with a maximum solar zenith angle (SZA_max) not exceeding 65° were considered, following the threshold used in the G18 sampling strategy for estimating N_d. All data were filtered through cloud fraction of ice clouds (ICF < 0.01) to minimise their contribution to the scene albedo, which restricted the analysis primarily to low-altitude clouds (with cloud top pressure usually above ∼ 680 hPa; not shown). The resulting subset of cases considered in this study is pictured in Fig. 1.

Figure 1 Number of days at each location with valid CERES albedo and MODIS cloud properties (CF, LWP and N_d) across the 2003–2021 years of data included in this study.

The main method of the study is a 3D joint histogram (kernel) constructed from global observational data, where mean albedo is calculated for discrete bins of CF, LWP, and N_d. Once constructed, this kernel functions as a reference distribution that maps combinations of these three cloud properties to expected albedo values, without requiring additional variables beyond CF, LWP, and N_d. This differs from radiative kernels such as that in Duran et al. (2025), which diagnose the TOA shortwave radiative response rather than albedo. Our approach is purely diagnostic and conceptually closer to the cloud-feedback kernels of Zelinka et al. (2012), as it is derived from observational relationships without requiring model perturbations or radiative transfer calculations.

The kernel was constructed as follows. The daily gridded 1 × 1° data were binned into discrete intervals based on CF, LWP, and N_d. Each dimension was divided into a predefined set of ranges, varying from 0 to 1 (linear scale) with 50 bins for CF, from 1 to 1000 g m⁻² (logarithmic scale) with 40 bins for LWP, and from 1 to 300 cm⁻³ (logarithmic scale) with 30 bins for N_d. For each bin, the average albedo (α_avg) was then calculated as a multi-year mean value of all 1 × 1° daily datapoints across the globe that fall into the same bin of CF, LWP, and N_d.

For each daily gridded MODIS and CERES observation, for each 1 × 1° gridbox, the difference between the average albedo in the given CF–LWP–N_d bin and the CERES albedo at the given location and time, expressed as Δα, was calculated as:

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }\mathit{\alpha }={\mathit{\alpha }}_{\mathrm{avg}}-{\mathit{\alpha }}_{\mathrm{CERES}}\end{array}$$

The relative difference in estimated albedo was calculated as:

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{rel}}=\left({\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{avg}}-{\mathit{\alpha }}_{\mathrm{CERES}}}{{\mathit{\alpha }}_{\mathrm{avg}}}}\right)\cdot \mathrm{100}\end{array}$$

3 Results

3.1 Geographic distribution of biases in reconstructed albedo

The percentage of correct albedo estimates depends on the adopted accuracy threshold, which in the present analysis is expressed both in terms of the absolute (Fig. 2a) and the relative (Fig. 2b) value of $\left|\mathrm{\Delta }\mathit{\alpha }\right|$. The results show that when the accuracy threshold is set to $\left|\mathrm{\Delta }\mathit{\alpha }\right|$ = 0.05, the reconstruction of albedo is correct in more than 80 % of the cases. However, if the stricter threshold of 0.02 is applied, this fraction decreases to 40.91 %. Although a deviation of 0.02 (corresponding to 2 percentage points) may be considered relatively small, the analysis presented in Fig. 2b demonstrates that such a change in absolute terms translates into approximately 10 % in relative terms.

Figure 2 Ratio of pixels with correctly estimated albedo as a function of the absolute difference $\left|\mathrm{\Delta }\mathit{\alpha }\right|$ (a) and the relative difference $\left|\mathrm{\Delta }{\mathit{\alpha }}_{\mathrm{rel}}\right|$ (b) between the estimated (bin-averaged) and observed CERES albedo. The red dashed lines in panel (a) are the thresholds of $\left|\mathrm{\Delta }\mathit{\alpha }\right|$ equal to 0.02, 0.05, and 0.1. In panel (b) they are the relative thresholds of 10 %, 25 % and 50 %.

The biases in albedo estimates are not distributed symmetrically around zero. On the contrary, their distribution shows a distinct left-skewed asymmetry (Fig. 3). The majority of discrepancies are concentrated around Δα ≈ 0.02, which corresponds to about 10 % relative change. At the same time, cases of underestimate exceeding 50 % relative change are clearly observed, whereas such extreme overestimates are practically absent.

Figure 3 Distribution of the absolute difference (Δα) (a) and the relative difference (Δα_rel) (b) between the estimated (bin-averaged) and observed CERES albedo. The dashed black line marks zero difference, and the red dashed lines indicate the chosen accuracy thresholds (±0.02 for Δα, ±10 % for Δα_rel).

The spatial distribution of the percentage of underestimated and overestimated albedo values, based on the threshold of 0.02, is illustrated in Fig. 4a–b. Both maps reveal well-marked zonal structures. Underestimates of Δα < −0.02 are particularly frequent at latitudes poleward of ∼ 40° in both hemispheres, whereas they occur relatively rarely in tropical regions. The opposite tendency is observed for overestimates (Δα > 0.02). A clear underestimate of albedo is visible over the regions dominated by marine stratocumulus clouds, with the most distinct example on the west coast of South America highlighted by a red rectangle in all panels of Fig. 4. Underestimates are also apparent in mid-latitudes, within regions influenced by the storm tracks of extratropical cyclones. In contrast, overestimates are most commonly observed along the Intertropical Convergence Zone (ITCZ) over the Pacific and Atlantic oceans, and to some extent over the Indian Ocean. Figure 4c presents the geographical distribution of absolute albedo bias (cases with either underestimated or overestimated albedo by more than 0.02). Generally, albedo biases are mostly evenly distributed spatially, regardless of the number of datapoints at a given location (Fig. 1). In addition to the features of a probable meteorological background, the maps show the presence of some artifacts, which may be related either to the geometry of the observations, such as the solar zenith angle or to the specific features of the retrieval algorithms.

Figure 4 Geographical distribution of the percentage of cases with underestimated albedo (a), and overestimated albedo (b), and absolute albedo bias exceeding 0.02 (c). Grey colour represents land. The red rectangle marks marine stratocumulus region used in Fig. 6a.

3.2 Drivers of biases in the reconstructed albedo

The zonal patterns shown in Fig. 4 become even more apparent when cases of under- and overestimates are separated into individual 5° latitude bands (Fig. 5). On average, the global percentage of underestimates over the analysed multi-year period amounted to 27.7 %. This value varies significantly with latitude: from less than 10 % around the equator to more than 50 % at latitudes higher than 50° in both hemispheres. Similarly, overestimates occurred in 32.1 % of cases globally; although the latitudinal variability was slightly less pronounced than in the case of underestimates, the values still ranged from 15 %–20 % at high latitudes to around 50 % in the equatorial zone.

Figure 5 Distribution of (a) under- and (b) overestimated cases of albedo across 5° latitude bands. The boxes show the 25th–75th percentile range, with the line inside marking the median. The whiskers extend to the minimum and maximum values. Blue dots indicate the mean, and black dots mark outliers.

The biases observed over marine stratocumulus regions were examined in greater detail in the context of the transition from thick stratocumulus decks close to the South American coastline towards shallow cumulus clouds occurring further westward over the ocean (the region marked in red in Fig. 4). To investigate this, the relationship between EIS and Δα was analysed for the test period 2007–2019 using Aqua satellite data. Figure 6a presents the joint EIS–Δα histogram, normalised by each EIS column. The results show a clear dispersion of Δα values along a characteristic curve: negative values dominate for EIS below 0 K, positive values are typical for EIS between approximately 0 and 15 K, and negative values appear again for EIS above 15 K. EIS values exceeding 15 K most likely then correspond to cases of thick stratocumulus, which in Fig. 4a are observed by CERES to be brighter than calculated from the CF–LWP–N_d values using the kernel approach; clouds with lower EIS are instead associated with shallow cumulus, which often had an overestimated albedo (Fig. 4b). EIS values below 0 K might represent either situations with very limited cloudiness or cases of convective clouds, that is, conditions without a well-defined temperature inversion.

Figure 6 (a) The joint EIS–Δα histogram for the 10–30° S, 85–65° W region (Aqua only, 2007–2019), marked with a red rectangle on Fig. 4. (b) The joint SZA_max–Δα histogram for the global ocean (2003–2021). In both plots each column is normalised so that it sums to 1, showing conditional probabilities P(Δα | EIS) (a) and P(Δα | SZA_max) (b). Red dashed lines indicate Δα = −0.02 and Δα = 0.02, black dashed line – Δα = 0.00.

Figure 6b shows the joint histogram of SZA_max and Δα, normalised by the explanatory variable in each case. The dependence between these two quantities is particularly distinct and appears to explain, to a large extent, the distribution of albedo biases presented in Fig. 3. Numerous cases of overestimates around Δα ≈ 0.02 are typical for SZA_max below around 35°, whereas above ∼ 40° a clear relationship between SZA_max and Δα exists, with Δα decreasing significantly as SZA_max increases. For solar zenith angles above 50°, differences in albedo estimates reach values as high as 0.05–0.07. This explains the significant number of strong underestimates also visible in Fig. 3, and might also account for the underestimates observed at high latitudes in Fig. 4a.

3.3 Improving albedo estimates

There is a clear relationship with the geometry of observations, as well as expected differences resulting from the variable properties of clouds. These deviations are not fully captured by the mean cloud properties. In this part of the study several possible methodological adjustments were tested with an aim of improving the accuracy of albedo estimates (Table 1).

Table 1 Ratio of underestimated (Δα < −0.02) and overestimated (Δα > 0.02) cases of albedo for the considered methodological modifications. Values in brackets indicate the difference with respect to the original method.

As a first step, modifications related to cloud fraction were considered. Since it can be expected that for low cloud fraction the varying brightness of the underlying ocean surface (resulting e.g. from the presence of atmospheric aerosol, waves, or the angle of solar rays) may have a larger influence on the albedo, only cases with CF greater than 0.95 were selected for analysis (methodological modification no. I). Methodological modification no. II involved performing an estimation using a much larger number of bins (1000) for CF, in order to verify that a smaller number of bins (50) was sufficient to capture the characteristic U-shaped distribution of CF (with very small or nearly complete cloud cover occurring most frequently, while intermediate values appear relatively rarely).

Two methodological adjustments related to the SZA_max were then tested. In the first case, all pixels with SZA_max ≥ 40° were excluded in order to eliminate situations where the relationship between SZA_max and Δα becomes approximately linear (modification no. III). The second correction consisted of determining the mean Δα value within each 1° SZA_max bin, and subsequently applying this correction to every estimated albedo case (modification no. IV).

The next group of modifications was aimed at testing possible corrections related to specific local conditions not captured by a global-mean analysis. These included: calculating α_avg separately for each 5° latitude band (modification no. V), calculating α_avg separately for each 5° latitude–longitude grid cell (modification no. VI), and calculating α_avg separately for each 1° latitude–longitude grid cell (modification no. VIIa). Additionally, a combined adjustment of both filtering out the SZA_max ≥ 40° pixels and computing α_avg separately for 1° grid cell (modification no. VIII) was tested.

In all of the above modifications, the kernel was constructed using CF, N_d, and LWP values from the full annual time series. The last two modifications were related to the temporal aggregation used in the kernel construction. In modification IX, the kernel is instead constructed from monthly-averaged CF, N_d, and LWP, while in modification X it is constructed from daily-averaged values of CF, N_d, and LWP. When evaluating the difference between estimated and observed albedo, the reconstructed albedo is compared against observations corresponding to the same month (modification IX) or the same day (modification X), respectively. Table 1 shows the ratio of under- and overestimated cases of albedo across the globe for each tested methodological modification.

Computing α_avg within 5° latitude bands yielded only a modest reduction in errors, and the total share of incorrect estimates decreased only from 59.1 % to 50.0 %. At the 1° grid (pixel) level, the error rate dropped nearly sixfold, to 11.2 %. Since the kernel was built for the ocean, the underlying surface albedo is most likely not affecting this improvement, and the other sources of variability that were not considered by the CF–LWP–N_d–α function seem to impact the estimates – possibly other cloud properties and regimes, such as cloud morphology, that are more location-dependent. Both modifications accounting for the SZA_max – the exclusion from the analysis of pixels with SZA_max ≥ 40° (modification III), as well as implementing a correction of mean Δα within each SZA_max interval (modification IV) – resulted in the improvement that was greater for overestimates than for underestimates. Figure 7a–b shows the histogram of Δα after applying modification no. IV. Many of the corrected biases correspond to overestimates previously identified in the ITCZ region (not shown).

Figure 7 Distribution of the absolute difference (Δα) and the relative difference (Δα_rel) between the estimated (bin-averaged) and observed CERES albedo for modifications no. IV (a–b) and no. VIII (c–d).

Restricting the analysis to cases with CF > 0.95 did not improve the estimates; instead, it substantially increased the number of overestimates. Increasing the number of CF bins had a very limited effect on the results – the improvement was negligible, showing that 50 bins is already sufficient to represent the non-linearity in the CF-albedo relationship.

The most effective improvements in the estimates were achieved when α_avg was calculated separately for individual grid cells (modification VIIa). However, the results of this modification were significantly affected by bin sparsity, which in many cases have only been filled once. As a result, the reconstructed albedo is effectively drawn from the same datapoint, that is later used in the comparison, meaning that the resulting bias estimates do not represent an independent test of a time-averaged kernel. In order to reduce strong undersampling of the kernel, for modification VIIa two alternative bin configurations were examined. While modification VIIa retained the original 50 × 40 × 30 bin configuration, modifications VIIb and VIIc used coarser bin structures (10 × 10 × 10 and 5x5x5, respectively) to improve bin population. As shown in Table 1, these coarser bin configurations reduced the sparsity problem but yielded only marginally better accuracy than the global kernel method, suggesting that simply widening bins does not fully resolve the trade-off between spatial resolution and kernel robustness.

Modification VIII – a combination of the two methods, IV and VIIa – produced the most accurate results: the overall error percentage decreased to 5.4 %, i.e. more than tenfold compared to the original value; moreover, the share of under- and overestimates became nearly equal (Fig. 7c–d).

When applying modifications IX and X, an improvement in the method is observed, with the ratio of correct albedo reconstructions increasing from 40.9 % for no temporal aggregation (original method without modifications) to 57.7 % when the kernel is constructed from the monthly-averaged CF, N_d, and LWP and to 60.6 % when it is constructed from the daily-averaged CF, N_d, and LWP (Table 1). However, it should be noted that the daily-averaged version (modification X) may also be subject to bin undersampling, as the number of datapoints available on individual days is considerably smaller than over longer averaging periods. Therefore, the daily results should be interpreted with similar caution as the finest spatial resolution case. Generally, the tests demonstrated that among the investigated modifications, only those accounting for regional and seasonal variability of albedo provided a substantial improvement in the estimates.

4 Discussion

These results show that the reconstructed albedo of a scene of clouds, estimated using the mean albedo for each CF–LWP–N_d bin, exhibits systematic biases, with underestimates prevailing at higher latitudes, in known stratocumulus regions and storm tracks, while overestimates are most pronounced in the tropics and the ITCZ. These contrasting biases may reflect differences in cloud morphology and sub-grid heterogeneity: in stratocumulus regions, nonlinear radiative effects and unresolved horizontal variability can lead to an underestimation of scene albedo when using mean cloud properties, whereas in the tropics, highly heterogeneous convective cloud fields may cause mean properties to overrepresent optically thick elements, resulting in albedo overestimates. While these biases are linked to other properties, such as solar zenith angle (Fig. 6b), mean cloud properties are only able to reconstruct the albedo of a scene of clouds accurately when regions are treated separately (Table 1).

As discussed in the results section, the finest spatial resolution modification (1° × 1°) suffers from significant bin undersampling, with many bins filled only once over the study period, limiting its validity as an independent test of the kernel. Among the tested alternatives, the 5° × 5° modification possibly represents the most reliable high-resolution approach, maintaining adequate bin population while capturing meaningful regional variation in the CF–LWP–N_d–albedo relationship. Future developments might address undersampling through alternative methodologies (e.g., machine learning approaches), potentially enabling robust use of finer spatial resolution data.

Regime-dependent retrieval biases could play a role in these biases in the reconstructed albedo. The N_d retrieval has significant biases in broken cloud fields (Grosvenor et al., 2018), which could contribute to the relationship between Δα and EIS. However, as this N_d bias exists even in regionally-specific studies, the accuracy of the local reconstructions (Table 1) suggests that retrieval biases are unlikely the main factor.

This study focuses on the ocean to reduce the impact of surface albedo variations, but other sources of variation in the clear-sky albedo (such as variations in aerosol loading) might also lead to biases in Δα. The patterns of biases in the Δα field (Fig. 4) are not consistent with the patterns of high aerosol optical depth and the bias remains even when considering only high cloud fraction cases (minimising the impact of surface/clear-sky albedo errors), suggesting that aerosol variability is unlikely to explain these results.

Even limiting this study to the global oceans (restricting surface albedo variability) and accounting for solar zenith angle variations (which corrects for latitude and time of year; Wall et al., 2022), only half of reconstructed albedo values are within 10 % of the observed value. This suggests that further properties of the cloud field may be necessary to accurately reconstruct the scene albedo using a single global relationship. The relationship to EIS (Fig. 6a) hints at an importance of cloud regime. The sub-gridbox LWP and cloud optical depth distributions are very different for shallow cumulus (at low EIS) and stratocumulus (at high EIS) (Bretherton et al., 2019). Given the non-linear relationship between LWP and cloud albedo (Platnick and Twomey, 1994), these regional variations could create local biases in the reconstructed albedo. If this additional parameter changes in response to aerosol loading, such as during the transition between open and closed celled stratocumulus (Rosenfeld et al., 2006), this could lead to an additional component of the aerosol forcing not covered by the usual decomposition into Twomey (N_d), LWP and CF adjustments common in observation-based studies (Bellouin et al., 2020). As this effect is not represented in climate models, the closure achieved in previous studies (e.g. Gryspeerdt et al., 2020) could be misleading.

The high accuracy/low Δα achieved when using locally-specific kernels identifies a near-term pathway to address this issue, demonstrating that with the correct explanatory variables, an accurate albedo reconstruction is possible. However, it also highlights the need to identify these additional factors that control the albedo of a cloud scene, beyond the gridbox mean cloud properties.

5 Conclusions

The two dominant factors determining liquid marine clouds albedo (α) are cloud fraction (CF) and liquid water path (LWP), with aerosol loading further impacting albedo's variability, giving additional benefits from using N_d as a further explanatory variable (Bellouin et al., 2020; Engström et al., 2015).

This study uses near-global (60° S–60° N) MODIS and CERES retrievals for marine regions from Terra and Aqua satellites spanning over the years 2003–2021 to build a simplified method to calculate albedo based only on CF, LWP and N_d values and examine the spatial differences in the sensitivity of albedo to cloud properties.

It was demonstrated that the percentage of datapoints in which the reconstructed albedo biases (relative to the measured CERES albedo) were > ±0.02 (a relative bias around ±10 %) can be as high as ∼ 60 % (Fig. 3). Moreover, it was shown that the CF–LWP–N_d–α relationship varies significantly across the globe. Underestimates are particularly frequent in regions of known stratocumulus regimes – on the west coast of South America, California and South Africa – as well as along the midlatitude storm tracks. In contrast, overestimates occur mainly in the tropics. Clear zonal patterns in the accuracy of albedo estimates suggested a potential relationship to the solar zenith angle of satellite observations. The average Δα for SZA_max up to 30–35° was found to be largely constant at a 0.02 overestimate. For the SZA_max > ∼ 40° however, Δα decreases significantly with the increase of SZA_max, reaching about −0.05 for SZA_max higher than 60°.

Several modifications of the initial method were tested in an attempt to improve the albedo estimates. Generally, correcting for SZA_max significantly decreases the number of overestimates with a modest decrease in underestimates. The highest apparent number of correct estimates (∼ 94.6 % for $\left|\mathrm{\Delta }\mathit{\alpha }\right|$ ≤ 0.02) occurs when the average albedo for each CF, LWP and N_d bin is calculated at a 1° grid resolution. However, at this resolution the CF–LWP–N_d bins are often filled with only a single datapoint, and the resulting accuracy therefore does not represent an independent test of a time-averaged kernel. At the finest (1° × 1°) spatial resolution, the reconstructed albedo biases must therefore be interpreted with caution, as the kernel is frequently undersampled and the apparent reduction in bias partly reflects self-matching of individual datapoints. Among the tested spatial resolutions, the 5° × 5° modification provides a more reliable high-resolution estimate, maintaining adequate bin population while still substantially outperforming the coarser modifications, suggesting that if bin undersampling could be resolved, the 1° × 1° results would indeed represent the true optimal performance.

Although it was shown that there are some geographical patterns in the CF–LWP–N_d–α relationship and possible modifications were suggested to improve the albedo estimates, significant uncertainties remain. Restricting the analysis to more limited regions, and thus accounting for local conditions, leads to a reduction in the reconstructed biases, suggesting that there are explanatory variables that could be used, beyond mean CF, N_d and LWP. However, the factors driving this regional variation are not yet clear, and a wider study further addressing this will be essential for a more complete understanding.