## Abstract

The Clouds and the Earth’s Radiant Energy System (CERES)–partial radiative perturbation [PRP (CERES-PRP)] methodology applies partial-radiative-perturbation-like calculations to observational datasets to directly isolate the individual cloud, atmospheric, and surface property contributions to the variability of the radiation budget. The results of these calculations can further be used to construct radiative kernels. A suite of monthly mean observation-based inputs are used for the radiative transfer, including cloud properties from either the diurnally resolved passive-sensor-based CERES synoptic (SYN) data or the combination of the *CloudSat* cloud radar and *Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations* (*CALIPSO*) lidar. The *CloudSat*/*CALIPSO* cloud profiles are incorporated via a clustering method that obtains monthly mean cloud properties suitable for accurate radiative transfer calculations. The computed fluxes are validated using the TOA fluxes observed by CERES. Applications of the CERES-PRP methodology are demonstrated by computing the individual contributions to the variability of the radiation budget over multiple years and by deriving water vapor radiative kernels. The calculations for the former are used to show that an approximately linear decomposition of the total flux anomalies is achieved. The observation-based water vapor kernels were used to investigate the accuracy of the GCM-based NCAR CAM3.0 water vapor kernel. Differences between our observation-based kernel and the NCAR one are marginally larger than those inferred by previous comparisons among different GCM kernels.

## 1. Introduction

Understanding the variability and response of Earth’s radiation budget to perturbations requires isolating individual contributions from each component of the atmosphere and surface. Isolating radiative effects has been strongly motivated by the need to evaluate feedbacks—processes that can amplify or dampen the response to a radiative forcing—in global climate models (GCMs). A feedback is defined as the change in top-of-the-atmosphere (TOA) radiative flux for a given change in surface temperature. To isolate the change in TOA flux caused by a particular variable, Wetherald and Manabe (1988) introduced the concept of partial radiative perturbation (PRP) calculations. A PRP calculation substitutes a single variable from the perturbed climate into the control climate to compute the resulting perturbation to the TOA flux. Alternatively, the radiative kernel technique (Held and Soden 2000; Soden and Held 2006; Soden et al. 2008; Shell et al. 2008) separates the radiative response and the perturbation. A radiative kernel is constructed by computing the TOA radiative effect from a small perturbation to a variable in the control climate. When multiplied by the perturbation in the variable itself, the kernel then gives the resulting perturbation to the TOA flux.

Radiative kernels require far fewer overall calculations than the PRP method and can be consistently and easily applied across different climate models. Many extensions and adaptations of the radiative kernel technique have been made (e.g., Jonko et al. 2012; Sanderson and Shell 2012; Zelinka et al. 2012; Held and Shell 2012; Block and Mauritsen 2013) and kernels are an indispensable tool for assessing climate feedbacks in GCMs (e.g., Vial et al. 2013). However, radiative kernels become increasingly inaccurate when applied to climates dissimilar from that in which they were created (Jonko et al. 2012; Block and Mauritsen 2013). Therefore, differences between the control climate of a GCM and the observed climate can introduce errors. Despite this, only a limited effort has been made to validate the accuracy of GCM radiative kernels relative to observation-based ones (Gordon et al. 2013; Yue et al. 2016).

GCM-derived radiative kernels have also been applied to observations in an attempt to constrain climate feedbacks simulated by GCMs (Dessler et al. 2008; Dessler and Wong 2009; Dessler 2010; Masters 2012; Dessler and Loeb 2013; Dessler 2013; Gordon et al. 2013; Zhou et al. 2013, 2014; Ceppi et al. 2016). Such studies are challenging since it is not yet clear how feedbacks in the current climate relate to those under a long-term global warming scenario (Trenberth et al. 2010; Dessler 2010, 2013; Zhou et al. 2015; Chung et al. 2012; Gordon et al. 2013; Colman and Hanson 2017). Regardless, there is a concerted effort to determine how observations can best be used to constrain climate feedbacks and sensitivity (Loeb et al. 2016; Forster 2016).

Given the widespread use of GCM-derived radiative kernels, the exercise of validating them versus observations is a worthwhile one. Additionally, for examining feedbacks in the observational record, it would be attractive to perform calculations based purely on observations. The concepts underlying PRP and radiative kernel calculations also have value in decomposing the observed variability of Earth’s radiation budget. While data from the Clouds and the Earth’s Radiant Energy System (CERES; Wielicki et al. 1996) instruments provide observations of the total radiative fluxes, knowledge of the contributions from individual variables can provide valuable insight into the process being studied. All of this motivates the work presented in this study. Here we were apply PRP-like calculations to observational datasets to directly isolate the individual contributions to the variability of radiation budget. The results of these calculations can then be further used to construct radiative kernels.

Observations are mostly taken from the various datasets used in processing CERES edition 4 data, and we therefore refer to the dataset as CERES-PRP. While similar calculations using observations have been performed in other studies (e.g., Loeb et al. 2012; Crook and Forster 2014), such work has been limited to a fixed set of parameters since the desired decomposition of the radiation budget depends on the question being investigated. Here we develop a more general approach and perform PRP-like calculations in a consistent framework to accommodate the calculation of flux perturbations due to any desired combination/subset of variables. This flexibility allows the CERES-PRP methodology to generate datasets at will and tailored to the specific question being investigated. To achieve this goal, CERES-PRP needs to be computationally efficient, and this is mainly accomplished through the use of monthly mean inputs for the radiative flux calculations. Given the nonlinearity of using time-averaged inputs, monthly averages are carefully constructed and the computed fluxes are validated relative to observed values. In particular, a clustering method is developed to allow for the use of active sensor cloud properties from the *CloudSat* cloud radar (Stephens et al. 2008) and *Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations* (*CALIPSO*) lidar (Winker et al. 2009, 2010).

While the application of radiative kernels takes a trivial amount of effort, the initial computations are time consuming. This often leads to GCM feedbacks being evaluated using kernels from outdated models (e.g., Vial et al. 2013) and the typical practice of using a single yearlong control run to compute kernels (e.g., Soden et al. 2008; Shell et al. 2008). The approach of CERES-PRP helps diminish these computational barriers. This also enables timely updates as new (and presumably improved) observational datasets become available and allows for calculations to be repeated using different sources of data to help quantify uncertainties, as is done in this study.

Descriptions of the inputs and radiative transfer model are given in sections 2 and 3, respectively. Section 4 reviews the PRP and kernel methods and discusses their implementation in CERES-PRP. Validation of the computed fluxes are provided in section 5, including assessing the efficacy of using monthly mean inputs and the *CloudSat*/*CALIPSO* clustering method, along with comparisons to the CERES-observed fluxes. A demonstration of decomposing the radiation budget is shown in section 6, and those calculations are used to test the linearity of the perturbations in section 6a. CERES-PRP calculations are then used to examine the accuracy of the widely used GCM-based NCAR CAM3.0 water vapor kernel (Shell et al. 2008) in section 6b.

## 2. Input datasets

Inputs for the radiative transfer calculations are mostly taken from the various observational datasets used in the processing of CERES edition 4 data products. Monthly mean values are used for all inputs, and vertical profiles are interpolated to 30 fixed pressure levels. An additional level is set at the monthly mean surface pressure for each grid box. All inputs are regridded, if necessary, to either a 1° × 1° (latitude × longitude) or 2° × 5° grid, depending on the source of the cloud properties as explained below.

Profiles of temperature, water vapor, and ozone along with surface pressure and skin temperature are taken from either the Goddard Earth Observing System (GEOS), version 5.4.1, reanalysis (Rienecker et al. 2008) or the Atmospheric Infrared Sounder (AIRS), version 6 level 3, standard product (Chahine et al. 2006; Tian et al. 2013). When using the AIRS data, the upper atmosphere above 100 (1) mb for water vapor (temperature and ozone) is filled in with GEOS data, since the AIRS product does not provide retrievals at these pressure levels. Methane concentrations are also taken from the AIRS, version 6 level 3, standard product (Xiong et al. 2008). Carbon dioxide concentrations are obtained from the AIRS, version 5 level 3, carbon dioxide product (Chahine et al. 2005). For months before when AIRS data are available (i.e., before September 2002), global means from the NOAA Earth System Research Laboratory (ESRL) Global Monitoring Division data are used (Conway et al. 1994; Dlugokencky et al. 2009). Other trace gases—nitrous oxide, CFC-11, CFC-12, and HCFC-22—are obtained from global mean NOAA ESRL data as well.

The Model of Atmospheric Transport and Chemistry (MATCH) MODIS aerosol assimilation product (Collins et al. 2001) provides aerosol optical depths (at 0.55 *μ*m) and vertical distributions for seven different types of aerosols. The spectral surface emissivities are set to fixed values based on the International Geosphere–Biosphere Programme (IGBP) surface classification system (Wilber et al. 1999). Spectral surface albedos are determined using surface-type-based lookup tables (Jin et al. 2004; Rutan et al. 2009) for the spectral shapes that are scaled by the broadband surface albedo determined from the CERES surface albedo history (SAH) map (Rutan et al. 2009). The TOA solar irradiance is taken from the Solar Radiation and Climate Experiment (SORCE) level 3 total solar irradiance (TSI) data product (Kopp and Lean 2011).

### a. Cloud properties

Cloud properties—cloud fraction, base pressure, top pressure, phase, optical depth, ice effective diameter, and liquid effective radius—are taken from two main sources. The CERES synoptic (SYN) edition 4 product provides diurnally resolved cloud properties by combining retrievals from *Terra* MODIS, *Aqua* MODIS, and geostationary satellites (Minnis et al. 2011a,b; Doelling et al. 2016; P. Minnis et al. 2018, unpublished manuscript). The second source of cloud properties is the combined CERES–*CALIPSO*–*CloudSat*–MODIS (C3M) product Kato et al. (2010). Additionally, CERES-PRP calculations are performed using the CERES *Aqua* single scanner footprint (SSF) product MODIS cloud retrievals (Minnis et al. 2011a,b) and the subset of the CERES *Aqua* SSF cloud properties collocated with the *CALIPSO* and *CloudSat* ground tracks (SSFC3M). These additional two sets of SSF cloud properties are used to investigate sampling differences between the SYN and C3M cloud properties.

The choice of cloud properties in CERES-PRP determines the nominal spatial resolution. For SYN and SSF, a 1° × 1° (latitude × longitude) grid is used. Since *CALIPSO* and *CloudSat* have no cross-track sampling, a 2° × 5° (latitude × longitude) grid is used in the same manner as the *CALIPSO* level 3 products (Winker et al. 2013). When using C3M cloud properties, no calculations are made poleward of the 81°N or 81°S grid box, since *CloudSat* and *CALIPSO* do not sample these latitudes.

Clouds are organized into monthly mean “cloud conditions,” that is, sets of nonoverlapping clouds that, when weighted by the cloud fraction, can be summed together to obtain the total flux in each grid box. This linear separation allows for the isolation of flux perturbations from each cloud condition. The SYN cloud properties are already organized as such since the dataset only contains single-layer clouds averaged into four categories: clouds with tops below 700 mb, tops between 500 and 700 mb, tops between 300 and 500 mb, and tops above 300 mb. We extend the nominal 4 SYN categories to 12 by further separating each cloud-top category into 3 different optical depth categories: optical depths less than 1, between 1 and 5, and greater than 5. The addition of these three optical depth categories helps to improve the linearity of the radiative response to cloud perturbations. The residuals from nonlinear effects are explored further in section 6a.

As described by Kato et al. (2010), clouds in the C3M product are obtained by merging similar *CloudSat* and *CALIPSO* cloud profiles contained within a single CERES footprint. This merging creates a maximum of 16 cloud conditions per CERES footprint, each with up to six overlapping layers. This concept is extended to the monthly time scale using *k*-means clustering (Arthur and Vassilvitskii 2007). The *k*-means clustering is applied separately to sets of properties from cloud conditions with the same number of overlapping layers (from one to six layers). We seek to partition a set of observations in sets

where is the chosen number of clusters. The cloud base, top, and natural log of the optical depth are taken as the observational set, for example, , where *x* are vectors containing all C3M single-layer cloud conditions that occur in a single grid box during a particular month. For multiple overlapping cloud layers, the observational set becomes , where is the number of cloud layers.

The *k*-means algorithm finds cluster centers , that is, the mean point in , by minimizing

where CF is the cloud fraction. The cloud properties vectors *x* are normalized by their standard deviations so that all three cloud properties contribute equally to the Euclidean distance. The factor is included in Eq. (2) to make cloud conditions with a larger spatial coverage (i.e., cloud fraction) have more influence over the cluster center. Separate clusters are found for each 2° × 5° grid box and month of the year (January–December). Cluster definitions remained fixed for each individual month in the time series; for example, clusters for January are obtained by using observations from all Januaries during the C3M time period (July 2006–December 2010).

The number of clusters is determined separately for each grid box and month by determining the fractional occurrence of the number of overlapping cloud layers (i.e., frequency of 1-layer cloud conditions, 2-layer, …, 6-layer). A total number of clusters *N* is chosen and allocated using the fractional occurrence of the number of cloud layers so . In this study, we chose *N* = 40 so, for example, if for a particular grid box and month, 50% of the C3M cloud conditions had 1-layer clouds, 30% had 2-layer clouds, and 20% had 3-layer clouds then *n*_{1}, *n*_{2}, and *n*_{3} would be set to 20, 12, and 8, respectively. The choice of *N* = 40 represents a balance of accuracy and computational time. Errors as a function of the number of clusters are explored further in section 5a.

Once the cluster centers are determined, they are applied to each month of C3M data by using a nearest-neighbor search and the Euclidean distance to match each individual cloud condition to its closest cluster. The matches for each cluster are then averaged, weighted by the cloud fraction, to obtain the monthly mean C3M cloud properties for up to 40 different cloud conditions.

By organizing the SYN and C3M observations into these cloud conditions, we seek to organize clouds “radiatively” to minimize the nonlinear response of fluxes to perturbations within the definitions of each category. This is essentially the same approach used by Zelinka et al. (2012), who derived cloud radiative kernels from cloud-top and optical depth classifications. Here our SYN clouds are organized into 12 categories, while the C3M cloud conditions are much more complex with 40 categories that vary both spatially and seasonally.

### b. Diurnal cycle

For most of the intended applications of the CERES-PRP calculations, diurnally resolved fluxes are desired. While varying the solar zenith angle in the radiative transfer calculation is simple enough, supplying diurnally varying observations for those calculations is challenging since more accurate global datasets are typically from instruments on polar-orbiting platforms. The CERES SYN product’s novel use of geostationary observations addresses this shortcoming by providing hourly cloud retrievals. To reduce the computation time, we subset the hourly SYN data into 3-hourly data and average to obtain 3-hourly monthly mean cloud properties. The monthly mean diurnal cycle of the meteorology (i.e., temperature, water vapor, and ozone) is also taken into account using 3-hourly GEOS reanalysis data. The remainder of the datasets we desire diurnal information from (AIRS, C3M, and *Aqua* SSF) are all in sun-synchronous polar orbits and only provide single daytime/nighttime views of a fixed location. Monthly mean values from these datasets are computed separately for daytime and nighttime observations.

To mix and match datasets with different diurnal resolutions, the day fraction, that is, the fraction of the time when the sun is above the horizon, at each grid point and for each hour in the 3-hourly data is used to compute a weighted averaged of the daytime and nighttime monthly means. In the case when both datasets only have single daytime and nighttime observation (e.g., a calculation that combines C3M and AIRS), multiple daytime solar zenith angles are used to allow the sun to vary, while the daytime properties stayed fixed. In this work, we use four solar zenith angles that evenly sample the distribution of solar zenith angles that occur in each grid box during that month when C3M and AIRS are combined.

CERES-PRP can also interpolate and/or average the 3-hourly datasets as needed for investigating the effects of sampling. In this study, we will make use of the SYN data that have been interpolated to *Aqua*’s overpass times and then averaged into single daytime and nighttime observations, which mimics the sampling of C3M/AIRS/SSF.

## 3. Radiative transfer model

Radiative transfer calculations are performed using the NASA Langley Fu–Liou radiative transfer model. A summary of the inputs needed are given in Table 1, and a full description of the model can be found in Rose et al. (2013). In brief, the model is based on Fu and Liou (1992, 1993) with updates to the shortwave (SW) and longwave (LW) gaseous absorption (Kato et al. 1999; Kratz and Rose 1999). The effects of ice cloud crystals are included using Fu (1996) and Fu et al. (1998b) and liquid cloud droplets using Hu and Stamnes (1993). Aerosol-scattering properties for each of the seven MATCH aerosol types are determined using Hess et al. (1998) or Tegen and Lacis (1996) and Sinyuk et al. (2003) for dust. The distribution of cloud optical depths are assumed to be a gamma distribution that is estimated by taking the linear and logarithmic mean optical depth (Barker 1996; Oreopoulos and Barker 1999; Kato et al. 2005). These linear and logarithm mean optical depths of each layer are used in a gamma-weighted 2-stream shortwave solver (Kato et al. 2005). The logarithmic mean cloud optical depth is used in the longwave flux computation using a hybrid 2-/4-stream approach (Fu et al. 1998a).

## 4. PRP calculations

Wetherald and Manabe (1988) introduced the concept of PRP calculations to isolate the change in TOA flux caused by a particular variable. The application of a PRP calculation is straightforward: A single variable from the perturbed climate is substituted into the control climate and a radiative transfer calculation is performed to determine the effect on the flux. In essence, PRP calculations are just a rebranding of finite differences. The effect on the flux due to some perturbation of variable *x* is

where the fluxes *F* are at some level of the atmosphere (e.g., the TOA or surface). In the traditional PRP approach, the instantaneous values of variables are used rather than the time-averaged monthly means used in this study. Therefore, in this adaptation, is the deseasonalized anomaly in the monthly mean value *x* relative to the climatological monthly mean , that is, . The climatological monthly mean is the mean value over all months in the time series; for example, for January is obtained by averaging values from all Januaries in the full time series. The climatological monthly mean values of all other variables relevant to the radiative transfer besides (i.e., ) are used in both flux calculations. The superscript *f* denotes that Eq. (3) is a forward finite difference, which has a truncation error of . This difference is also referred to as a “one sided” PRP calculation (Soden et al. 2008). The subscript *C* denotes that the flux perturbation is relative to the climatological monthly mean base state . Alternatively, the effect of the perturbation can be computed using the backward finite difference,

where the superscript *b* is used to denote the backward finite difference and the subscript *M* indicates that the flux perturbation is relative to the monthly mean base state .

The problem with both Eqs. (3) and (4) is that the truncation errors are the same order as the perturbation itself. This results in calculating not only a flux perturbation due to , but also one that results from a decorrelation between the perturbed variable and nonperturbed variables. This additional undesirable term can dominate the perturbation calculation (Soden et al. 2008). To minimize this, an improved estimate of the radiative perturbation can be computed using the centered finite difference,

whose truncation error is improved to second order. The centered difference is equivalent to averaging the forward [Eq. (3)] and backward [Eq. (4)] differences and is also referred to as the “two sided” PRP (Colman and McAvaney 1997). As an alternative to the two-sided PRP technique, the radiative kernel technique (Held and Soden 2000; Soden and Held 2006; Soden et al. 2008; Shell et al. 2008) uses Eq. (3) and fixes as a small perturbation in order to minimize the truncation error to . Taking the ratio

gives a radiative kernel that can then be multiplied by a larger caused by, for example, a change in the climate to obtain the corresponding radiative effect.

### a. Implementation in CERES-PRP

CERES-PRP seeks to decompose the total flux anomaly using finite differences (i.e., PRP calculations) rather than using the radiative kernel technique. This allows for more flexible calculations to be performed, since the radiative kernel technique requires that the flux perturbation and the perturbation itself be linearly related for the ratio in Eq. (6) to be valid. If one, for example, wanted to compute the effect on the TOA flux from a perturbation in water vapor, separate flux calculations for each pressure level are needed to construct a kernel that gives the impact of each level’s perturbation on the TOA flux. This kernel profile is then multiplied by the observed/simulated water vapor anomaly profile and summed to give the total effect on the TOA flux. In contrast, in the finite-difference framework, CERES-PRP can substitute the entire water vapor profile at once and compute the combined impact on the TOA flux in one step. This flexibility also allows for computing radiative perturbations that result from any combination of the variables in Table 1, for example, different variables (e.g., temperature and water vapor), specific dimension(s) of a variable (e.g., stratospheric water vapor or an individual aerosol type), or total effect of clouds by perturbing all cloud properties. Additionally, CERES-PRP can examine the resulting flux perturbations at the TOA, surface, or any other level of interest in the atmosphere.

While an effort is made to facilitate using monthly mean inputs for the flux calculations and quantifying the uncertainties this causes (section 5), Eq. (5) introduces a flux calculation with yet another time average that uses climatological monthly mean inputs . This diminishes the accuracy of Eq. (5), since the fluxes computed relative to the monthly () and climatological means will have different time-average-related biases. To avoid this, we first note that the centered difference can also be computed by averaging

and

Eqs. (7) and (8) compose a “wide” two-sided PRP calculation, since these equations contain values beyond the nominal range to *x*, from to . Since Eq. (7) also relies on a flux calculation that uses climatological time-averaged inputs, it is more desirable to compute the centered difference by averaging Eqs. (4) and (8). However, when applying Eq. (8), it is possible for the substitution to be unphysical and therefore not usable in the radiative transfer model, for example, a perturbed surface albedo below (above) 0 (1) or an ice effective diameter outside the radiative transfer model’s parameterization range. Additionally, for any of the finite differences, it is possible to encounter unphysical values when perturbing cloud boundaries, for example, a perturbed cloud top that is below its corresponding base or a perturbed cloud base below the surface.

CERES-PRP follows the logic outlined in Fig. 1 to obtain an estimate of the centered difference that avoids unphysical perturbations and minimizes the use of Eqs. (3) and (7). CERES-PRP starts by attempting to calculate the forward difference using the month base state [Eq. (8)]. If that perturbation is unphysical, then the forward difference using the climatological base state is computed instead [Eq. (3)], which, with the exception of cloud boundary perturbations, is guaranteed to contain physically valid values. In practice, we find that Eq. (3) is used in lieu of Eq. (8) from 0% of calculations for temperature perturbations to a maximum of 12% of calculations for surface albedo perturbations. When perturbing cloud boundaries, it is still possible for the evaluation of Eq. (3) to be unphysical. In these cases, the cloud boundary perturbation is adjusted to be as close as possible to the true perturbation without being unphysical and both Eqs. (8) and (3) are evaluated. We find that these adjustments are needed for approximately 4% of calculations. The same hierarchy of decisions is repeated for the backward difference. All valid calculations are then averaged together to obtain the final estimate of the centered difference .

### b. CERES-PRP radiative kernels

While computing the flux perturbations directly allows for flexibility, computing flux perturbations indirectly using a radiative kernel has a number of advantages. By precomputing the radiative effect, a radiative kernel can easily and consistently be applied (e.g., to assess feedbacks among GCMs) without the need for additional radiative transfer calculations. Instead of using a fixed perturbation, CERES-PRP computes kernels as a by-product of the centered difference calculations using the observed anomaly of the variable:

Multiple calculations are performed as needed to avoid nonlinearity in the ratio. For the example of water vapor, as described above, separate CERES-PRP calculations would be performed to obtain the effect on the TOA flux from each individual water vapor level [i.e., compute Eq. (9) for each level]. A complete kernel can then be obtained by compiling the from each level. We expect to be small enough to give an acceptably small truncation error since the error is minimized by using the centered difference and the typically small size of perturbations in the observational record. In section 6b, we compare the size of the perturbation used in the CERES-PRP water vapor kernel to the standard perturbations used in GCM-based kernels.

As described above, the cloud inputs in CERES-PRP are organized into cloud conditions, nonoverlapping profiles whose fluxes can be linearly combined using cloud fraction. Therefore, CERES-PRP can compute radiative kernels for individual cloud properties by perturbing each cloud layer separately. A cloud fraction kernel can be obtained directly from the monthly mean fluxes by simply dividing the fluxes for each cloud condition by 100%. Radiative kernels, specifically temperature, water vapor, cloud fraction, and cloud-base kernels, computed by CERES-PRP, are used in edition 4 of the Energy Balanced and Filled (EBAF)-Surface product (Kato et al. 2018) to assist in adjusting computed surface fluxes for various bias corrections.

## 5. Validation of computed fluxes

In this section, the efficacy of using monthly mean inputs for radiative transfer calculations is tested along with the C3M cloud clustering method. A flux closure exercise is performed by comparing computed monthly fluxes and their variability to edition 4 of the CERES EBAF-TOA observations (Loeb et al. 2018a). EBAF-TOA improves the accuracy of the CERES-derived fluxes via a one-time adjustment to the net flux using an in situ–derived value from Johnson et al. (2016) and infers clear-sky fluxes from partly cloudy CERES footprints. Four main dataset combinations are used to explore the impact of varying the cloud and meteorological properties: 1) SYN cloud properties and GEOS meteorological properties (SYN + GEOS), 2) C3M cloud properties and GEOS meteorological properties (C3M + GEOS), 3) SYN cloud properties and AIRS meteorological properties (SYN + AIRS), and 4) C3M cloud properties and AIRS meteorological properties (C3M + AIRS). Given the datasets used in CERES-PRP, three different time periods can be examined: 1) the *Terra* period, March 2000–February 2016, 2) the *Aqua* period, September 2002–February 2016, and 3) the *CloudSat*–*CALIPSO* period, July 2006–December 2010. Although both *CloudSat* and *CALIPSO* operated together after December 2010, we limit our analysis to before *CloudSat* began daytime-only operations. Unless otherwise noted, all the following results are for the *CloudSat*–*CALIPSO* period, since this allows for comparisons among all datasets used. In general, results derived using the *CloudSat*–*CALIPSO* period are very similar to those from either the *Terra* or *Aqua* period. Any differences are noted where applicable.

### a. Flux errors from monthly mean inputs

To expedite the CERES-PRP calculations, monthly mean inputs are used for the radiative transfer calculations despite the uncertainty from the nonlinear relationship between fluxes and the inputs used to compute them. An attempt to lessen this uncertainty is made by including a monthly mean representation of the diurnal cycle (section 2b) and, in the shortwave, using a gamma-weighted 2-stream solver to better capture the effective distribution of cloud optical depths within a month (section 3). Our SYN + GEOS calculation is essentially a monthly mean version of the computed fluxes given in the CERES SYN product, which are calculated using the same radiative transfer model (section 3) but at the nominal SYN hourly resolution. Therefore, we quantify the error due to monthly mean meteorological properties by comparing the clear-sky fluxes between the monthly SYN + GEOS calculation and computed fluxes averaged from the SYN product. Similarly, to assess the effect of averaging the SYN cloud properties, we compare the cloud radiative effect (CRE; the difference between the all-sky and clear-sky fluxes) as well. Since the C3M product only provides instantaneous fluxes at the *Aqua* overpass times, we perform our own diurnally averaged calculation by repeating the flux calculation in each of the C3M’s nominal cloud groups for every 3-hourly monthly mean solar zenith angle and day fraction. The effect of averaging the C3M cloud properties before performing the radiative transfer can then be assessed by comparing to the C3M + GEOS monthly calculation.

The global mean root-mean-square errors (RMSEs) and bias errors due to the use of monthly mean meteorological and cloud properties averaged over all months in the time series are given in Figs. 2 and 3 for the TOA and surface, respectively. Biases (RMSEs) are less than about 1.7 (3.0) W m^{−2} except for clear-sky longwave errors at the surface. Relative to their respective reference calculations, errors in the TOA fluxes are all less than 6%. At the surface, both the use of monthly mean meteorological properties and C3M clouds have errors less than 2%, while using monthly mean SYN clouds causes errors less than 8%. Overall, we consider the errors in Figs. 2 and 3 to be a reasonable trade-off for the large decrease in computational time achieved by using monthly mean inputs.

As described in section 2a, the monthly mean C3M cloud properties are determined using a *k*-means clustering method that requires choosing the total number of clusters. To determine this, monthly averages of the instantaneous all-sky fluxes from the C3M product are compared to those computed using the monthly mean C3M cloud inputs. An instantaneous flux calculation is performed by both interpolating the monthly mean GEOS properties to the times of *Aqua* overpasses and by setting the day fraction and solar zenith angle to the monthly mean values from the C3M product. This comparison is repeated while varying the total number of clusters from 5 to 320.

Figure 4 shows the global mean RMSE and bias errors relative to the all-sky fluxes in the C3M product averaged over all months. Remarkably, using a mere five clusters in our relatively simple clustering method reproduces the true fluxes fairly well with errors less than about 4 W m^{−2} (about 1.6%). In this study, we chose to use 40 C3M cloud clusters, which reduced errors to less than 1.5 W m^{−2} (less than 0.75%), as a balance of reducing errors and improving computational time. Compared to the C3M product, using our 40 cloud clusters reduces the number of calculations needed for monthly mean fluxes by about a factor of 30.

### b. Comparison to CERES EBAF

Comparison of the global mean fluxes from CERES EBAF-TOA and our computed fluxes from the four different dataset combinations is shown in Fig. 5. The largest biases, regardless of datasets used, are generally found in the clear-sky comparison (Fig. 5b). This could in part be due to the differing definitions of clear sky. In EBAF-TOA, clear-sky fluxes are averaged from MODIS pixels determined to be cloud free via a cloud mask. For the CERES-PRP calculation, the clear-sky fluxes are computed regardless of the presence of clouds by removing clouds and recomputing the fluxes. From Kato et al. (2013), we expect this difference in clear-sky definition to cause a 0.24 and −1.25 W m^{−2} global mean bias in the shortwave and longwave fluxes, respectively. The biases in Fig. 5b are larger than this, so we do not expect the clear-sky definition to be the dominant cause of the differences. The use of AIRS instead of GEOS reduces clear-sky longwave flux errors by about a factor of 2. The persistent clear-sky shortwave bias suggests that CERES-PRP could benefit from an improved representation of the surface albedo, although some of this bias can be attributed to the use of monthly mean inputs as shown above (section 5a).

For the all-sky fluxes (Fig. 5a), better agreement with EBAF-TOA is achieved using the SYN cloud properties than C3M. This is somewhat counterintuitive as C3M is expected to provide more accurate cloud properties and showed smaller errors from monthly averaging (section 5a). To more directly isolate the choice of cloud properties, the comparison of CRE is given in Fig. 5c. There it is shown that using C3M does indeed provide a better representation of clouds with better agreement to EBAF-TOA. Consistent with Kato et al. (2011), using C3M reduces both the shortwave and longwave upwelling flux at the TOA. Taking Figs. 5a–c as a whole shows that the better agreement in all-sky fluxes when using the SYN cloud properties is due to a fortunate cancellation of clear-sky and SYN cloud biases.

One of the goals of this work is to decompose the variability of the radiation budget. In that respect, it is more important that the variability of EBAF-TOA is reproduced rather than the absolute value of the fluxes. Figure 6 shows the correlation coefficients between the EBAF-TOA zonal-mean time series and the computed fluxes using the four different dataset combinations. The choice of GEOS (solid lines) or AIRS (dashed lines) has little effect on the agreement between the time series. Using the SYN cloud properties produces fluxes that correlate well with EBAF-TOA (*r* typically above 0.90) with lower correlations toward the poles, particularly for the Southern Hemisphere shortwave fluxes. At all latitudes, fluxes computed with the SYN cloud properties have a higher correlation with EBAF-TOA than when using C3M cloud properties.

The poorer C3M correlations in Fig. 6 are likely due to a difference in sampling since, relative to SYN cloud properties, C3M lacks a full representation of the diurnal cycle and is spatially limited to nadir-only sampling. To correct the C3M calculation for these effects, fluxes are computed using SYN cloud properties that have been interpolated to *Aqua*’s overpass times and then averaged into single daytime and nighttime observations. The ratio of this interpolated/averaged SYN calculation to the nominal SYN calculation gives a multiplicative factor for each grid point that can be used to correct for the lack of the cloud diurnal cycle in the C3M calculations. In a similar manner, the ratio of the fluxes computed with the SSF cloud properties to those computed using SSF cloud properties from only the MODIS pixels along the C3M ground track is determined. This ratio is used to correct for the nadir-only sampling of C3M. The correlations with EBAF-TOA of the zonal-mean time series using these scaled C3M calculations are shown in Fig. 7. In Fig. 7, all calculations use the GEOS properties since results using AIRS are very similar (not shown). Figure 7 shows that scaling the C3M-computed fluxes by the SYN cloud diurnal cycle causes little change in the correlations. Scaling to correct for the nadir-only sampling generally improves the correlations. The effect of nadir-only sampling is also evident by the more similar correlations of unscaled C3M and SYN near the poles, since the polar orbits of *CloudSat*/*CALIPSO* sample these latitudes more frequently. However, even after correcting for both sampling effects, the correlations remain smaller than those obtained using the SYN clouds.

Instead of scaling using SSF-based calculations, the additional noise from the nadir-only sampling could be reduced by applying more averaging to the C3M data. However, while increasing the spatial resolution by a factor of 10 [i.e., using a C3M grid size of 10° × 10° (latitude × longitude)] did improve correlations, the increased averaging only gives about half the improvement as scaling by the SSF-based calculations (not shown). Additionally, we also found similar results to those in Fig. 7 when using the raw C3M cloud groups, that is, without the clustering method applied, indicating that our clustering method has little impact on modifying the variability (not shown).

## 6. Perturbation calculations

The methodology outlined in section 4 is now used to perform calculations to demonstrate the application of CERES-PRP to decomposing the variability of the radiation budget and in deriving/validating radiative kernels. The calculations for the former application are also used to test the linearity of the perturbations.

The analysis in this section is informed by the error analysis in the previous section. For decomposing the monthly variability of the radiation budget, using the SYN clouds rather than C3M clouds is preferred as suggested by the correlations to EBAF-TOA (Fig. 6). Both sets of meteorological properties (GEOS or AIRS) produce fluxes with similar monthly variability as EBAF-TOA. For the purpose of computing radiative kernels, the error in the perturbation itself is less of a concern since the computed flux perturbation is normalized by . More important for radiative kernels is the accuracy of the variables used in the flux calculations. As previously discussed, C3M clouds (AIRS meteorological properties) give fluxes with smaller CRE (clear-sky flux) biases (Fig. 5). Therefore, we expect the combination of C3M + AIRS to provide the most accurate radiative kernels.

### a. Linearity

A CERES-PRP calculation is performed where each variable in Table 1 is perturbed individually for SYN + GEOS during the *Terra* period. The global mean flux anomalies are shown in Figs. 8 and 9 for the all sky and clear sky, respectively. In these figures, for simplicity, the anomalies from the individual cloud properties have been summed together to show the total effect of clouds. The total effect of all trace gases (ozone, carbon dioxide, methane, nitrous oxide, CFC-11, CFC-12, and HCFC-22) along with the incoming solar irradiance are summed as well. In the shortwave clear sky (Fig. 9a), the variability is driven mostly by a combination of aerosols and the surface albedo, with the latter causing a decreasing trend in the reflected shortwave. The longwave clear sky (Fig. 9b) is mostly an interplay between temperature and water vapor with a strong anticorrelation between the two. Also apparent is the steady, albeit small, decrease in the outgoing longwave radiation from greenhouse gases. For the all-sky shortwave fluxes (Fig. 8a), clouds dominate the variability. The all-sky longwave anomalies (Fig. 8b) are composed of a more complicated balance of contributions from temperature, water vapor, and cloud anomalies. Insights from figures like these can aid in understanding the underlying processes associated with variations in Earth’s radiation budget. However, an in-depth analysis of the features in Figs. 8 and 9 is beyond the scope of this current study. A more detailed application of similar CERES-PRP calculations is presented in Loeb et al. (2018b), who used the CERES-PRP calculations to investigate recent changes in Earth’s energy budget along with changes during and after the global warming “hiatus.”

The calculations in Figs. 8 and 9 are used to test the linearity of the CERES-PRP perturbations. For CERES-PRP perturbation calculations to be useful, interactions between individual perturbations must be small. Therefore, we expect the monthly flux anomalies computed in the usual fashion, that is, , to be approximately equal to the sum of the individual perturbations, , where *i* are the individual perturbed variables. Figure 10 shows this nonlinear error, that is, the difference between and for SYN + GEOS and the *Terra* period (similar results are obtained if the *Aqua* or C3M periods are used). The boxplots in Fig. 10 summarize the relative error found across all grid boxes and all months in the time series. The error is narrowly distributed and generally quite small with medians near zero and interquartile ranges within about ±10%. The 95% confidence intervals do show there are instances of more substantial errors, but overall the CERES-PRP calculations provide an approximately linear decomposition of the total flux anomalies.

A similar analysis was performed for C3M + GEOS in Fig. 11. Like the SYN + GEOS assessment, the median C3M + GEOS errors are less than a few percent. However, the interquartile ranges are larger, in particular for the all-sky fluxes. This is likely due to the increased complexity in organizing the C3M cloud conditions. Using AIRS instead of GEOS with either the SYN or C3M cloud properties produced similar nonlinear errors (not shown).

### b. Water vapor radiative kernel

CERES-PRP calculations were used to validate the NCAR CAM3.0 water vapor kernel (Shell et al. 2008). While other GCM-derived water vapor kernels exist, we chose to focus on the NCAR kernel, as this set of kernels has been shown to better reproduce TOA flux changes (Vial et al. 2013). Our analysis mirrors that of Soden et al. (2008), who documented differences in water vapor kernels from three different GCMs. These differences among GCMs are compared to the differences between the CERES-PRP observation-based kernels and the NCAR kernel.

Water vapor kernels are derived from the four different dataset combinations used throughout our analysis. As described above, a kernel is obtained by using separate CERES-PRP calculations to get the effect on the TOA flux from each individual water vapor level [i.e., compute Eq. (9) for each level]. A final kernel is then obtained by compiling the from each level. The anomalies of water vapor are computed in log space [i.e., ], since the absorption of radiation responds approximately logarithmically to the water vapor concentration (Soden et al. 2008). The anomalies used to compute the kernel are relative to the C3M period, and kernels are not sensitive to the chosen time period (not shown). Data from 2007 are used to compute monthly mean kernels. This follows the typical practice of using a single year-long control run to compute kernels (e.g., Soden et al. 2008; Shell et al. 2008). The impact of interannual variability is examined further at the end of this section.

For CERES-PRP kernels, the anomalies are in the same units as those expected by the radiative transfer model, that is, for water vapor, the mass mixing ratio. However, GCM-based water vapor kernels use a nonconstant standard water vapor perturbation defined as the specific humidity increase corresponding to a 1-K temperature perturbation at a constant relative humidity, giving the kernel units of W m^{−2} K^{−1} (Soden et al. 2008). To ease comparisons, the CERES-PRP kernels are converted to match the standard GCM convention by first converting the water vapor mass mixing ratio perturbations to specific humidity perturbations. This specific humidity kernel is then multiplied by the change in the logarithm of the specific humidity increase corresponding to a 1-K temperature perturbation at a constant relative humidity for the given month. This gives a CERES-PRP kernel for the same standard perturbation used in the GCM-based kernels.

As discussed in section 4, we assume that the observed anomalies used for the CERES-PRP kernels are small enough to give acceptably small truncation errors. For these water vapor kernel calculations, we find that, on average, the observed specific humidity anomalies are slightly smaller than those given by applying the standard GCM perturbation to our monthly observations (not shown). Therefore, we expect truncation errors in the flux partial derivatives to also be slightly smaller as well. This error is also further minimized by the use of the centered difference in the CERES-PRP calculations.

Figures 12 and 13 show the annual-mean zonal-mean CERES-PRP kernels for the TOA net (downward minus upward) shortwave and longwave fluxes, respectively. Figure 14 shows the kernels integrated vertically over the troposphere along with global mean values of the kernels. The tropopause is defined as 100 mb at the equator, decreasing linearly with latitude to 300 mb at the poles. Global means are computed only over the latitudes sampled by the C3M data (i.e., excluding latitudes poleward of 81°N or 81°S). The four dataset combinations produce similar kernels in many respects. The response to water vapor perturbations is much larger in the longwave than the shortwave by about a factor of 4–5. The shortwave kernel is most sensitive in regions where increased reflection enhances the water vapor absorption, that is, over higher surface albedos near the poles and regions of low clouds. The longwave fluxes are most sensitive to water vapor perturbations in the tropical upper troposphere. Compared to using GEOS (SYN) properties, AIRS (C3M) properties produce slightly less sensitive kernels. The reduced sensitivity from the C3M cloud properties is expected, since its larger amount of cloud, particularly high cloud, masks any underlying water vapor changes.

Figure 14 also shows the vertically integrated values of the NCAR kernel, which is more (less) sensitive in the longwave (shortwave) than any of the CERES-PRP kernels. More detailed annual-mean zonal-mean comparisons between the CERES-PRP C3M + AIRS and the NCAR kernel are given in Fig. 15. The C3M + AIRS kernel is chosen as the reference since this combination of datasets produced the smallest errors when compared to EBAF-TOA, as discussed earlier. Figure 15 shows that while the relative errors can be quite large in certain regions, this typically only occurs where the absolute values of the kernel are small. Overall, the structures of the kernels are similar. In the upper troposphere, the observation-based C3M + AIRS longwave kernel is less (more) sensitive in the tropics (extratropics). Note that it is possible for differences to be caused by discrepancies in the radiative transfer models themselves (Soden et al. 2008), something that is not explored in this work.

Differences between the sum of the shortwave and longwave NCAR and CERES-PRP C3M + AIRS kernels are within about 10%–40% for annual-mean zonal-mean values, 13% for vertical-integrated annual-mean zonal-mean values, and 8% for global means. This reduction of errors when averaging demonstrates that errors in the NCAR kernel are more of a concern if one, for example, wanted to examine local water vapor feedbacks rather than simply the global mean feedback. In contrast, Soden et al. (2008) reported that GCM water vapor kernels typically agree to within 10%–20% for annual-mean zonal-mean values, 10% for vertical-integrated annual-mean zonal-mean values, and 5% for global means. Therefore, our observation-based validation gives larger errors than has been inferred from different GCMs. However, these differences are still somewhat similar to past assessments, especially for global mean values.

Kernels were constructed using data from 2007 for the above comparisons as it is typical to use a single year to derive kernels (e.g., Soden et al. 2008; Shell et al. 2008). While previous studies have recognized that ignoring the effects of interannual variability may be a source of error, none have quantified its impact. To explore this further, the SYN + GEOS observations from 2005 to 2009 were used to compute water vapor kernels. Over this time period, the annual-mean global mean values of the longwave (shortwave) kernel varies from 1.154 to 1.176 (0.276 to 0.281) W m^{−2} K^{−1}. None of the global mean differences between pairs of years are statistically significant. However, when each year is compared to the NCAR water vapor kernel, global mean values of each year are significantly different. Therefore, while there does exist some effect of interannual variability on the kernels, we do not expect this to explain a significant portion of the differences between the CERES-PRP and NCAR kernels.

## 7. Summary and conclusions

The CERES-PRP methodology applies partial radiative perturbation (PRP)-like calculations (i.e., finite differences) to observational datasets to directly isolate the individual contributions to the variability of the radiation budget. Calculations are performed in a consistent framework to accommodate the calculation of flux perturbations due to any desired combination/subset of variables. This flexibility allows the CERES-PRP methodology to generate datasets at will and tailored to the specific question being investigated. The results of these calculations can further be used to construct radiative kernels.

For the sake of computational efficiency, monthly mean inputs are used for the radiative flux calculations. The nonlinear errors from these time-averaged inputs are reduced by including a monthly mean representation of the diurnal cycle and using a gamma-weighted 2-stream shortwave solver to better capture the effective distribution of cloud optical depths within a particular month. A *k*-means-based clustering method was applied to cloud properties from the *CloudSat* cloud radar and *Cloud–Aerosol Lidar and Infrared Pathfinder Satellite Observations* (*CALIPSO*) lidar to reorganize the profiles into monthly means suitable for accurate radiative transfer calculations. Errors from using monthly mean inputs were assessed by comparing to the averages of fluxes computed at their native temporal resolutions. Overall, we find the use of monthly mean inputs reasonably accurate, reproducing the true TOA fluxes to within 6%.

Uncertainties in the computed fluxes were explored through comparisons to the CERES EBAF-TOA observed fluxes. Using the Atmospheric Infrared Sounder (AIRS) meteorological properties (i.e., water vapor, temperature, and ozone) instead of those from the Goddard Earth Observing System (GEOS) reanalysis improves agreement in the clear-sky longwave fluxes by about a factor of 2. Using the C3M cloud properties produces fluxes that agreed more favorably to the EBAF-TOA cloud radiative effect (CRE) than the CERES synoptic (SYN) MODIS and geostationary-based cloud properties. Comparisons of the EBAF-TOA zonal-mean time series showed that fluxes computed using the SYN cloud properties agreed well with correlation coefficients above 0.90 for most latitudes. Using C3M cloud properties gives considerably poorer correlations due to *CloudSat*/*CALIPSO*’s limitation to nadir-only sampling.

The CERES-PRP methodology was demonstrated by decomposing the variability of the radiation budget and by deriving/validating water vapor radiative kernels. The calculations for the former were used to test the linearity of the CERES-PRP perturbation calculations. This nonlinear residual was found, on average, to be very close to zero. Errors are slightly larger when C3M cloud properties are used, likely because of the increased complexity in organizing the C3M cloud conditions. CERES-PRP water vapor kernels were used to validate the widely used NCAR CAM3.0 water vapor kernel (Shell et al. 2008). This observation-based validation showed that errors in the NCAR kernel are larger than those given in Soden et al. (2008), who documented differences in water vapor kernels among three different GCMs. Overall, these errors are still somewhat similar to past assessments, especially for global mean values.

Although we mostly focused on the development and validation of the methodology in this study, the decomposition of the variability of the radiation budget provided by CERES-PRP calculations can be a useful aid in understanding the observed radiation budget. Loeb et al. (2018b) used the CERES-PRP calculations to investigate recent changes in Earth’s energy budget along with changes during and after the global warming “hiatus.” CERES-PRP also provides a means for examining feedbacks and other radiative relationships (e.g., radiative kernels) in the observational record. Radiative kernels (temperature, water vapor, cloud fraction, and cloud-base kernels) computed by CERES-PRP are used edition 4 of the Energy Balanced and Filled (EBAF)-Surface product (Kato et al. 2018) to assist in adjusting computed surface fluxes for various bias corrections. While this study only focused on the water vapor kernel, work is underway to make comparisons for the full set of radiative kernels and determine how observation-based kernels impact estimates of GCM feedbacks.

## Acknowledgments

This research has been supported by the NASA CERES and *CALIPSO* projects. Part of this work was performed while T. J. Thorsen was supported by a NASA Postdoctoral Program Fellowship. The CERES datasets (EBAF-TOA, SYN, SSF, C3M, SAH) are available online (https://ceres.larc.nasa.gov). AIRS data were obtained from the GES DACC (https://disc.gsfc.nasa.gov). The NOAA-ESRL data are available online (https://www.esrl.noaa.gov/gmd/ccgg). The SORCE TSI data were obtained online (http://lasp.colorado.edu/home/sorce/data/tsi-data). The NASA Langley Fu–Liou radiative transfer model can be downloaded (https://www-cave.larc.nasa.gov/cgi-bin/lflcode/accesslfl.cgi). Those interested in kernels or other applications of CERES-PRP are encouraged to contact the corresponding author.

## REFERENCES

*K*-means++: The advantages of careful seeding.

*Proc. 18th Annual ACM-SIAM Symp. on Discrete Algorithms*, Philadelphia, PA, Society for Industrial and Applied Mathematics, 1027–1035.

_{2}

_{2}

_{4}burden

_{2}forcing. Part I: Adapting the linear radiative kernel technique to feedback calculations for a broad range of forcings

**19**, 6263, https://doi.org/10.1175/JCLI9028.1.

## Footnotes

For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).