# Challenges in constraining anthropogenic aerosol effects on cloud radiative forcing using present-day spatiotemporal variability

^{a}Atmospheric Sciences and Global Change Division, Pacific Northwest National Laboratory, Richland, WA 99352;^{b}Institute for Climate and Global Change Research, Nanjing University, 210023 Nanjing, China;^{c}School of Atmospheric Sciences, Nanjing University, 210023 Nanjing, China;^{d}Collaborative Innovation Center of Climate Change, 210023 Nanjing, China;^{e}Institute for Atmospheric and Climate Science, ETH Zurich, 8092 Zurich, Switzerland;^{f}National Center for Atmospheric Research, Boulder, CO 80305;^{g}Information Technology Division, Norwegian Meteorological Institute, 0313 Oslo, Norway;^{h}Atmospheric, Oceanic and Planetary Physics, Department of Physics, University of Oxford, Oxford OX13PU, United Kingdom;^{i}Department of Environmental Science and Analytical Chemistry, Stockholm University, SE-106 91 Stockholm, Sweden;^{j}Bert Bolin Centre for Climate Research, Stockholm University, SE-106 91 Stockholm, Sweden;^{k}Research Institute for Applied Mechanics, Kyushu University, Fukuoka 816-8580, Japan

See allHide authors and affiliations

Edited by John H. Seinfeld, California Institute of Technology, Pasadena, CA, and approved January 21, 2016 (received for review September 26, 2015)

### This article has a Correction. Please see:

## Abstract

A large number of processes are involved in the chain from emissions of aerosol precursor gases and primary particles to impacts on cloud radiative forcing. Those processes are manifest in a number of relationships that can be expressed as factors dlnX/dlnY driving aerosol effects on cloud radiative forcing. These factors include the relationships between cloud condensation nuclei (CCN) concentration and emissions, droplet number and CCN concentration, cloud fraction and droplet number, cloud optical depth and droplet number, and cloud radiative forcing and cloud optical depth. The relationship between cloud optical depth and droplet number can be further decomposed into the sum of two terms involving the relationship of droplet effective radius and cloud liquid water path with droplet number. These relationships can be constrained using observations of recent spatial and temporal variability of these quantities. However, we are most interested in the radiative forcing since the preindustrial era. Because few relevant measurements are available from that era, relationships from recent variability have been assumed to be applicable to the preindustrial to present-day change. Our analysis of Aerosol Comparisons between Observations and Models (AeroCom) model simulations suggests that estimates of relationships from recent variability are poor constraints on relationships from anthropogenic change for some terms, with even the sign of some relationships differing in many regions. Proxies connecting recent spatial/temporal variability to anthropogenic change, or sustained measurements in regions where emissions have changed, are needed to constrain estimates of anthropogenic aerosol impacts on cloud radiative forcing.

Radiative forcing of climate change through interactions between liquid clouds and anthropogenic aerosol arises through a chain of processes from emissions of primary particles and aerosol precursor gases *E*, to establishment of a balance between production and removal of cloud condensation nuclei (*CCN*), to effects of the *CCN* on droplet number concentration *N*_{d}, to effects of *N*_{d} on cloud radiative forcing *R*.

This chain can be expressed mathematically for a single-layer liquid cloud*R* is the “clean-sky” shortwave cloud radiative forcing, i.e., the shortwave cloud radiative forcing calculated as a diagnostic with aerosol optical depth set to zero (1). Note that this formalism allows feedbacks such as cloud effects on *CCN*, so the terms should not be interpreted as only the response of the numerator to changes in the denominator.

Cloud radiative forcing can be expressed as the product of cloud fraction *C* and the clean-sky cloud radiative forcing for the cloudy fraction of the sky, *R*_{c}, so the first term on the right-hand side (RHS) of Eq. **1** can be expressed as*R*_{c} depends almost entirely on the cloud optical depth *τ*, the second term on the RHS of Eq. **2** becomes*τ* and *N*_{d} can be decomposed into contributions from changes in droplet effective radius *r*_{e} and cloud liquid water path *L* using the common expression for cloud optical depth *τ* ∝ *L*/*r*_{e} (2),**1**−**4** gives**5** is often called the first indirect, Twomey (3), or cloud albedo effect, and the first and second terms are together called the second indirect, Albrecht (4), or cloud lifetime effect. The cloud lifetime effect was originally associated with changes in cloud fraction, but, as can be seen from Eq. **5**, it involves changes in both cloud fraction and liquid water path. There is, however, an important distinction between Eq. **5** and the Twomey mechanism, which assumes constant *L*: Eq. **5** makes no such assumption.

In model estimates, these influences are determined from the difference between pairs of simulations [say, preindustrial (PI) and present day (PD)] by each atmosphere model, with the only difference in configuration being the anthropogenic emissions,*R* (about 1 W⋅m^{−2}) is much smaller than *R* (about 50 W⋅m^{−2} in the global mean) and the other numerators and denominators in Eq. **6** cancel, it is essentially a statement of an identity.

Uncertainty in estimates of the effective radiative forcing through aerosol−cloud interactions (ERFaci) can arise through uncertainty in each of the terms in Eq. **6**. Although one might argue that it is only the final uncertainty that matters, efforts to reduce the uncertainty in ERFaci can be most effectively focused if the uncertainty in each term is known. Quantification of the uncertainty from each term has not been attempted before.

Uncertainty can arise from uncertainty in the value of model parameters (parametric uncertainty), in the limitations of the formulations of the physical processes represented in the model (structural uncertainty), and in the numerical representation of the physical processes (numerical uncertainty). Parametric uncertainty can be quantified by simultaneously varying the values of uncertain model parameters within the range of their uncertainty (5, 6). Structural uncertainty is commonly estimated by comparing the values of each term from different models against observations. Additional uncertainty is evident through lack of agreement between simulated and observed estimates of the terms.

Uncertainty in each term can be reduced if observations can be used effectively to constrain the values. However, if radiative forcing is estimated over a period predating reliable measurements (e.g., beginning with preindustrial conditions), the necessary observations are not available. Although spatial and temporal variations in terms over recent periods (such as the satellite era) have been used to estimate several of the terms in Eqs. **6**−**8**, Penner et al. (7) showed that, at least for the relationship between *N*_{d} and aerosol optical depth (a proxy for *CCN*) from one model, the relationship estimated from recent variations is not a useful constraint on the preindustrial to present-day relationship. Stier (8) has recently shown that aerosol optical depth (AOD) is not always a good proxy for *CCN*, but we avoid this issue by using *CCN* rather than AOD.

This study addresses several related issues. First, following Schulz et al. (9), we estimate the structural uncertainty in each of the terms in Eqs. **6**−**8** using differences between simulations by a suite of nine Aerosol Comparisons between Observations and Models (AeroCom) atmosphere models with present-day and preindustrial emissions to determine which terms contribute most to the diversity of the estimated radiative forcing. Second, for each model, we determine how well the terms estimated from present-day variations match the terms determined from differences between results for PI and PD emissions. Third, we discuss alternate methods of constraining the terms.

## Structural Uncertainty for Each Term

Each term contributes to the structural uncertainty in the aerosol radiative forcing. Fig. 1*A* shows normalized values of each term in Eq. **6**, except emissions, for all nine models, after averaging the numerators and denominators globally over five simulated years (see *Methods*). Because the emissions of primary and precursor anthropogenic aerosol mass are the same for all models in this study, we show the ln*CCN* change rather than its relationship with the emissions change. Emissions are also a source of uncertainty (10), but that source is not considered in this study.

The threefold range in the global mean radiative forcing is driven by diversity in all factors.

Because satellite observations provide a strong constraint on the present-day cloud radiative forcing, and because the difference between the “dirty sky” and clean-sky cloud radiative forcing is less than 1 W⋅m^{−2} (much smaller than the present-day shortwave cloud radiative forcing) (1), one might expect little diversity in the simulated cloud radiative forcing, the first term in Eq. **6**. Indeed, there is little diversity when all clouds are sampled. However, according to Fig. 1*A*, there is considerable diversity when only warm clouds are sampled (normalized values ranging between 0.7 and 1.3), because the climate model calibration process focuses on all clouds rather than on warm clouds only.

The relationship between *R* and *N*_{d} has somewhat more diversity, with normalized values between 0.7 for SPRINTARS (Spectral Radiation Transport Model for Aerosol Species) and HadGEM3-U.K.CA (Hadley Center Global Environmental Model with United Kingdom Chemistry and Aerosols) and 1.4 for CAM5.3_CLUBB (Community Atmosphere Model 5.3 with CLouds Unified By Binormals). This diversity is driven by diversity in several terms. Fig. 1*B* shows the same relationship, not normalized by the multimodel mean, and, following Eq. **8**, contributions to that relationship from the relationships of low cloud fraction and from in-cloud radiative forcing *R*_{c} with *N*_{d}. The large value for the *R*−*N*_{d} relationship for CAM5.3_CLUBB is evidently explained by the much larger contribution of cloud fraction changes to the cloud radiative forcing response for CAM5.3_CLUBB. This can be understood by noting that CAM5.3_CLUBB treats aerosol effects on shallow cumulus clouds (11), whereas most other models do not. However, it is not clear why CAM5.3_CLUBB_MG2 (CAM5.3 with CLUBB and second generation Morrison & Gettelman cloud microphysics), which also treats aerosol effects on shallow cumulus clouds, does not produce a large contribution of cloud fraction changes to the cloud radiative forcing response.

For almost all other models, Fig. 1*B* shows that the relationship of the in-cloud radiative forcing *R*_{c} to *N*_{d} exceeds the relationship of cloud fraction to *N*_{d}, but the diversity across models is comparable.

The relationship of *R*_{c} to *N*_{d} is particularly small for CAM5.3_CLUBB. This is surprising, because Fig. 1*C* shows that the relationship of *τ* to *N*_{d} is relatively strong for CAM5.3_CLUBB. That implies, according to Eq. **7**, that the relationship of *R*_{c} to *τ* is much smaller for CAM5.3_CLUBB. Indeed, it is an order of magnitude smaller for CAM5.3_CLUBB than for all other models (see Dataset S1). This unexpected result is difficult to explain. Potential explanations are the sublinear dependence of cloud albedo on cloud optical depth, a tendency for cloud changes to occur over snow or ice, or the dominant influence of variations in solar zenith angle on cloud radiative forcing. However, the mean cloud optical depth is no larger for CAM5.3_CLUBB than for other models, the spatial distribution of the cloud optical depth change for CAM5.3_CLUBB is not predominantly over bright surfaces, and weighting cloud fraction and cloud optical depth by the incoming solar flux does not affect the relationship between *R*_{c} and *τ* for CAM5.3_CLUBB.

Fig. 1*C* shows the relationship between *τ* and *N*_{d}. The relationship differs considerably across the nine models, with a range of a factor of 2, because the relationship depends on cloud microphysical processes (12), which are treated differently in different models. The low values of the relationship for the SPRINTARS and HadGEM3-U.K.CA models clearly contribute to their low estimates of the relationship between *R* and *N*_{d}, although the cloud fraction relationship to *N*_{d} is also small. The KK version of the SPRINTARS model produces a much larger value of the *τ*−*N*_{d} relationship, as might be expected because it uses the Khairoutdinov and Kogan (13) autoconversion scheme used by all versions of CAM5.3, which also produce larger values. To confirm this, we use Eq. **8** to decompose the relationship into contributions from changes in *L* and *r*_{e}. The diversity in these two terms from global and annual mean changes in the numerators and denominators is also shown in Fig. 1*C*. Note that, because these terms are added rather than multiplied, their values have not been normalized by the means across the models. However, the two terms don’t strictly add because the finite differences in Eq. **8** approximate the differential form.

From Fig. 1*C*, it is obvious that the diversity in the relationship between *τ* and *N*_{d} is dominated by the difference between the SPRINTARS and HadGEM3-U.K.CA models and the other models, and that the difference in the *τ*−*N*_{d} relationship between the SPRINTARS and HadGEM3-U.K.CA models and the other models is dominated by the difference in the relationship between *L* and *N*_{d}. Because the relationship between *L* and *N*_{d} for SPRINTARSKK (SPRINTARS with Khairoutdinov and Kogan autoconversion scheme) is about halfway between the relationship for SPRINTARS and the rest of the models, we can conclude that, to a significant extent, the much weaker relationship from SPRINTARS and HadGEM3-U.K.CA is due to their use of a different autoconversion schemes. However, because the liquid water path relationship from SPRINTARSKK is about half of the relationship of the rest of the models that use the same autoconversion scheme, there must be some other model differences that also contribute to the weaker relationship in the SPRINTARS and HadGEM3-U.K.CA models.

Note that the relationship between *L* and *N*_{d} from CAM5.3_MG2 is nearly indistinguishable from that simulated by the other versions of CAM5.3. This result is surprising because previous work (11, 14, 15) suggested that the prognostic treatment of rain in CAM5.3_MG2 would produce a weaker relationship between *L* and *N*_{d} than the diagnostic treatment in the other models.

It is also noteworthy that, for all models except SPRINTARS, SPRINTARSKK, and HadGEM3-U.K.CA, the relationship between *r*_{e} and *N*_{d} is weaker than the relationship between *L* and *N*_{d}. The metrics used in Fig. 1*C* provide a convenient method of quantitatively comparing the magnitudes of these two different mechanisms producing aerosol effects on cloud optical depth.

Note also that the relationship between *r*_{e} and *N*_{d} is substantially weaker for ECHAM6-HAM (Hamburg version of European Center for Medium-range Weather Forecasting model) than the other models, mostly because *r*_{e} changes for ECHAM6-HAM are smaller.

Returning to Fig. 1*A*, the relationship between *N*_{d} and *CCN* varies widely (more than twofold) across the models. The CAM5.3 versions produce fairly consistent relationships, but the two SPRINTARS versions produce much weaker relationships, which contribute substantially to the smaller estimates of radiative forcing. Although all models except ECHAM6-HAM and HadGEM3-U.K.CA use the same droplet nucleation scheme, SPRINTARS and ECHAM6-HAM apply a lower bound on *N*_{d} (16, 17) that clearly substantially limits the droplet number sensitivity to *CCN* for clean conditions such as the preindustrial (18). Differences in updraft velocity and in natural emissions could also contribute to the diversity (19), but additional experiments without lower bounds would be required to determine the contribution of the lower bound to the diversity.

The response of boundary layer *CCN* concentration to anthropogenic emissions is surprisingly similar across all models except HadGEM3-U.K.CA, with agreement to within 20% for the other models. The agreement would likely be worse for supersaturations higher than the chosen supersaturation (0.3%), but supersaturations are higher than 0.3% only for strong updrafts or very clean conditions (20). Some of the agreement is undoubtedly due to the use of common treatments of aerosol emissions, properties, and processes in the CAM5.3 versions, with only CAM5.3_PNNL (Pacific Northwest National Laboratory version of CAM5.3) using a different representation of cloud effects on the aerosol. However, the representations of the aerosol and its lifecycle in SPRINTARS and ECHAM6-HAM are quite different, with external mixtures of components assumed in SPRINTARS and internally mixed modes in ECHAM6-HAM. The 60% larger anthropogenic *CCN* change for HadGEM3-U.K.CA, which also uses an internally mixed modal treatment, could be due to differences in size distributions of primary anthropogenic emissions (21⇓–23) (which affects the number emitted), in the treatments of particle nucleation (24), or in the representations of clouds, which affects vertical transport, wet removal, and aerosol lifetime (25). However, although the *CCN* change for HadGEM3-U.K.CA is by far the largest, the aerosol radiative forcing by HadGEM3-U.K.CA is among the smallest, mostly because of its small sensitivity of cloud radiative forcing to changes in *N*_{d}.

*CCN* measurements are sparse and seldom used to constrain aerosol models, and previous studies have found considerable dependence of *CCN* concentration on uncertainty in model parameters (5, 26) and processes (nucleation, wet removal) (24, 27, 28), so much greater diversity was expected. This shows that model diversity is not always a reliable measure of structural uncertainty (29). Comparison with more models using different treatments of the aerosol lifecycle is needed (24).

## Relationships from PD Variations vs. PI−PD Changes

Given the diversity of the relationships, constraints are needed to guide model development and reduce uncertainty in estimates of the radiative forcing. Unfortunately, the preindustrial observations needed to constrain the sensitivities are not available. One could estimate relationships from regressions on spatial and temporal variability during the recent decades when data from satellite, aircraft, and surface-based remote sensing are available (30⇓⇓⇓⇓⇓⇓⇓⇓⇓⇓–41). However, Penner et al. (7) have shown that, at least for the relationship between *N*_{d} and aerosol optical depth from one model, estimates from recent variability underestimate the relationship from anthropogenic emissions.

We have explored this issue for all relationships estimated from global mean anthropogenic change and from global spatial and temporal variability. We have found that these two methods of estimating the relationships differ considerably, differing even in sign for some relationships and models. However, Grandey and Stier (42) have found 20% biases in estimates of the relationships for regions larger than 15° and 50% errors for regions larger than 30°. We therefore focus on the relationships at subglobal scales, specifically for the 14 regions defined by Quaas et al. (43) and used by Penner et al. (7). Rather than use aerosol optical depth as a proxy for *CCN*, we use *CCN* concentration to avoid ambiguity from the use of a proxy.

Fig. 2 explores this issue for all relationships using all models in this study. Relationships estimated from the anthropogenic change (PD−PI) are compared with estimates from regressions over spatial (within each region) and temporal (3 h to 5 y) variability of the 3-h PD simulation data.

For the relationship between *N*_{d} and *CCN* (Fig. 2*A*), the two estimate methods yield similar results for some regions for all models, but poor agreement for other regions for almost all models. Six of the models yield a value less than 0.1 for the relationship estimated from PD variability for one region, but values exceeding 0.6 when estimated from the PI to PD change for that region. The difference in the estimates could reflect saturation in droplet number for PD conditions. Only CAM5.3-CLUBB-MG2 and HadGEM3-U.K.CA yield agreement within a factor of 2 for all regions. Most models have much less interregional diversity in the relationship derived from variability than from anthropogenic change, which suggests estimates of the relationship from variability are poor constraints on the relationships from anthropogenic change. Penner et al. (7) also found spatial variability in the bias, but the estimate from variability mostly underestimated the value from anthropogenic change.

For the relationship between *τ* and *N*_{d} (Fig. 2*B*), the estimate from variability is less than the estimate from anthropogenic change for all models and nearly all regions. Even the sign of the relationship differs for at least one of the regions for most models. Note that negative slopes for PD variability have also been found using cloud-resolving model simulations (44), which suggests that those simulations are not necessarily effective constraints on the *τ* response to anthropogenic aerosol.

For *L* and *N*_{d} (Fig. 2*D*) the results are similar to those for *τ* and *N*_{d}, but with more regions with more negative PD slopes. The sign of the relationship differs for about half of the regions and for at least one region for all but one model.

The two methods yield much more consistent estimates of the relationship between *r*_{e} and *N*_{d} (Fig. 2*B*) for all models except SPRINTARS and SPRINTARSKK, which yield opposite signs for several regions. For most models and most regions the estimates from present day spatial/temporal variability can be used to constrain the anthropogenic change.

For *R* and *τ* (Fig. 2*E*), the signs agree for nearly all regions and models, but, for many regions and models, the relationship estimated from variability exceeds that from anthropogenic change by more than a factor of 2; spatial variations in solar zenith angle within those regions could explain some of this difference. For six of the nine models, the relationship from variability in at least one region is less than 0.05, whereas the estimate from anthropogenic change is larger than 0.25.

For cloud fraction and *N*_{d} (Fig. 2*F*), the signs of the relationship from variability and anthropogenic change differ for most regions and models, with a positive value from anthropogenic change and a negative value for variability from many regions, and vice versa in other regions. Processes other than aerosol activation, such as turbulence or large-scale motion, likely play a relatively larger role when only spatial/temporal variability is driving the relationship.

## Emerging Constraints on Estimates of Past Forcing

Given the poor agreement between the two methods of estimating most relationships for most models, other ways of using recent measurements based on spatial/temporal variability are needed to constrain the anthropogenic influence of aerosol on clouds and the Earth’s energy balance.

One method would be to use recent trends in regions where emissions have changed substantially during the period when reliable measurements are available. For example, Cherian et al. (45) used measurements of trends in the downward solar radiance at European sites from the period 1990–2005, when SO_{2} emissions declined threefold (46), to constrain global estimates of aerosol radiative forcing since the preindustrial era. Although such an analysis is highly informative, it does not provide guidance on removing biases in models that overestimate or underestimate the downward solar trend over Europe, which could be due to errors in any of the factors that produce the cloud radiative forcing change or the clear-sky change, as well as natural variability in cloud cover. Removing those biases is necessary if climate models are to be used for simulations of future climate change. Additional data characterizing each of the factors and components are needed. Some of the necessary data (*L*, *r*_{e}, aerosol optical depth) are available from 1990, but reliable estimates of *N*_{d}, *τ*, and *R* are not available for years before 2001, when the Earth Observing System satellite constellation was launched.

Although most of the reduction in emissions from Europe had already occurred by 2001, emissions from east Asia continued to rise through 2007 (47). This presents an opportunity to constrain the factors and components that contribute to aerosol radiative forcing, if the aerosol signature exceeds radiance changes due to natural variability in clouds.

A third opportunity is the large tropospheric emission of SO_{2} from the Bárðarbunga volcano on Iceland between 29 August 2014 and 27 February 2015. A group led by James Haywood is studying this promising case.

Another approach is to develop metrics from variability that can constrain the anthropogenic sensitivity of selected factors and components. For example, Wang et al. (48) have shown that S_{pop} ≡ −dln(POP)/dlnAI, a measure of the sensitivity of the probability of precipitation (POP) to aerosol [expressed as Aerosol Index (AI), the product of the aerosol optical depth and the Angstrom exponent], is highly correlated with the simulated anthropogenic change in *L* to *CCN*. Because S_{pop} can be determined from recent measurements, both satellite (48) and ground-based (49), the S_{pop} measurements can be used to constrain the relationship between anthropogenic changes in *L* and *CCN*.

Ref. 48 considered only three fundamentally different models, and then adjusted parameters in one model to produce a wider range in results. Fig. 3 applies the ref. 48 analysis to the nine aerosol models in this study. As in ref. 48, the two variables are well correlated, but land values (not considered in ref. 48) are less well correlated. Interestingly, the SPRINTARS model produces values of S_{pop} closest to the ocean mean value (0.12) estimated from satellite by ref. 48. The strongly negative values of S_{pop} produced by HadGEM3-U.K.CA are suggestive of aerosol invigoration of convective clouds, but HadGEM3-U.K.CA neglects that influence, and the focus of this analysis on low warm clouds should preclude that mechanism.

A third approach is to separate clouds further into different cloud regimes, e.g., sorting by large-scale vertical velocity and lower-troposphere static stability (50, 51). Such stratification might improve consistency between estimates of relationships from anthropogenic change and temporal variability.

## Conclusions

We have found that uncertainty in anthropogenic aerosol effects on cloud radiative forcing arises from uncertainty in several relationships, and that estimates of several of those relationships from recent spatial and temporal variability are not necessarily relevant constraints on the relationship from anthropogenic change. Because few measures of preindustrial aerosol are available, this manifestation of equifinality (52) presents considerable challenges for constraining estimates of anthropogenic aerosol radiative forcing. Constraining *R _{c}* using the observed present-day relationship between

*N*and AOD and the strong correlation between

_{d}*R*and the relationship between

_{c}*N*and AOD (43) might not be justified.

_{d}Fortunately, the S_{pop} parameter, which can be characterized from recent measurements, correlates well with the anthropogenic relationship between cloud liquid water path and aerosol. However, it is disconcerting that SPRINTARS, the model producing an S_{pop} most consistent with the estimate from satellite retrievals, apparently does so because of an unrealistically large lower bound on *N*_{d}, and that HadGEM3-U.K.CA produces negative values of S_{pop}. Moreover, Lebo and Feingold (53) showed that cloud regime can influence the sign of the relationship between S_{pop} and the influence of *CCN* on *L*.

Further analysis differentiated by cloud regime in regions influenced by recent trends in emissions could constrain other terms affecting ERFaci (50, 51). Such analyses could be used to constrain the selection of model parameter values that affect the estimate of ERFaci. Constraining the cloud fraction response is more challenging because the influence of natural variability on cloud fraction must be overcome.

Uncertainty in estimates of radiative forcing by cloud−aerosol interactions due to the choice of parameter values and numerical integration methods is an emerging area of investigation (5, 6, 10, 54). Further investigations of parametric sensitivity involving cloud lifetime effects and dependence on numerical integration methods are underway.

## Methods

A total of nine global aerosol models participated in this study. Salient features of each model are summarized in Table 1. Five of the nine models are versions of CAM5.3, and two are versions of SPRINTARS. CAM5-PNNL differs from CAM5 in the treatment of cloud effects on the aerosol (15).

All models were driven by the same AeroCom emissions for years 1850 and 2000 (55). All simulations were nudged toward winds analyzed by operational forecast centers; some were also nudged toward analyzed temperatures, but this was discouraged because moist convection simulated in some models is sensitive to temperature nudging (56). Nudging greatly limits natural variability in the aerosols and clouds, permitting robust estimates from simulations of 5 y (57).

To focus the analysis on liquid clouds, all variables in Eqs. **6**−**8** are sampled from 3-h history only when cloud top temperature exceeds −10 °C. This filters out both cold clouds and warm clouds obscured by cold clouds above. Droplet number is at cloud top. CCN concentration at 0.3% supersaturation is at 1 km above ground. The ratios in Eqs. **6**−**8** are determined after temporal and spatial averaging of the preindustrial to present-day change in each variable. For in-cloud values of *R*_{c} and *τ*, the grid cell mean values are averaged over time and space and then divided by the mean of the total cloud fraction. To permit comparison of the relative contribution of the diversity from each factor in Eq. **6** to the total diversity, the values of the terms for each model are normalized by the multimodel mean of the term,

where the overbars denote temporal and spatial mean and subscripts *j* and *k* denote different models. Terms in Eqs. **7** and **8** are not normalized because they add rather than multiply.

To characterize structural uncertainty, we follow a method used by Schulz et al. (9) to diagnose the factors driving uncertainty in aerosol effective radiative forcing through cloud−radiation interactions (formerly called aerosol direct radiative forcing). Structural uncertainty can be quantified as the SD *σ* across all models. If the factors *x*_{i} driving a product *y* are uncorrelated across all models, one can show that

where the overbar represents the multimodel mean and *N* is the number of factors. In practice, the factors can be negatively correlated for radiative forcing (9).

For the relationships estimated from present-day spatial and temporal variability shown in Fig. 2, regressions were formed after binning by *L* and lower tropospheric stability (LTS) to isolate aerosol effects from thermodynamics (48), and then averaging the regression over the joint pdf of *L* and LTS. The regressions are performed after binning by equally sampled bins in the denominator of each term using the 3-h model data.

## Acknowledgments

Peter Caldwell prodded our thinking on factorization, Rob Wood coaxed us to explore the role of cloud fraction, and Nicolas Bellouin provided helpful comments. Reviewer comments were also helpful. The Pacific Northwest National Laboratory (PNNL) is operated for the Department of Energy (DOE) by Battelle Memorial Institute under Contract DE-AC06-76RLO 1830. Work at PNNL was supported by the US DOE Decadal and Regional Climate Prediction using Earth System Models program and by the US DOE Earth System Modeling program. Work of M.W. and S.Z. performed at Nanjing University was supported by the One Thousand Young Talent Program, Jiangsu Province Specially-Appointed Professor Grant, and the National Natural Science Foundation of China (41575073). A portion of this research was performed using PNNL Institutional Computing resources. The ECHAM6-HAM model was developed by a consortium composed of ETH Zurich, Max Planck Institut für Meteorologie, Forschungszentrum Jülich, University of Oxford, the Finnish Meteorological Institute, and the Leibniz Institute for Tropospheric Research, and is managed by the Center for Climate Systems Modeling (C2SM) at ETH Zurich. D.N. acknowledges support by the Austrian Science Fund (J 3402-N29, Erwin Schrödinger Fellowship Abroad). C2SM at ETH Zurich is acknowledged for providing technical and scientific support. This work was also supported by a grant from the Swiss National Supercomputing Centre under Project ID s431. D.G.P. and P.S. acknowledge support from the United Kingdom (UK) Natural Environment Research Council Grant NE/I020148/1. P.S. and Z.K. acknowledge funding from the European Research Council (ERC) under the European Union’s Seventh Framework Programme (FP7/2007–2013) ERC project ACCLAIM (Grant Agreement FP7-280025). The development of modal version of the GLObal Model of Aerosol Processes (GLOMAP-mode) within Hadley Center Global Environmental Mode (HadGEM) is part of the United Kingdom Chemistry and Aerosols (UKCA) project, which is supported by both National Environmental Research Council (NERC) and the Joint Department of Energy & Climate Change/Department for Environment, Food & Rural Affairs Meteorology Office Hadley Centre Climate Programme. We acknowledge use of the Met Office and NERC MONSooN high performance computing system, a collaborative facility supplied under the Joint Weather and Climate Research Programme, a strategic partnership between the Met Office and the NERC. Simulations by SPRINTARS were executed with the supercomputer system SX-9/ACE of the National Institute for Environmental Studies, Japan. SPRINTARS is partly supported by the Environment Research and Technology Development Fund (S-12-3) of the Ministry of the Environment, Japan and Japan Society for the Promotion of Science KAKENHI Grants-in-Aid for Scientific Research 15H01728 and 15K12190. Computing resources for CAM5-MG2 simulations were provided by the Climate Simulation Laboratory at National Center for Atmospheric Research (NCAR) Computational and Information Systems Laboratory. NCAR is sponsored by the US National Science Foundation.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: Steve.Ghan{at}pnnl.gov.

Author contributions: S.G. and M.W. designed research; S.G., M.W., S.F., A.G., Z.K., U.L., H.M., D.N., D.G.P., P.S., T.T., H.W., and K.Z. performed research; M.W., S.F., A.G., J.G., Z.K., U.L., H.M., D.N., D.G.P., P.S., T.T., H.W., and K.Z. contributed new reagents/analytic tools; M.W. and S.Z. analyzed data; and S.G. wrote the paper.

This paper results from the Arthur M. Sackler Colloquium of the National Academy of Sciences, “Improving Our Fundamental Understanding of the Role of Aerosol–Cloud Interactions in the Climate System,” held June 23−24, 2015, at the Arnold and Mabel Beckman Center of the National Academies of Sciences and Engineering in Irvine, CA. The complete program and video recordings of most presentations are available on the NAS website at www.nasonline.org/Aerosol_Cloud_Interactions.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Data deposition: The model history is available at the AeroCom archive at aerocom.met.no/data.html.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1514036113/-/DCSupplemental.

## References

- ↵
- Ghan SJ

- ↵
- ↵
- ↵
- Albrecht BA

- ↵
- ↵
- ↵
- Penner JE,
- Xu L,
- Wang M

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Takemura T,
- Nozawa T,
- Emori S,
- Nakajima TY,
- Nakajima T

- ↵
- ↵
- Hoose C, et al.

- ↵
- ↵
- Ghan SJ, et al.

*J*. Adv Model Earth Syst 3(1):M10001. - ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Ghan SJ, et al.

- ↵
- ↵
- Posselt R,
- Lohmann U

- ↵
- ↵
- Wang M, et al.

- ↵
- Martin GM,
- Johnson DW,
- Spice A

- ↵
- Ramanathan V,
- Crutzen PJ,
- Kiehl JT,
- Rosenfeld D

- ↵
- ↵
- Chameides WL, et al.

- ↵
- ↵
- Quaas J,
- Boucher O,
- Breon FM

- ↵
- Quaas J,
- Boucher O,
- Bellouin N,
- Kinne S

- ↵
- Twohy CH, et al.

- ↵
- Wilcox EM,
- Roberts G,
- Ramanathan V

- ↵
- McComiskey A,
- Feingold G

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Wang M, et al.

- ↵
- Mann JA, et al.

- ↵
- Gryspeerdt E,
- Stier P

- ↵
- ↵
- Caballero R,
- Huber M

- ↵
- ↵
- Yan H, et al.

- ↵
- ↵
- ↵
- Kooperman GJ,
- Pritchard MS,
- Ghan SJ,
- Sommerville RCJ,
- Russell LM

- Gettelman A, et al.

- Kärcher B,
- Lohmann U

- Lohmann U,
- Kärcher B,
- Hendricks J

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Earth, Atmospheric, and Planetary Sciences