Separating radiative forcing by aerosol–cloud interactions and rapid cloud adjustments in the ECHAM–HAMMOZ aerosol–climate model using the method of partial radiative perturbations

Using the method of offline radiative transfer modeling within the partial radiative perturbation (PRP) approach, the effective radiative forcing by aerosol–cloud interactions (ERF_aci) in the ECHAM–HAMMOZ aerosol climate model is decomposed into a radiative forcing by anthropogenic cloud droplet number change and adjustments of the liquid water path and cloud fraction. The simulated radiative forcing by anthropogenic cloud droplet number change and liquid water path adjustment are of approximately equal magnitude at −0.52 and −0.53 W m⁻², respectively, while the cloud-fraction adjustment is somewhat weaker at −0.31 W m⁻² (constituting 38 %, 39 %, and 23 % of the total ERF_aci, respectively); geographically, all three ERF_aci components in the simulation peak over China, the subtropical eastern ocean boundaries, the northern Atlantic and Pacific oceans, Europe, and eastern North America (in order of prominence). Spatial correlations indicate that the temporal-mean liquid water path adjustment is proportional to the temporal-mean radiative forcing, while the relationship between cloud-fraction adjustment and radiative forcing is less direct. While the estimate of warm-cloud ERF_aci is relatively insensitive to the treatment of ice and mixed-phase cloud overlying warm cloud, there are indications that more restrictive treatments of ice in the column result in a low bias in the estimated magnitude of the liquid water path adjustment and a high bias in the estimated magnitude of the droplet number forcing. Since the present work is the first PRP decomposition of the aerosol effective radiative forcing into radiative forcing and rapid cloud adjustments, idealized experiments are conducted to provide evidence that the PRP results are accurate. The experiments show that using low-frequency (daily or monthly) time-averaged model output of the cloud property fields underestimates the ERF, but 3-hourly mean output is sufficiently frequent.

1 Introduction

Following Boucher et al. (2014), it has become common to distinguish between radiative forcing (RF) by aerosol–cloud interactions (ACIs) and rapid adjustments to this radiative forcing, with the sum of the forcing (RF_aci) and the rapid adjustments denoted as the effective radiative forcing (ERF_aci). In liquid-water clouds, RF_aci arises from the increased availability of cloud condensation nuclei (CCN) in a polluted atmosphere leading to higher droplet number N_d and smaller effective radius r_e at constant cloud liquid water path ℒ (Twomey, 1977). The atmosphere responds to the higher-N_d, lower-r_e clouds by various processes occurring on short timescales, which leads to adjustments in other cloud properties, including ℒ and cloud vertical and horizontal geometric extent.

Physically, the most important adjustment mechanisms are suppression of precipitation formation (Albrecht, 1989) and enhanced cloud-edge dry-air entrainment and droplet evaporation in the smaller-droplet clouds (Ackerman et al., 2004). The former mechanism, in isolation, would lead to an increase in cloud condensate and cloud fraction and thus a negative adjustment to the radiative forcing; the latter, in isolation, would lead to a decrease in cloud condensate and thus a positive adjustment to the radiative forcing. Since no mechanism occurs in isolation in a coupled system (Stevens and Feingold, 2009), the question of whether the net adjustment is positive or negative is both difficult and unresolved (e.g., Mülmenstädt and Feingold, 2018; Gryspeerdt et al., 2019a).

One method by which aerosol effective radiative forcing is estimated is to calculate the difference in top-of-atmosphere (TOA) radiative fluxes in general-circulation models (GCMs) between runs with fixed sea-surface temperature and present-day or preindustrial aerosol concentrations or emissions (e.g., Lohmann et al., 2010; Forster et al., 2016). This requires a GCM that includes the relevant aerosol–radiation and aerosol–cloud interaction mechanisms. By representing the optical properties of aerosol in the radiative transfer, the radiative forcing by aerosol–radiation interactions (RF_ari) can be estimated; by representing the N_d dependence on aerosol activation during cloud formation, RF_aci can be estimated; and by representing precipitation suppression and enhanced evaporation in smaller-r_e clouds, the adjustments to RF_aci can be estimated. (We exclude ice and mixed-phase cloud processes, which introduce further complications, from this discussion.) These processes occur on scales far below the resolved scale, so their representation in the GCM requires parameterization. Thus, the model is only imperfectly (if at all) aware of subgrid-scale variability in the process rates and feedbacks between the processes; relies on imperfect base-state cloud properties (e.g., Penner et al., 2006); and only considers effects that are amenable to parameterization, meaning that precipitation suppression is included in many models but enhanced evaporation is not (e.g., Salzmann et al., 2010; Michibata et al., 2016; Zhou and Penner, 2017). Based on these considerations, a prevalent view is that, from the standpoint of achieving GCM fidelity, ACIs are more difficult than aerosol–radiation interactions, and ACI adjustments are more difficult than the ACI forcing. On top of this, the usual concerns about parametric uncertainty apply, so that the overall uncertainty on GCM estimates of ERF_aci is large (Boucher et al., 2014). Nevertheless, the ability of GCMs to produce a global estimate will assure their star will continue to shine brightly in the firmament of ERF estimation methods until competing methods overcome their own significant drawbacks.

In general, one might argue that knowing the uncertainty on each term in a sum is a good first step towards attacking the uncertainty on the total; certainly, this is consistent with the GCM paradigm of building up the total forcing from parameterizations for each of the contributing processes, even if it is less applicable to top-down estimates from the historical evolution of the climate system. Thus, we write the effective radiative forcing by aerosol as

$$\begin{array}{}\text{(1)}& F=F_{\text{ari}}+F_{N_{\mathrm{d}}}+F_{\mathcal{L}}+F_{f_{\mathrm{c}}},\end{array}$$

where F_ari is the RF_ari, $F_{N_{\mathrm{d}}}$ is the RF_aci due to the increase in N_d, and F_ℒ and $F_{f_{\mathrm{c}}}$ are the ℒ and cloud fraction (f_c) adjustments to the RF_ari and RF_aci. Other adjustments, e.g., due to rapid changes in land surface temperatures or atmospheric temperature and humidity profiles, have been estimated as small in a previous study (Heyn et al., 2017). Each of these terms maps fairly well onto a parameterization in the GCM: RF_ari is parameterized in the radiative transfer, $F_{N_{\mathrm{d}}}$ is parameterized in a droplet activation scheme, the ACI part of F_ℒ is parameterized in the precipitation microphysics (and, if enhanced evaporation becomes tractable in the future, that component will presumably be parameterized in the turbulence scheme; e.g., Guo et al., 2011; Neubauer et al., 2014), and $F_{f_{\mathrm{c}}}$ is parameterized in the cloud cover scheme (although in our model the response to the perturbation is indirect, via relative humidity changes subsequent to precipitation rate changes). The only component that emerges from the model dynamics rather than from an explicit parameterization is the adjustments of temperature and moisture profiles that entail further adjustments to aerosol–cloud interactions and to aerosol–radiation interactions (formerly known as “semi-direct effect”); both of these terms are small (Heyn et al., 2017; Stjern et al., 2017). A natural first step would then be to ask how large each of these terms is, and a natural second step would be to ask how much uncertainty each contributes to the total. One of the benefits of such a decomposition would be that it would provide a more solid footing for – or falsify – the notion that models agree fairly well on the supposedly simpler problem of RF_aci and not well at all on the supposedly harder problem of the adjustments. However, performing the decomposition is quite difficult in practice. Ghan (2013) addresses the issue of separating F_ari from $F_{N_{\mathrm{d}}}+F_{\mathcal{L}}+F_{f_{\mathrm{c}}}$ using the model's intrinsic knowledge of the anthropogenic perturbations. It would be desirable to separate the three ERF_aci components analogously, using the intrinsic model knowledge of the time-varying, three-dimensional aerosol perturbation and resulting perturbation of cloud properties. However, methods that attempt to do so need to contend with the problem that ERF_aci is diagnosed from two separate runs that cannot easily share fields online, so that double radiation calls as in Ghan (2013) are not feasible; furthermore, even if double radiation calls could be made, this would be of little help in estimating adjustments, which, by definition, do not respond to the aerosol perturbations instantaneously. One method to estimate ERF_aci components is the approximate partial radiative perturbation (APRP) decomposition (Taylor et al., 2007; Zelinka et al., 2014). APRP decomposes the cloud property changes into changes in area fraction, cloud albedo, and cloud absorption; the change in area fraction maps well onto the f_c adjustment in the forcing–adjustment framework, but the APRP cloud albedo change includes both the effect of the anthropogenic N_d change and the ℒ adjustment. Another method is to deactivate the parameterized precipitation suppression (and, if models include a parameterization of the enhanced evaporation, deactivate that as well); however, the model with parameterized adjustments will produce a different climate than the model without (Penner et al., 2006). Due to the complications arising in methods that directly use the model state, less direct methods have been developed that idealize the cloud field as a globally homogeneous single layer (Ghan et al., 2016) or use similar regression-based statistical techniques as satellite studies (Gryspeerdt et al., 2019b).

In this work, we apply the method of partial radiative perturbations (PRPs; Wetherald and Manabe, 1988; Colman and McAvaney, 1997; Colman, 2003; Klocke et al., 2013) to the ERF decomposition problem. PRP falls in the category of methods that directly use the intrinsic model knowledge of the time-varying, three-dimensional aerosol perturbation and resulting perturbation of cloud properties. The starting point for PRP is a perturbed and an unperturbed model run. One then introduces the fields of the perturbed run into the unperturbed run, one at a time, and reruns the radiative transfer scheme on the partially perturbed state to derive the resulting change in radiative fluxes. In our application, the two runs are simulations with an atmospheric GCM with prescribed climatological, seasonally varying sea-surface temperature (SST) and sea ice cover (SIC) distributions with present-day and preindustrial aerosol emissions, nudged to present-day large-scale upper-level winds to reduce the internal variability without overconstraining the behavior of lower-tropospheric warm cloud and to allow significant changes in cloud property to emerge after a shorter integration time than would otherwise be required (e.g., Kooperman et al., 2012; Zhang et al., 2014). The perturbed fields are N_d, ℒ, and f_c; the corresponding changes in radiative fluxes are $F_{N_{\mathrm{d}}}$, F_ℒ, and $F_{f_{\mathrm{c}}}$.

In Sect. 2, we describe the PRP method and ECHAM–HAMMOZ model in detail; in Sect. 3, we use PRP to estimate the ERF components in the ECHAM–HAMMOZ model and determine whether the adjustments are a simple proportional response to the forcing.

2 Methods

We first give a brief formal description of the PRP method; we then describe the ECHAM–HAMMOZ model configurations to which we will apply the method.

2.1 Partial radiative perturbations

We denote the shortwave TOA flux as Q and the longwave flux as R (all-sky, positive downward in both cases). For the purposes of this analysis, the radiative flux in each spectral range is considered a function of the cloud properties N_d, ℒ (or the vertically resolved analogue q_l), and f_c. The dependence of the fluxes on other climate state variables – water vapor mixing ratio q; ice-water particle size, shape, and mixing ratio q_i; aerosol properties; radiatively active gases; surface properties; and incoming solar radiation – is implicit:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }Q(\mathit{\lambda },\mathit{\varphi },t)=Q\left(N_{\mathrm{d}}\right(\mathit{\lambda },\mathit{\varphi },p,t),q_{\mathrm{l}}(\mathit{\lambda },\mathit{\varphi },p,t),f_{\mathrm{c}}(\mathit{\lambda },\mathit{\varphi },p,t\left)\right),\text{(3)}& {\displaystyle }& {\displaystyle }R(\mathit{\lambda },\mathit{\varphi },t)=R\left(N_{\mathrm{d}}\right(\mathit{\lambda },\mathit{\varphi },p,t),q_{\mathrm{l}}(\mathit{\lambda },\mathit{\varphi },p,t),f_{\mathrm{c}}(\mathit{\lambda },\mathit{\varphi },p,t\left)\right).\end{array}$$

Let ${\mathit{x}}^A=\mathit{\{}{N}_{\mathrm{d}}^A,{\mathcal{L}}^A,f_{\mathrm{c}}^A\mathit{\}}$ and ${\mathit{x}}^B=\mathit{\{}{N}_{\mathrm{d}}^B,{\mathcal{L}}^B,{f_{\mathrm{c}}}^B\mathit{\}}$ denote the cloud properties in runs A and B. We then define forward and backward PRP as inserting one cloud property at a time from one run into the cloud field of the other and recalculating the radiative fluxes:

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\delta }}_{A\to B}{Q}_{\mathit{\xi }}& {\displaystyle }=Q\left(\mathit{\right\{}{x}_{\mathit{\xi }}^A,x_{\mathit{\zeta }\ne \mathit{\xi }}^B\mathit{\left\}}\right)-Q\left({\mathit{x}}^B\right),\text{(5)}& {\displaystyle }{\mathit{\delta }}_{B\to A}{Q}_{\mathit{\xi }}& {\displaystyle }=Q\left(\mathit{\right\{}{x}_{\mathit{\xi }}^B,x_{\mathit{\zeta }\ne \mathit{\xi }}^A\mathit{\left\}}\right)-Q\left({\mathit{x}}^A\right),\end{array}$$

where δ denotes the difference in TOA flux when cloud property ξ is substituted from run A into run B or run B into run A, respectively, and Q (or R, analogously) is recalculated using the offline version of the model's radiative transfer scheme. Forward–backward PRP is simply the average of the two, taking into consideration that reversing direction reverses the sign of the radiative-flux perturbation (e.g., Klocke et al., 2013):

$$\begin{array}{}\text{(6)}& {\mathit{\delta }}_{A\leftrightarrow B}{Q}_{\mathit{\xi }}={\displaystyle \frac{{\mathit{\delta }}_{A\to B}{Q}_{\mathit{\xi }}-{\mathit{\delta }}_{B\to A}{Q}_{\mathit{\xi }}}{\mathrm{2}}}.\end{array}$$

When A denotes the preindustrial (PI)-emissions run and B denotes the present-day (PD)-emissions run, the components of ERF_aci correspond to

$$\begin{array}{}\text{(7)}& {\displaystyle }{F}_{N_{\mathrm{d}}}& {\displaystyle }=\stackrel{\mathrm{‾}}{{\mathit{\delta }}_{\text{PI}\leftrightarrow \text{PD}}{Q}_{N_{\mathrm{d}}}+{\mathit{\delta }}_{\text{PI}\leftrightarrow \text{PD}}{R}_{N_{\mathrm{d}}}},\text{(8)}& {\displaystyle }{F}_{\mathcal{L}}& {\displaystyle }=\stackrel{\mathrm{‾}}{{\mathit{\delta }}_{\text{PI}\leftrightarrow \text{PD}}{Q}_{\mathcal{L}}+{\mathit{\delta }}_{\text{PI}\leftrightarrow \text{PD}}{R}_{\mathcal{L}}},\text{(9)}& {\displaystyle }{F}_{f_{\mathrm{c}}}& {\displaystyle }=\stackrel{\mathrm{‾}}{{\mathit{\delta }}_{\text{PI}\leftrightarrow \text{PD}}{Q}_{f_{\mathrm{c}}}+{\mathit{\delta }}_{\text{PI}\leftrightarrow \text{PD}}{R}_{f_{\mathrm{c}}}}.\end{array}$$

For other meanings of A and B, as in the additional experiments performed in Sect. 3.2–3.4, the equivalent expressions to Eqs. (7)–(9) describe pseudoforcing components rather than forcing components; we denote them as ${\stackrel{\mathrm{̃}}{F}}_{N_{\mathrm{d}}}$, ${\stackrel{\mathrm{̃}}{F}}_{\mathcal{L}}$, and ${\stackrel{\mathrm{̃}}{F}}_{f_{\mathrm{c}}}$.

In Eqs. (7)–(9),

$$\begin{array}{}\text{(10)}& \stackrel{\mathrm{‾}}{\mathcal{F}}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\mathcal{F}\left(t_i\right)\end{array}$$

indicates averaging over the time dimension of a field ℱ evaluated at the N time steps $\mathit{\{}{t}_{\mathrm{1}},\mathrm{\dots },t_N\mathit{\}}$. “Evaluated” can, itself, refer to a temporal average over the interval between evaluation time steps, as in a 3-hourly or daily mean, or it can refer to the instantaneous value of the field at that time step; when the distinction matters (because Q and R are not linear functions of their input variables), we will indicate the averaging interval as ${\stackrel{\mathrm{‾}}{\mathcal{F}}}^{\left(\mathrm{\Delta }t\right)}$. Thus, ${\stackrel{\mathrm{‾}}{\mathcal{F}}}^{\left(\text{inst}\right)}$ denotes the temporal mean of instantaneous model output, while ${\stackrel{\mathrm{‾}}{\mathcal{F}}}^{\left(\mathrm{3}\text{h}\right)}$ denotes the temporal mean of 3-hourly averaged model output.

2.2 Model description

We use several model runs performed with the ECHAM–HAMMOZ model, version echam6.1–ham2.2–moz0.9 (Neubauer et al., 2014). The model is based on the ECHAM atmospheric general-circulation model (Stevens et al., 2013), the HAM interactive aerosol module (Stier et al., 2005; Zhang et al., 2012), and the trace-gas chemistry module MOZ (Kinnison et al., 2007) (the latter is disabled in these runs). Of most direct relevance to our study, the parameterized processes contributing to warm-cloud–aerosol interactions are aerosol activation into cloud droplets according to Lin and Leaitch (1997), diagnostic warm rain processes (autoconversion and accretion) according to Khairoutdinov and Kogan (2000), and aerosol scavenging according to Croft et al. (2009, 2010). The stratiform cloud scheme is that of Lohmann and Roeckner (1996), extended to double-moment microphysics by Lohmann et al. (2007) and Lohmann and Hoose (2009), with the Sundqvist et al. (1989) critical-relative-humidity cloud cover parameterization.

To reduce internal variability and achieve low statistical uncertainty on the forcing components within a reasonable integration time, we use monthly varying but yearly repeating SST and SIC from the observed climatology and nudge the large-scale wind fields to the present-day ERA-Interim reanalysis (Dee et al., 2011) wind fields of the years 2000–2010 (in some sensitivity studies, only the year 2000 is used).

Estimates of radiative forcing are computed by performing a pair of model runs with present-day SST, SIC, and wind fields, as well as aerosol (precursor) emissions estimates for either the year 2000 or the year 1850. Emissions are from the AEROCOM-II ACCMIP dataset; in particular, anthropogenic emissions follow Lamarque et al. (2010).

To perform PRP on the model output, we have updated the offline version of the RRTM-based ECHAM6 radiative transfer code (Pincus and Stevens, 2013) that was originally used in Klocke et al. (2013). We use climatological aerosol mixing ratios (Kinne et al., 2013) rather than the model's time-varying aerosol fields for aerosol–radiation interactions; this reduces the technical complexity as well as the volume of model output needed to perform the PRP calculation. To the extent that aerosol overlying cloud is a small effect, this choice mainly affects our estimate of the f_c adjustment (Ghan, 2013; Zelinka et al., 2014), which, unlike the N_d forcing and ℒ adjustment, is straightforward to compute without the PRP machinery; comparison to Gryspeerdt et al. (2019b) shows that the f_c adjustment estimate is not strongly affected by this simplification.

When clouds are absent (or the cloud fraction is very low) in one run and present in the other, perturbing cloud properties can yield unrealistically large or small q_l or N_d (and thus r_e); this decorrelation problem is well known from the application to climate feedbacks (Colman and McAvaney, 1997) in the context of the correlation between water vapor and cloudiness. We allow the radiative transfer code to resolve the conflicting cloud properties in the same way as it does when the cloud microphysics and cloud cover schemes produce conflicting cloud properties; in particular, r_e can only vary within the limits of the cloud optics lookup table used by the model ($\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}$ to $\mathrm{32}\times {\mathrm{10}}^{-\mathrm{6}}$ m). Appendices A1 and A2 describe tests we performed to verify that forward–backward PRP ERF_aci estimates are not impacted by the decorrelation problem.

3 Results

Since the components of the ERF_aci have not been diagnosed before in ECHAM–HAMMOZ by any method, we begin by presenting their global-mean values and geographic distributions in Sect. 3.1. In Sect. 3.2, we investigate whether rapid adjustments to the Twomey forcing are proportional to $F_{N_{\mathrm{d}}}$ in terms of their spatial patterns. Section 3.3 investigates the sensitivity of the PRP results to the treatment of model columns containing ice and mixed-phase clouds. In Sect. 3.4, we determine how much temporal averaging is permissible before the PRP estimate becomes inaccurate. Sections A1 and A2 discuss whether PRP diagnoses the ERF_aci components correctly in the presence of decorrelation effects (i.e., effects of introducing one cloud property from one run into an uncorrelated cloud field in another run).

3.1 What are the ERF_aci components in ECHAM–HAMMOZ?

Using the PRP decomposition, Eqs. (7)–(9), we can diagnose the contributions to the ERF_aci from PD and PI-emission fixed-SST model runs. This is shown in Fig. 1 and Table 1. The longwave effective forcing due to warm-cloud–aerosol interactions is small, as could be expected. No single forcing or adjustment dominates the shortwave ERF_aci; the global-mean forcing $F_{N_{\mathrm{d}}}$ and ℒ adjustment F_ℒ are of comparable magnitude at −0.52 and −0.53 W m⁻². The cloud-fraction adjustment $F_{f_{\mathrm{c}}}$ at −0.31 W m⁻² is the smallest of the components, consistent with CMIP5 models but not with other AeroCom models (Zelinka et al., 2014; Ghan et al., 2016; Gryspeerdt et al., 2019b). Our decomposition agrees with the results of Gryspeerdt et al. (2019b) to within 0.1 W m⁻² but disagrees with Ghan et al. (2016) (who estimate $F_{\mathcal{L}}/F_{N_{\mathrm{d}}}\approx \mathrm{5}$, in contrast to our estimate of ≈1), which both use a very similar configuration of the same version of ECHAM–HAMMOZ. Research to understand the differences between our decomposition and others is underway. All results in this section are based on 3-hourly mean output, ${\stackrel{\mathrm{‾}}{F}}^{\left(\mathrm{3}\text{h}\right)}$. This is a commonly used model output configuration, albeit at the expensive end of the spectrum from the standpoint of storage space requirements. We will justify this choice in Sect. 3.4.

Table 1 ERF_aci components in ECHAM–HAMMOZ estimated by forward–backward PRP. The total ERF also includes the ice-phase ACI effects (−0.59 W m⁻² in the SW, 0.88 W m⁻² in the LW), RF_ari (−0.17 W m⁻² in the SW), and a negligible surface-albedo contribution (−0.01 W m⁻²), estimated for a very similar model run in Gryspeerdt et al. (2019b). The sum of the components thus balances at approximately the 0.2 W m⁻² level, a relative error similar to the 0.1 W m⁻² estimated uncertainty on the ERF_aci components.

Figure 1 ERF_aci components in ECHAM–HAMMOZ estimated by forward–backward PRP.

The geographic patterns of all components exhibit similar features. In the Northern Hemisphere, fairly strong forcing prevails over both oceans and over most of the continents, with the exception of desert regions, northern Asia, and the Arctic; over the continents, a plume of high forcing components over China, extending eastwards into the Pacific Ocean, of a magnitude far greater than over Europe and North America, is especially pronounced. In the Southern Hemisphere, on the other hand, sizable forcing components are largely limited to the subtropical southern Pacific and southern Atlantic in the vicinity of the persistent stratocumulus decks; smaller local maxima in the forcing components also exist in the outflow regions of the midlatitude westerlies downwind of South America, Africa, and Australia.

These geographic patterns result from a convolution of the distributions of susceptible clouds and N_d perturbations. Figure 2 shows the ERF_aci sensitivity, defined as the forcing or adjustment strength per e-folding of the N_d burden (with the bar denoting temporal averaging over the length of the run),

$$\begin{array}{}\text{(11)}& S_{\mathit{\xi }}={\displaystyle \frac{F_{\mathit{\xi }}}{\mathrm{\Delta }\ln {\stackrel{\mathrm{‾}}{\mathcal{N}}}_{\mathrm{d}}}},\end{array}$$

where the N_d burden is defined as

$$\begin{array}{}\text{(12)}& {\mathcal{N}}_{\mathrm{d}}=\int N_{\mathrm{d}}\mathrm{d}z;\end{array}$$

Equation (11) is similar to the factorization of Bellouin et al. (2019). The sensitivity for all ERF_aci components is greater over ocean than over land. Over land, the most susceptible clouds can be found over tropical Africa, southeastern Asia, and Central America; the anthropogenic N_d perturbation over South America is too small to determine reliable sensitivities. Over ocean, the largest sensitivities occur in the midlatitudes, near the Equator, and over the eastern ocean boundaries; again, N_d perturbations in the Southern Hemisphere are too small to estimate sensitivities reliably. The regions of maximum sensitivity are the regions where the base-state ℒ is high (not shown), which is intuitive in light of the strong influence of ℒ on cloud albedo.

Figure 2 Sensitivity of ERF_aci components to anthropogenic 𝒩_d change (shown only where $\mathrm{\Delta }\ln {\stackrel{\mathrm{‾}}{\mathcal{N}}}_{\mathrm{d}}>\mathrm{0.05}$).

Figure 3 Anthropogenic ln 𝒩_d change.

Figure 3 shows the anthropogenic perturbation in 𝒩_d. The perturbation is generally stronger over the continents than over ocean; over ocean, it is generally strongest in the subtropics. Thus, the highest perturbations tend to coincide with the least susceptible clouds, which explains why the ERF_aci components exhibit far less variability over a wide latitudinal range over the Northern Hemisphere oceans than either the sensitivity or the N_d perturbation, and why the remote Northern Hemisphere oceans are comparable to the polluted continental regions in ERF_aci strength. The exceptions to this pattern are the near-shore eastern ocean boundaries, where a high sensitivity and reasonably strong N_d perturbation coincide, and eastern China, where a reasonably high sensitivity and very strong N_d perturbation coincide. The location where ECHAM–HAMMOZ simulates both the strongest N_d forcing and the strongest ℒ adjustment is over land in China, extending downwind into the northwestern Pacific Ocean.

In observational studies or observationally constrained modeling studies, it is common to define susceptibilities analogously to Eq. (11) based on PD variability in cloud and aerosol variables and then multiply those susceptibilities by wholly or partially (Bellouin et al., 2013; Kinne, 2019) model-derived estimates of anthropogenic aerosol perturbations. (In the terminology we adopt here, “sensitivity” is a change in cloud property or cloud radiative effect in response to a climatological change in an aerosol variable, whereas “susceptibility” is a change in response to an instantaneous change in an aerosol variable.) Our results show that the most susceptible oceanic clouds in this version of ECHAM–HAMMOZ occur where continental pollution intrudes on relatively clean conditions over the eastern ocean boundaries. Some observational (or observationally constrained) studies agree with this result (Quaas et al., 2008; Alterskjær et al., 2012; Engström et al., 2015; Gryspeerdt et al., 2016; Andersen et al., 2017; Christensen et al., 2017; McCoy et al., 2017), while others disagree (Lebsock et al., 2008; Chen et al., 2014; Gryspeerdt et al., 2017). Of those studies that are not restricted to oceanic clouds, some agree with our finding of strong forcing due to relatively susceptible clouds over China (Gryspeerdt et al., 2017; McCoy et al., 2017) and some do not (Quaas et al., 2008; Alterskjær et al., 2012; Gryspeerdt et al., 2016).

3.2 Are the adjustments proportional to the forcing?

An intriguing aspect of the ACI problem is whether the adjustments may be described approximately as a proportional response to the forcing (Gryspeerdt et al., 2019a). On the one hand, we do not necessarily expect proportionality in the physical atmosphere, since the processes responsible for the adjustments carry memory of the cloud evolution over various timescales; the parameterized cloud processes in GCMs share this feature, at least in principle, since the anthropogenic N_d perturbation seen by the precipitation parameterization at one time step could be the result of a CCN perturbation at some point in the past, carried to another point in space by advection, and influenced by any of the other parameterized cloud processes. On the other hand, complex systems often exhibit simple emergent behaviors (e.g., Mülmenstädt and Feingold, 2018, and references therein). If the adjustments were to follow proportionally from the forcing, one consequence for the ACI problem would be that the total ERF_aci uncertainty should not be estimated as the uncertainty on the sum of uncorrelated RF_aci and adjustments but rather take the correlation between the forcing and adjustments into account, which would result in a smaller ERF_aci uncertainty estimate.

Figure 4 ERF_aci adjustments relative to the Twomey forcing.

In this study, we can test for proportionality in terms of the geographic distribution using the spatial variability in the temporal-mean ERF_aci components. Figure 4 shows that the zonal mean of the ratio between ℒ adjustment and N_d forcing is relatively stable around unity between the southern and northern midlatitudes, with fairly small interhemispheric differences except in the Southern Ocean. The picture is somewhat different for the ratio between f_c adjustment and N_d forcing, which is more latitudinally variable and more different between the Northern and Southern Hemisphere. Figure 5 reinforces these conclusions, showing that F_ℒ and $F_{N_{\mathrm{d}}}$ are fairly tightly correlated with a regression relationship remarkably close to one to one, while the relationship between $F_{f_{\mathrm{c}}}$ and $F_{N_{\mathrm{d}}}$ is much looser.

Figure 5 Correlation plots between the temporal-mean Twomey forcing and the adjustments; color indicates the number of grid boxes within each $\mathrm{0.05}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}\times \mathrm{0.05}\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ bin; the red line is a linear least-squares regression; the blue line is a generalized additive model regression (Wood, 2011), with 95 % confidence interval shaded in light blue; and the dashed gray line is the one-to-one line.

One interpretation of these results is that the ERF_aci components share a geographic pattern due to the fact that large effects result from the coincidence of large anthropogenic aerosol sources and susceptible clouds; the shared geographic pattern then leads to an approximately proportional relationship that breaks down farther from the source regions or where a different mixture of cloud processes dominates the cloud response (e.g., the Southern Ocean). The cloud cover scheme, which diagnoses f_c from the grid mean relative humidity, to some extent decouples f_c from the other cloud properties, which attenuates the influence of N_d on f_c. Nevertheless, the vagueness of this argument is unsatisfactorily mismatched against the precision of the F_ℒ–$F_{N_{\mathrm{d}}}$ relationship, which suggests a deeper mechanism at play, e.g., that precipitation acts as a common sink process for both N_d and ℒ.

Further evidence for proportionality comes from Gryspeerdt et al. (2019b), who find an intermodel proportional relationship between global-mean forcing and rapid adjustments.

3.3 How should we treat columns containing ice?

In attempting to diagnose warm-cloud ACI forcing components, the question of how ice-containing clouds should be handled arises. We can conduct the following set of experiments to determine the range of forcing strengths associated with different thermodynamic-phase treatments.

1. Perturb cloud properties in all cloudy model levels.

2. Perturb cloud properties in any liquid-containing cloudy model levels.

3. Perturb cloud properties in liquid-only cloudy model levels (default).

4. Perturb cloud properties in liquid-only cloudy columns.

5. Perturb cloud properties in liquid-only cloudy columns, correcting the result by their temporal occurrence fraction.

Table 2 Dependence of diagnosed ERF_aci components on treatment of thermodynamic phase.

Table 2 summarizes the results. (For reasons of efficiency, we performed these sensitivity experiments on daily-mean output, shortwave flux only. We did not perform experiment 2 because we expect the result to lie between experiments 1 and 3, whose separation is already in the noise.) We conclude that how we choose to treat mixed-phase and ice clouds makes little difference in ECHAM–HAMMOZ, so long as we do not restrict ourselves to columns containing only warm clouds. In the latter case, correcting the forcing by the temporal occurrence fraction of liquid-only columns in each model latitude–longitude box approximately recovers the results when ice-containing columns are retained; however, there is some indication of diverging trends in $F_{N_{\mathrm{d}}}$ (which decreases in magnitude as the ice filtering becomes more restrictive) and F_ℒ (which increases in magnitude as the ice filtering becomes more restrictive). The ice-free column requirement is often made in passive remote sensing studies to prevent contamination from ice clouds overlying warm clouds and uncertainties in multilayer cloud retrievals.

3.4 Does temporal averaging bias the results?

As Table 3 shows, longer averaging periods underestimate the forcing, but the differences between instantaneous output (the model time step is 7.5 min, but we sample every 3 h to reduce the data volume) and 3 h averages are minimal. Multimodel ensembles which archive 3 h average output or 3 h subsampled instantaneous output of column cloud properties, e.g., AeroCom and CFMIP2, are therefore amenable to treatment by the PRP method.

Table 3 Dependence of diagnosed ERF_aci components on temporal averaging.

4 Conclusions

We have presented the first decomposition of the ACI effective forcing in ECHAM–HAMMOZ into a Twomey forcing component and rapid adjustments of ℒ and f_c. In ECHAM–HAMMOZ, no single component dominates: $F_{N_{\mathrm{d}}}=-\mathrm{0.52}$ W m⁻², $F_{\mathcal{L}}=-\mathrm{0.53}$ W m⁻², and $F_{f_{\mathrm{c}}}=-\mathrm{0.31}$ W m⁻²; the Twomey forcing and ℒ adjustment are approximately equally strong, and the f_c adjustment is somewhat weaker, as in many other models. The global ERF is dominated by the Northern Hemisphere forcing. Within the Northern Hemisphere, the strongest forcing components occur over land in China in $F_{N_{\mathrm{d}}}$ and F_ℒ. As expected, the stratocumulus sheets over the eastern ocean basins also show strong responses in both hemispheres, as do the midlatitude North Atlantic and North Pacific.

The temporal-mean spatial patterns of $F_{N_{\mathrm{d}}}$ and F_ℒ are highly correlated, suggesting an effective proportionality in the ℒ adjustment to the Twomey forcing even though the precipitation-suppression mechanism by which the ℒ adjustment is parameterized in the model has inherent memory that could decouple it from the Twomey effect. The spatial patterns of the temporal-mean $F_{N_{\mathrm{d}}}$ and F_ℒ, while sharing some of the same gross features, have a much less tight relationship than $F_{N_{\mathrm{d}}}$ and F_ℒ.

In our study of ECHAM–HAMMOZ, the forcing components are fairly insensitive to how we treat columns containing both ice and liquid cloud condensate. Requiring that columns be free of ice and then correcting for the temporal fractional occurrence of ice cloud, a technique that is often necessary in observational studies, largely reproduces the results we obtain when we do not filter out such columns, albeit possibly causing an overestimate of the ℒ adjustment and an underestimate of the N_d forcing. (In interpreting the bearing of these results on analyses of satellite cloud retrievals, note that these studies do not necessarily apply the ice-free requirement at the coarse GCM scales of the present work, depending on whether they use gridded level 3 data or the level 2 native resolution of the retrieval algorithms.)

Through idealized sensitivity studies presented in the Appendix, we have showed that PRP is a viable method for accurately decomposing ERF_aci into a N_d forcing and ℒ and f_c adjustments. This is the case despite large artifacts that occur due to the decorrelated cloud property fields; the forward–backward technique advocated by Colman and McAvaney (1997) is capable of removing these artifacts.

PRP directly uses the intrinsic model knowledge of the time-varying, three-dimensional aerosol perturbation and resulting perturbation of cloud properties to diagnose the ERF_aci components and their spatial patterns. This makes it a useful tool for providing context to other less-resource-intensive decomposition methods (e.g., Ghan et al., 2016; Gryspeerdt et al., 2019b) or to intercomparison studies (e.g., Pincus et al., 2016; Smith et al., 2018) despite its demand for high-frequency vertically resolved model output.

Appendix A: Validation of the PRP method for ACI decomposition

A1 What is the effect of decorrelating the cloud properties?

Consider the results of forward and backward PRP plotted separately for the PD–PI experiment in Fig. A1. Not only are the magnitudes grotesque, but taken at face value, they would imply a positive forcing in one direction and a negative forcing in the other. Furthermore, the spatial patterns bear no resemblance to that expected for ERF_aci. In this section, we investigate the consequences of these features for the ERF_aci decomposition.

Figure A1 Forward (a) and backward (b) PRP estimates of the ERF_aci components. Note the significantly wider color scale than in Fig. 1.

Figure A2 Forward (a), backward (b), and forward–backward (c) PRP performed on the same-climate–different-weather case. Note the significantly wider color scale than in Fig. 1.

Any given atmospheric property is often correlated with many others. Substituting cloud properties one at a time breaks these correlations. For example, since ECHAM–HAMMOZ parameterizes precipitation suppression by aerosol, we expect a positive correlation between N_d and ℒ within a model run. If we substitute ℒ from another run, the mechanistic link between N_d and ℒ through precipitation suppression, by which higher N_d at a given point in time leads to higher ℒ at later times, is broken, and, therefore, the correlation between N_d and ℒ is altered.

We estimate the strength of such decorrelation effects by performing two model runs with (almost exactly) the same model physics, both nudged to the same large-scale dynamics and with the same fixed SST; the only difference between the runs is that a parameter in the Khairoutdinov and Kogan (2000) formulation for the autoconversion rate, tuned for ECHAM–HAMMOZ,

$$\begin{array}{}\text{(A1)}& Q_{\text{aut}}=\mathit{\gamma }{q}_{\mathrm{l}}^{\mathit{\alpha }}{\left({\displaystyle \frac{N_{\mathrm{d}}}{\mathrm{1}{\mathrm{cm}}^{-\mathrm{3}}}}\right)}^{-\mathit{\beta }},\end{array}$$

has been slightly perturbed from β=1.79 to ${\mathit{\beta }}^{\prime }=\mathrm{1.79}+{\mathrm{10}}^{-\mathrm{5}}$ (α=2.47 and $\mathit{\gamma }=\mathrm{4}\times \mathrm{1350}{\mathrm{s}}^{-\mathrm{1}}$ are unchanged). Even over short integration times (a year), these model runs will converge on the same climate, with nearly identical forcing components. (The small perturbation in β does not result in a significant change in model sensitivity.) However, at any given elapsed integration time and geographic location, the cloud properties in the two runs are essentially uncorrelated. We refer to this pair of runs as the same-climate–different-weather experiment. Knowing that the true climatological TOA flux difference between this pair of runs is zero, we can use these runs to estimate decorrelation effects between any other decorrelated pair of runs, including the PD and PI emissions runs.

We find that decorrelation effects cause the PRP method to estimate enormous TOA flux perturbations when we perform forward or backward substitution of any single cloud property; this is shown in panels (a) and (b) of Fig. A2.

Unlike forward PRP or backward PRP individually, forward–backward PRP is unaffected by decorrelation, both in the global mean and locally in the temporal mean: panel (c) of Fig. A2 shows that the fluctuations in ΔQ rapidly (i.e., within a year) average to zero. This confirms that the Colman and McAvaney (1997) prescription is successful at minimizing the spurious effects of decorrelation.

A2 Does PRP give the right answer?

The preceding section provides evidence that strong decorrelation effects do not lead to a spurious offset in forward–backward PRP results. Next, we show that decorrelation effects also do not lead to spurious scale factors. To do so, we scale N_d and ℒ by a globally constant factor of 1.1 at all time steps and scale f_c by 0.99. We use PRP to diagnose the forcing associated with each of these perturbations; the results are shown in the first three rows of Table A1. We can then estimate the strength of the decorrelation effects by performing PRP on the ${\mathit{\beta }}^{\prime }=\mathrm{1.79}+{\mathrm{10}}^{-\mathrm{5}}$ run and the scaled-$\mathit{\{}{N}_{\mathrm{d}},\mathcal{L},f_{\mathrm{c}}\mathit{\}}$ β=1.79 run. This is shown in the middle three rows of Table A1; the correct results are recovered to good approximation, with generally small attribution to incorrect ERF_aci components (the largest is −0.05 W m⁻² incorrectly diagnosed as f_c adjustment in the ℒ×1.1 experiment) and generally small differences between the actual and diagnosed ERF_aci components (the largest difference is a diagnosed $F_{\mathcal{L}}=-\mathrm{0.48}$ W m⁻² when the correct value is −0.53 W m⁻²). The final test is an experiment in which all cloud properties are perturbed simultaneously and the clouds are decorrelated by using a ${\mathit{\beta }}^{\prime }=\mathrm{1.79}+{\mathrm{10}}^{-\mathrm{5}}$ baseline run. The results are shown on the last line of Table A1; the correct ERF_aci components are recovered in the presence of the confounding effects of decorrelation and of perturbing multiple cloud properties simultaneously with 0.1 W m⁻² or better accuracy, with the largest discrepancy being the diagnosed $F_{f_{\mathrm{c}}}=\mathrm{0.14}$ W m⁻² when the correct value is 0.24 W m⁻²).

Thus, we find that forward–backward PRP can diagnose the forcing components correctly in the presence of decorrelations, in addition to diagnosing the absence of forcing correctly in the same-climate–different-weather case.

A3 Are the results sensitive to choosing grid-mean or in-cloud perturbations?

Perturbing in-cloud or grid-mean N_d and ℒ would be equivalent in the limit in which TOA flux perturbations are linear in the cloud properties. While individual model columns do not satisfy this linearity requirement, the temporal mean apparently exhibits sufficient effective linearity that the choice of in-cloud or grid-mean perturbations has little effect on the ERF component estimate; compare Tables 1 and A2.

Table A1 ERF_aci components resulting from idealized perturbations to N_d, ℒ, and f_c; estimate of the same ERF_aci components by forward–backward PRP in the presence of decorrelation effects.

Table A2 ERF_aci components calculated by PRP on f_c and in-cloud N_d and q_l.