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

Cloud radiative kernels (CRK) built with radiative transfer models have been widely used to analyze the cloud radiative effect on top of atmosphere (TOA) fluxes, and it is expected that the CRKs would also be useful in the analyses of surface radiative fluxes, which determines the regional surface temperature change and variability. In this study, CRKs at the surface and TOA were built using the Rapid Radiative Transfer Model (RRTM). Longwave cloud radiative effect (CRE) at the surface is primarily driven by cloud base properties, while TOA CRE is primarily decided by cloud top properties. For this reason, the standard version of surface CRK is a function of latitude, longitude, month, cloud optical thickness (τ) and cloud base pressure (CBP), and the TOA CRK is a function of latitude, longitude, month, τ and cloud top pressure (CTP). Considering that the cloud property histograms provided by climate models are functions of CTP instead of CBP at present, the surface CRKs on CBP-τ histograms were converted to CTP-τ fields using the statistical relationship between CTP, CBP and τ obtained from collocated CloudSat and MODIS observations. For both climate model outputs and satellites observations, the climatology of surface CRE and cloud-induced surface radiative anomalies calculated with the surface CRKs and cloud property histograms are well correlated with those calculated from surface radiative fluxes. The cloud-induced surface radiative anomalies reproduced by surface CRKs and MODIS cloud property histograms are not affected by spurious trends that appear in Clouds and the Earth’s Radiant Energy System (CERES) surface irradiances products.


Introduction
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 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. (2012bZelinka et al. ( , 2016Zelinka et al. ( , 2020 attributed the changes in ΔR cloud,TOA under global warming to changes in cloud amount, cloud optical depth and CTP. Zhou et al. (2013Zhou et al. ( , 2014 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.

Surface CRKs built on CBP-τ histograms
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 2 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.
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 is the surface radiative flux, is the surface radiative flux for cloudy sky, and is the clearsky 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.
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 irra-diance 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 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.
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 τ.

Surface, TOA and atmospheric CRK on CTP-τ histograms
The surface CRE can be calculated with CRKs on either CBP-τ or CTP-τ cloud fraction histograms: where C 1 and C 2 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 , then we have Combining Eqs. (5−7) we have We used the collocated CloudSat-MODIS data to calcu- late the probability function 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 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 values averaged over regions with similar cloud features (classified by low cloud fraction with CTP > 680 hPa). If regional mean 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
, 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.
We also calculated a set of TOA CRKs using the clearsky 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.

Calculation of the climatology and monthly anomalies of surface CRE
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, CRKderived 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 calcu- lated the CRE between 60°S−60°N in this study. 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. ∆R cloud,SFC Next, we calculated the cloud induced surface radiative flux anomalies ( , surface flux anomalies due to changes in cloud property anomalies) using the surface CRKs and CTP-τ cloud fraction histograms using the following equation: ∆C ∆R cloud,SFC where is the difference between the monthly cloud fraction and the climatological value of the corresponding month. The values of approximately equal the anomalies of CRE, and can be calculated by adjusting CRE anomaly with non-cloud kernels (Soden et al., 2008): where is the i'th non-cloud variable (air temperature, water vapor, surface temperature, etc.), is the clearsky surface kernel for the corresponding non-cloud variable, is the all-sky surface NCRK, is the clear-sky forcing, and 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 derived from surface CRKs are highly correlated with the values calculated from surface fluxes in climate models. The correlation between regional We also calculated the CRK-derived with MODIS cloud fraction histograms and surface CRKs using Eq. (10) and calculated the NCRF-derived 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 NCRKderived are also highly correlated for observational data, but the correlation coefficients of the observationbased 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 and NCRK-derived 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.  (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 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).

Conclusions and discussions
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-Cloud-Sat 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.