Calculating the Climatology and Anomalies of Surface Cloud Radiative Effect Using Cloud Property Histograms and Cloud Radiative Kernels

Clouds cover over two thirds of the Earth’s surface, and the radiative effects of clouds are important for the Earth’s global and regional energy budget (Trenberth et al., 2009). The clouds cool the Earth’s surface by reflecting shortwave (SW) solar radiation back to space and heat the surface by emitting downward longwave (LW) thermal radiation (Intergovernmental Panel on Climate Change, 2014). Meanwhile, the clouds heat the atmosphere by absorbing both SW radiation and LW radiation and cool the atmosphere by emitting LW radiation. The summation of cloud radiative effects (CRE) on the Earth’s surface and atmosphere is equal to the CRE on the top of atmosphere (TOA) fluxes, which can be directly measured by satellites.

The net CRE on surface, atmosphere and TOA depends on cloud properties, solar zenith angle, surface emissivity, and many other meteorological variables. When the global surface temperature changes in response to climate forcing, the cloud properties all over the world change subsequently, and the corresponding changes in cloud radiative effect would in turn change the global surface temperature, resulting in feedback to the Earth’s climate system (Intergovernmental Panel on Climate Change, 2014). Uncertainties about the magnitude of cloud feedback have been identified as the primary contributor to the uncertainty of equilibrium climate sensitivity (Webb et al., 2006; Andrews et al., 2012; Caldwell et al., 2016; Zelinka et al., 2020).

TOA fluxes are frequently used to analyze the effect of cloud feedback on global mean surface temperature changes, because long-term changes in global surface temperature are primarily driven by changes in total heat content of the climate system, and the latter is determined by net fluxes at TOA. The magnitude of cloud-induced radiative anomalies at TOA (ΔR_cloud,TOA) can be calculated with the partial radiative perturbation (PRP) method (Colman and McAvaney, 1997), or be calculated with a combination of non-cloud radiative kernels (NCRK) and CRE (Shell et al., 2008; Soden et al., 2008), but these methods do not directly link ΔR_cloud,TOA to the changes in cloud properties. Zelinka et al. (2012a) calculated a set of cloud radiative kernels (CRK) using the Fu-Liou radiative transfer model,

where τ is the cloud optical thickness, CTP is the cloud top pressure, CRK_TOA is the TOA CRK for a specific CTP-τ bin in a specific latitude with the surface albedo of α, C is the area fraction of the cloud in the CTP-τ bin, R_TOA is the net radiative fluxes at TOA, R_TOA,overcast is the net radiative fluxes at TOA when a specific type of cloud covers 100% of the sky, and R_TOA,clr denotes clear-sky radiative flux. Then the CRKs can be used to quantify ΔR_cloud,TOA together with cloud property histograms produced by climate models at each latitude, longitude and time step:

where i and j are bin numbers of τ and CTP, respectively. Using the TOA CRKs, Zelinka et al. (2012b, 2016, 2020) attributed the changes in ΔR_cloud,TOA under global warming to changes in cloud amount, cloud optical depth and CTP. Zhou et al. (2013, 2014) calculated TOA CRKs with meteorological variables from reanalyses, and explained how cloud properties and ΔR_cloud,TOA change in response to interannual global surface temperature anomalies. Yue et al. (2016) calculated TOA CRKs based on satellite observations, and their observation-based CRKs are generally consistent with CRKs calculated with radiative transfer models. Huang et al. (2019) analyzed the band-by-band contributions to the longwave cloud radiative feedback using spectrum-based CRKs. Wang et al. (2020) built a set of CRKs using the Beijing Climate Center radiative transfer model, and analyzed regional cloud radiation anomalies at TOA.

Although TOA fluxes have been successfully applied to analyze the relationship between global climate feedbacks and global mean surface temperature changes, it is inconvenient to diagnose changes in regional surface temperature using TOA fluxes, because the total heat content of regional climate system is not only determined by TOA radiative fluxes, but also significantly affected by horizontal heat advection. Instead, surface radiative fluxes are preferred to diagnose regional surface temperature changes (Colman, 2015). There are existing satellite-estimated surface radiation products that can be used to link the observed surface CRE changes to changes in cloud properties (Zhang et al., 2006; Gui et al., 2010; Young et al., 2018; Yang and Cheng, 2020), but these products are not convenient for attribution analyses, and could not be applied to analyses associated with climate models. Considering that the TOA CRKs have been widely used to diagnose global cloud feedback processes, it would be useful to develop surface CRKs to analyze regional cloud feedback processes.

In this paper, we developed a set of surface CRKs as a function of cloud base pressure (CBP) and τ with a radiative transfer model, and produced a set of surface, TOA and atmospheric CRKs on CTP-τ histograms. We calculated climatological CRE and cloud-induced radiative anomalies at surface using climate models and Moderate Resolution Imaging Spectroradiometer (MODIS) observations and compared them with values computed from surface fluxes.

In this study, the Rapid Radiative Transfer Model (RRTM, Mlawer et al., 1997) was used to build the CRKs. RRTM utilizes the Correlated-K approach to calculate fluxes and heating rates efficiently and accurately. We selected the 8-stream discrete ordinate radiative transfer solver (DISORT solver in RRTM) to perform SW scattering calculations and used the RRTM radiative transfer solver to perform LW calculations. We calculated the climatological surface temperature, surface albedo, air temperature, ozone concentration and humidity at each month of year, latitude, longitude and levels using 2000−18 data from ERA-interim (Dee et al., 2011) reanalysis, and inserted these values into the RRTM input files. The CO₂ concentration was set to be fixed at 400 ppm (close to present-day level), and the concentration of other trace gases was set to be same as that in an RRTM example file representing midlatitude summer. The background concentration of these trace gases has negligible impact on the CRKs.

Before the CRK calculations, we sorted the clouds into categories using criteria that are closely related to CRE. According to Figs. 1(a-b), the LW CRE at TOA, defined as the difference between all-sky TOA fluxes and clear-sky TOA fluxes, depends primarily on CTP, and is not sensitive to CBP; the SW CRE at TOA is primarily sensitive to τ. Therefore, the TOA CRKs are usually built on CTP-τ cloud fraction histograms. However, the LW CRE at the surface is primarily determined by CBP and is not sensitive to CTP (Figs. 1c-d), so we built the surface CRKs on CBP-τ cloud fraction histograms. Following Zelinka et al. (2012a), we divided the pressure level into 7 bins (50−180 hPa, 180−310 hPa, 310−440 hPa, 440−560 hPa, 560−680 hPa, 680−800 hPa, and 800−1000 hPa), and divided the optical depth into 7 bins (0−0.3, 0.3−1.3, 1.3−3.6, 3.6−9.4, 9.4−23, 23−60, 60−380). We calculated the CRE of ice and liquid clouds separately, so there are altogether 7 × 7 × 2 types of clouds for each latitude, longitude and month of the year.

Figure 1.  Relationships between TOA/surface LW CRE and CTP/CBP in a group of idealized experiments. In each idealized experiment, a homogenous cloud layer with the optical depth of 40 was inserted into an atmospheric column, where the meteorological parameters were from an RRTM example file representing midlatitude summer. The CTP and CBP in each sub-experiment were set to have different values, and the corresponding CRE was simulated with RRTM. Each point in the figure denotes value from an individual RRTM sub-experiment.

We calculated the clear-sky radiative fluxes at the surface and TOA at each grid box for each month of the year. The computation cost is quite high to calculate the SW radiation at each hour of the year and at each location of the earth, so we calculated the SW TOA and surface clear-sky radiative fluxes for each latitude, month of the year, surface albedo bin and solar zenith angle bin using RRTM, and calculated the solar zenith angle at each latitude, month of the year and hour of the day. Subsequently, we were able to calculate the mean SW radiative fluxes at each latitude, longitude and month of the year. LW radiation is independent of solar zenith angle, so the LW TOA and surface clear-sky radiative fluxes were directly calculated for each latitude, longitude and month of the year.

Then we calculated the radiative fluxes for overcast-sky conditions. Corresponding to each clear-sky case, we repeated 98 overcast-sky RRTM calculations to calculate the average CRE of each type of cloud. In each overcast-sky RRTM calculation, a single-layer cloud was inserted into the atmospheric column by setting liquid or ice water content to nonzero values at the level closest to the specified cloud level following Zelinka et al. (2012a). The effective cloud particle radius for liquid and ice clouds were set to be 10 and 30μm, respectively, following Zelinka et al. (2012a). For each cloud property bin, the log-scale central value of τ and central value of CBP were used in the calculations. Then we calculated the average radiative fluxes for each latitude, longitude, month of the year, τ and CBP bins. With these values, the surface CRK of each type of cloud for each month of the year was calculated as:

where $ {R}_{\mathrm{S}\mathrm{F}\mathrm{C}} $ is the surface radiative flux, $ {R}_{\mathrm{S}\mathrm{F}\mathrm{C},\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{c}\mathrm{a}\mathrm{s}\mathrm{t}} $ is the surface radiative flux for cloudy sky, and $ {R}_{\mathrm{S}\mathrm{F}\mathrm{C},\mathrm{c}\mathrm{l}\mathrm{r}} $ is the clear-sky surface radiative flux.

The global average surface CRKs for ice clouds and liquid clouds are shown in Figs. 2(a-c) and Figs. 2(d-f), respectively, and their differences are shown in Figs. 2(g-i). In general, thin-low clouds create a mildly positive net CRE on surface fluxes, while thick-high clouds have a strongly negative net CRE. All clouds cool the earth’s surface by blocking incoming solar radiation, so SW CRK at surface is negative for all types of clouds. The SW CRE is not sensitive to CBP, but is negatively correlated with τ, because thicker clouds scatter more solar radiation back to space. All clouds warm the earth’s surface by emitting downward LW radiation, so the global mean LW CRK at surface is positive for all types of clouds. Low clouds are usually warmer than high clouds, so low clouds emit more downward LW radiation to the surface, and the corresponding surface CRKs are more positive than those of high clouds. Particles in liquid water clouds are spherical, while the crystals in ice clouds are non-spherical, and the asymmetric factor of spherical liquid clouds is greater than that of ice clouds, so the SW CRK of liquid clouds is slightly less negative than ice clouds for a specific CBP-τ bin. The effective particle radius of ice clouds is greater than that of liquid clouds, and the overall LW emissivity of ice cloud particles is also greater than liquid clouds, so the LW CRK of liquid clouds is slightly less positive than ice clouds for a specific CBP-τ bin. These liquid and ice CRKs might be useful for study of the cloud phase shift feedback with separate ice and liquid histograms in the future. The cloud radiative effects of ice and liquid clouds can be calculated using the corresponding CRKs for ice and liquid clouds. For mixed-phase clouds, the cloud radiative effect can be calculated using the mean value of ice CRK and liquid CRK weighted by the optical thickness contributed from ice cloud particles and liquids cloud particles, respectively.

Figure 2.  Global and annual mean surface CRKs on CBP-τ cloud fraction histograms. (a−c) Surface CRKs for liquid clouds only. (d−f) Surface CRKs for ice clouds only. (g−i) Difference between the surface CRKs of liquid clouds and ice clouds. (j−l) Surface CRKs with combined liquid and ice clouds. The left, middle and right columns are for net CRKs, SW CRKs and LW CRKs, respectively.

Next, we combine the liquid and ice CRKs based on the cloud temperature of each latitude, longitude, month and altitude. To account for mixed-phase clouds, we calculated the CRK for each CBP, τ, latitude, longitude, and month bin as a weighted mean value of liquid and ice CRKs,

where p_liquid is the area proportion of liquid water cloud particles among all cloud particles in the corresponding CBP-τ bin. Ideally, the value of p_liquid for a thick cloud layer might be calculated for individual altitude levels, but the vertical distribution of cloud phase is highly uncertain for thick clouds, and the assumptions on cloud phase would result in biases for both cloud retrieval processes and CRK calculation processes. To reduce the uncertainty induced by cloud phase assumptions on estimated radiative effects, it is preferred that the cloud phase assumptions in kernel calculations are the same as that in the cloud retrieval processes, then the bias induced by the assumption in the inversion processes (satellite retrieval) would be compensated by a bias with the opposite sign in the forward processes (CRK calculations). The cloud phase algorithms for passive satellites usually retrieve the cloud top phase for thick clouds (for example, Baum et al., 2012), so the retrieval processes of τ in cloud property histograms usually contain the assumption that the cloud phase is the same as the cloud top. Therefore, we calculated the value of p_liquid using cloud top temperature and the statistical relationship between cloud temperature and liquid cloud probability in Hu et al. (2010). To estimate the cloud top temperature of each CBP-τ bin, we collocated the observations of CloudSat (Stephens et al., 2002) cloud profiling radar, which can penetrate cloud layers and measure the cloud base accurately, and Aqua MODIS level-2 cloud property products (Platnick et al., 2017) during 2010−12. The combined surface CRKs are shown in Figs. 2(j-l). The SW CRK of liquid clouds in a specific CBP-τ bin is less negative than that of ice clouds, so the combined SW CRK is slightly more negative for high clouds than that of low clouds (Fig. 2k). Other features of the combined surface CRKs are very similar to that of liquid and ice surface CRKs (Fig. 2j-l).

Figure 3 shows the net surface CRKs for different latitudes and different seasons. In tropical regions, the solar irradiance does not vary significantly over the seasonal cycle. In polar regions, the net CRK is positive for all types of clouds in winter, and is negative for most clouds in summer; also, the annual net CRKs in polar regions are more positive than those in the tropics.

Figure 3.  Surface net CRKs on CBP-τ cloud fraction histograms for different latitudes and seasons. (a−c) Annual mean surface CRKs. (d−f) Surface CRKs in June, July, and August. (g−i) Surface CRKs in December, January, and February. The left, middle and right columns are for net CRKs in the tropics, midlatitudes and polar regions, respectively.

Figure 4 shows the spatial distribution of CRKs for three illustrative CBP-τ bins, which represent cirrus, altostratus and stratocumulus, respectively. For a specific latitude, the SW CRKs over land are usually less negative than that over oceans due to differences in surface albedo. LW CRKs over plateaus are generally greater than those over plains because less downward cloud thermal radiation is absorbed by the atmosphere. The temperature also affects SW CRK by determining the cloud phase.

Figure 4.  The spatial distribution of CRKs for three illustrative CBP-τ bins. (a−c) represent cirrus bin (180 < CTP < 310 hPa, 1.3 < τ < 3.6), (d−f) are for altostratus bin (560 < CTP < 680 hPa, 3.6 < τ < 9.4), and (g−i) are for stratocumulus bin (800 < CTP < 1000 hPa, 9.4 < τ < 23). Missing values are marked with white color.

The CRKs on CBP-τ cloud fraction histograms could be used to calculate the climatology and anomalies of surface CRE if there are cloud property products on CBP-τ histograms. However, while there exist many CTP-τ cloud fraction histograms derived from observations and climate models, it is difficult to find any CBP-τ cloud fraction histograms at present. It is plausible to convert CTP-τ histograms to CBP-τ histograms using advanced statistical methods (e. g., deep learning algorithms), but it is more convenient to calculate CRE with CRKs on CTP-τ histograms, so we converted the surface CRKs from CBP-τ fields to CTP-τ fields using the statistical relationship between CTP, CBP and τ.

The surface CRE can be calculated with CRKs on either CBP-τ or CTP-τ cloud fraction histograms:

where C₁ and C₂ denote the cloud fraction of CBP-τ and CTP-τ histograms, respectively. We denote the probability that the cloud base pressure of a cloud layer with CTP=CTP_i and τ=τ_j equals to CBP_k is $ P\left({\tau }_{i},\mathrm{C}\mathrm{T}{\mathrm{P}}_{j},\mathrm{C}\mathrm{B}{\mathrm{P}}_{k}\right) $, then we have

Combining Eqs. (5−7) we have

We used the collocated CloudSat-MODIS data to calculate the probability function $ P\left({\tau }_{i},\mathrm{C}\mathrm{T}{\mathrm{P}}_{j},\mathrm{C}\mathrm{B}{\mathrm{P}}_{k}\right) $ at each longitude and latitude. The CloudSat and Aqua satellites passed by each location within several minutes, and we collocated the nearest MODIS pixel for each CloudSat pixel. Then CBP calculated from CloudSat cloud profile, CTP and τ obtained from MODIS products are used to build CTP-CBP-τ joint histograms. The probability function $ P\left({\tau }_{i},\mathrm{C}\mathrm{T}{\mathrm{P}}_{j},\mathrm{C}\mathrm{B}{\mathrm{P}}_{k}\right) $ is subsequently calculated from the CTP-CBP-τ joint histograms. Clouds of some CTP-τ bins were not observed in some regions, and the probability functions for these clouds and regions are missing. These missing values were filled up with $ P\left({\tau }_{i},\mathrm{C}\mathrm{T}{\mathrm{P}}_{j},\mathrm{C}\mathrm{B}{\mathrm{P}}_{k}\right) $ values averaged over regions with similar cloud features (classified by low cloud fraction with CTP > 680 hPa). If regional mean $ P\left({\tau }_{i},\mathrm{C}\mathrm{T}{\mathrm{P}}_{j},\mathrm{C}\mathrm{B}{\mathrm{P}}_{k}\right) $ was also missing (e. g., thick-low clouds over deserts), then the global average value was used to fill up the regional values.

Figures 5(a-c) shows $ {\mathrm{C}\mathrm{R}\mathrm{K}}_{\mathrm{S}\mathrm{F}\mathrm{C},\mathrm{C}\mathrm{T}\mathrm{P}} $, the surface CRK on CTP-τ histograms calculated from Eq. (8). The LW CRKs for thick-high clouds (cloud height determined by cloud top) are close to that of low-thick clouds, which is consistent with the fact that most thick-high clouds are deep convective clouds with low cloud base levels.

Figure 5.  Global and annual mean surface, TOA and atmospheric CRKs on CTP-τ cloud fraction histograms.

We also calculated a set of TOA CRKs using the clear-sky and overcast-sky TOA fluxes in the RRTM simulations with Eq. (1). In general, our TOA CRKs (Figs. 5d-f) are similar to the TOA CRKs in Zelinka et al. (2012a), but their LW CRKs are the same for all locations at a specific latitude, whereas the LW CRKs in our study are a function of both longitude and latitude. Hence, the TOA CRKs in this study contain more information than that of Zelinka et al. (2012a). Comparisons of TOA and surface CRKs indicate that clouds have similar effects on TOA and surface SW radiative fluxes, but the LW radiative effects of clouds on TOA and surface fluxes are very different. The temperature of high clouds is low, and the LW emission of these clouds is weak, so high clouds have more positive CRE on TOA fluxes and less positive CRE on surface fluxes. On the contrary, the temperature of low clouds is high, and the LW emission of the low clouds is strong, so they have less positive CRE on TOA fluxes and more positive CRE on surface fluxes.

The difference between the CRE at the TOA and the surface equals the CRE for the atmosphere, so we calculated the atmospheric CRKs by subtracting surface CRKs from TOA CRKs (Figs. 5g-i). The absolute values of SW atmospheric CRKs are small (Fig. 5h). A portion of SW radiation is absorbed by clouds, so the SW atmospheric CRK is positive for most clouds except for thin-high cirrus clouds. These cirrus clouds reduce the amount of SW radiation that is absorbed by water vapor in the lower troposphere, and since they absorb little SW radiation, their SW radiative effect has a slight cooling effect on the atmosphere. Cold high clouds absorb more LW radiation than their emitted thermal radiation, so their radiative effects warm the atmosphere. Warm low clouds emit more thermal radiation than their absorbed LW radiation, so they have a net radiative cooling effect on the atmosphere.

To evaluate the performance of the above surface CRKs, we calculated the surface CRE climatology during 1980−2014 in Coupled Model Intercomparison Project phase 6 (CMIP6, Eyring et al., 2016) Atmospheric Model Intercomparison Project (AMIP) simulations with the ISCCP CTP-τ cloud fraction histograms outputs from Cloud Feedback Model Intercomparison Project Observation Simulator Package (COSP, Bodas-Salcedo et al., 2011; Bony et al., 2011). The surface CRE climatology is calculated using the following equation:

where the bars denote climatological average during the reference period.

Figure 6 shows the climatological cloud cover, CRK-derived surface CRE, and actual surface CRE calculated from the CESM2 model in CMIP6-AMIP simulations (we choose CESM2 for illustration purposes because the model is widely used in climate studies). The net CRE is more negative over regions with large amounts of high clouds, the SW CRE is more negative over regions with greater total cloud fraction, and the LW CRE is more positive over regions with large amounts of low clouds. The CRE calculated from CRKs is very consistent with model-produced values, with a spatial correlation of 0.96, 0.99 and 0.98 for net CRE, SW CRE and LW CRE, respectively. Note that the ISCCP simulator contains no cloud fraction data in the winters of polar regions when there is no satellite retrieval of cloud optical thickness due to the lack of solar radiation, so we only calculated the CRE between 60°S−60°N in this study.

Figure 6.  Climatology of cloud fraction and surface CRE simulated by CESM2 in CMIP6 AMIP experiments (January 1980—December 2014). (a) Climatological high cloud (CTP <440 hPa) cover calculated from the CTP-τ cloud fraction histograms. (b) Middle cloud cover (440 < CTP < 680 hPa). (c) Low cloud cover (CTP > 680 hPa). (d−f) Surface CRE calculated with surface CRKs and CTP-τ cloud fraction histograms. (g−i) Surface CRE calculated with clear-sky and all-sky surface radiative fluxes. The spatial correlation coefficients between the CRK-derived CRE and actual values are shown in (g−i).

Figure 7 shows the observation-based cloud fraction calculated from the cloud fraction histograms of the MODIS level-3 C6.1 Monthly Global Gridded data (Hubanks et al., 2016), surface CRE calculated with MODIS cloud fraction histograms and CRKs, and surface CRE from Clouds and the Earth's Radiant Energy System (CERES) Energy Balanced and Filled (EBAF) v4.1 data (Loeb et al., 2018). Note that joint CTP-τ cloud fraction histograms for partially cloudy pixels are listed separately in MODIS level 3 data, and we added these partially cloudy pixels to the cloud fraction histograms. The climatological surface CRE calculated from MODIS+CRK is consistent with CERES values, with a spatial correlation of 0.92, 0.96 and 0.92 for net CRE, SW CRE and LW CRE between 60°S-60°N, respectively. These results indicate that the surface CRKs are can serve to explain the spatial distribution of climatological CRE for both climate models and observations.

Figure 7.  Climatology of cloud fraction and surface CRE from observations (July 2002−June 2020). The figure is the same as Fig. 6, except that the cloud fraction histograms are from MODIS C6.1 level 3 data, and the surface clear-sky and all-sky radiative fluxes are from CERES EBAF v4.1 data.

Next, we calculated the cloud induced surface radiative flux anomalies ($ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $, surface flux anomalies due to changes in cloud property anomalies) using the surface CRKs and CTP-τ cloud fraction histograms using the following equation:

where $ \Delta C $ is the difference between the monthly cloud fraction and the climatological value of the corresponding month. The values of $ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $ approximately equal the anomalies of CRE, and can be calculated by adjusting CRE anomaly with non-cloud kernels (Soden et al., 2008):

where $ {X}_{i} $ is the i’th non-cloud variable (air temperature, water vapor, surface temperature, etc.), $ {K}_{i,\mathrm{S}\mathrm{F}\mathrm{C},\mathrm{c}\mathrm{l}\mathrm{r}}^{} $ is the clear-sky surface kernel for the corresponding non-cloud variable, $ {K}_{i,\mathrm{S}\mathrm{F}\mathrm{C}}^{} $ is the all-sky surface NCRK, $ {G}_{\mathrm{S}\mathrm{F}\mathrm{C},\mathrm{c}\mathrm{l}\mathrm{r}}^{} $ is the clear-sky forcing, and $ {\mathrm{G}}_{\mathrm{S}\mathrm{F}\mathrm{C}} $ is the all-sky forcing. Considering that the change of direct forcing is small during the AMIP period, and it is difficult to estimate regional aerosol direct forcing on surface fluxes accurately, the last term of Eq. (11) was set to be zero in our calculations. We used the NCRKs from Huang et al. (2017) in the calculation of Eq. (11) and compared the values with the CRK-derived values. According to Figs. 8(a-h), net $ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $ derived from surface CRKs are highly correlated with the values calculated from surface fluxes in climate models. The correlation between regional CRK-derived $ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $ and values calculated from surface flux are high in most regions for all these models [Figs. 9(a-h)].

Figure 8.  Relationship between cloud-induced surface net radiation monthly anomalies (ΔR_sfc,cloud) derived with CRKs (y-axes) and the values derived with surface fluxes and NCRKs (x-axes). The colors denote the number of cases in each pixel. (a−h) are for individual CMIP6-AMIP models (January 1980−December 2014), and (i) is for CERES and MODIS observations (July 2002−June 2020). The black lines denote the equal-value lines, and the correlation coefficients are shown in each subplot.

Figure 9.  Correlation coefficients between regional CRK-derived and NCRK-derived net ΔR_sfc,cloud in (a−h) climate models (January 1980−December 2014) and (i) observations (July 2002−June 2020).

We also calculated the CRK-derived $ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $ with MODIS cloud fraction histograms and surface CRKs using Eq. (10) and calculated the NCRF-derived $ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $ during July 2002 to June 2020 with CERES radiative fluxes and ERA5 (Hersbach et al., 2020) meteorological fields (we are not using ERA-interim here because ERA-interim data ends in 2019). Figure 8i shows that the CRK-derived and NCRK-derived $ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $ are also highly correlated for observational data, but the correlation coefficients of the observation-based values are lower than that in climate models. One reason is that CERES does not directly measure the radiative fluxes at the surface, so that the uncertainty of CERES surface CRE is inevitable, so the correlation between CRK-derived $ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $ and NCRK-derived $ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $ is reduced by the inaccuracy of CERES surface fluxes. In addition, uncertainties regarding NCRKs and reanalysis data also reduce the correlations in Figs. 8i and 9i.

The uncertainty of CERES surface fluxes can be illustrated by the trend of CRE from July 2002 to June 2020. It is apparent from Figs. 10a-f that the trend of MODIS cloud fraction and CRK-derived $ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $ is generally continuous spatially, with no spurious signals. The SW $ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $ trends calculated with MODIS cloud data and surface CRKs are consistent with the SW CRE trends of CERES. However, geostationary satellite observations are used in the CERES surface irradiance products (Kato et al., 2018; Loeb et al., 2018), and the CERES net CRE and LW CRE trends are affected by spurious spatial discontinuities caused by the orbital tracks of geostationary satellites. According to Fig. 10i, the LW CRE trend in CERES is positive between 160°W and 90°W in the midlatitudes and is negative in the midlatitudes of 90°−180°E and 90°−30°W, indicating strong artifacts. In this case, the trends of $ {\Delta R}_{\mathrm{c}\mathrm{l}\mathrm{o}\mathrm{u}\mathrm{d},\mathrm{S}\mathrm{F}\mathrm{C}} $ calculated from CRKs and MODIS are more reliable than that calculated from CERES surface fluxes. These trends might be used to analyze the evolution of SST pattern, which has been demonstrated to be essential in the closure of the Earth’s energy budget (Zhou et al., 2021).

Figure 10.  Trends of cloud cover and surface fluxes during July 2002 and June 2020. (a−c) Trends of low, middle and high cloud cover calculated with MODIS CTP-τ histograms. (d−f) Trends of ΔR_sfc,cloud calculated with MODIS CTP-τ histograms and surface CRKs. (g−i) CRE trends calculated with CERES EBAF v4.1 data. Spurious trends associated with obit tracks can be found in (g) and (i). The unit is W m⁻² (10 yr)⁻¹.

In this study, we developed a set of surface CRKs from CBP-τ cloud fraction histograms, and a set of surface, TOA and atmospheric CRKs based on CTP-τ histograms. The surface CRKs from CBP-τ histograms and TOA CRKs from CTP-τ cloud fraction histograms were calculated from RRTM simulations, while the statistical relationship between CTP, CBP and τ from collocated CloudSat-MODIS data was used in the calculation of the surface and atmospheric CRKs from the CTP-τ histograms.

We calculated the climatological CRE and monthly ΔR_cloud anomalies with modeled and observed CTP-τ histograms, and these CRK-derived values agree well with NCRK-derived values. The CRKs can be used to attribute ΔR_cloud to changes in cloud properties, which would be useful in the analyses of surface cloud radiative feedback processes. Our separate ice/liquid CRKs can be used to analyze the cloud phase feedback with ice/liquid cloud fraction histograms. The CRKs can also be used to diagnose the CRE changes in surface radiative products (for example, the spurious surface LW CRE trend in CERE EBAF v4.1 data).

It is worth noting that Zhang et al. (2021) were also developing a set of surface and TOA CRKs using different methods when our analyses were carried out, so it is worthwhile to discuss the differences between these two studies. The kernels in Zhang et al. (2021) were calculated using the radiation code of a global climate model (GISS GCM), while the radiative fluxes in our study were calculated from a stand-alone radiative transfer model (RRTM). The cloud vertical structure of Zhang et al. (2021) was based on a combination of rawinsonde climatology and CloudSat-CALIPSO climatology, while the statistical relationships between CTP, CBP and τ in our study were calculated from collocated MODIS-CloudSat climatology. The surface CRKs in Zhang et al. (2021) were not evaluated with independent data sources, so our evaluation of surface CRKs in reproducing the climatology and anomalies of cloud radiative effects is unique. Our surface CRKs on CBP-τ histograms and separate CRKs for liquid/ice clouds are also unique. Our surface CRKs on CTP-τ histograms are generally consistent with Zhang et al. (2021), which reinforces the robustness of the surface CRKs in both studies.

The surface CRKs on CBP-τ histograms are more accurate than surface CRKs on CTP-τ histograms because the surface CRE is more affected by CBP than CTP, but CBP-τ histograms product can be hardly found at present. Our study encourages the production of CBP-τ cloud fraction histograms for both climate models and observations, which can be used to diagnose the climatology and changes of surface CRE more accurately than CTP-τ histograms.

Acknowledgements. We thank Dr. M. ZELINKA for valuable suggestions and thank Dr. Y. ZHANG for useful discussions. This work is supported by the National Natural Science Foundation of China (Grant No. NSFC 41875095 and 42075127). The CRK datasets are available online at Zenodo, doi: https://doi.org/10.5281/zenodo.4732640

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author (s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.