General circulation models simulate negative liquid water path–droplet number correlations, but anthropogenic aerosols still increase simulated liquid water path

Abstract

General circulation models' (GCMs) estimates of the liquid water path adjustment to anthropogenic aerosol emissions differ in sign from other lines of evidence. This reduces confidence in estimates of the effective radiative forcing of the climate by aerosol–cloud interactions (ERFaci). The discrepancy is thought to stem in part from GCMs' inability to represent the turbulence–microphysics interactions in cloud-top entrainment, a mechanism that leads to a reduction in liquid water in response to an anthropogenic increase in aerosols. In the real atmosphere, enhanced cloud-top entrainment is thought to be the dominant adjustment mechanism for liquid water path, weakening the overall ERFaci. We show that the latest generation of GCMs includes models that produce a negative correlation between the present-day cloud droplet number and liquid water path, a key piece of observational evidence supporting liquid water path reduction by anthropogenic aerosols and one that earlier-generation GCMs could not reproduce. However, even in GCMs with this negative correlation, the increase in anthropogenic aerosols from preindustrial to present-day values still leads to an increase in the simulated liquid water path due to the parameterized precipitation suppression mechanism. This adds to the evidence that correlations in the present-day climate are not necessarily causal. We investigate sources of confounding to explain the noncausal correlation between liquid water path and droplet number. These results are a reminder that assessments of climate parameters based on multiple lines of evidence must carefully consider the complementary strengths of different lines when the lines disagree.

1 Introduction

Aerosol–cloud interactions (ACIs) remain the greatest source of uncertainty in our estimates of anthropogenic perturbations in the Earth's energy budget . In liquid clouds, an anthropogenic aerosol perturbation essentially instantaneously alters the number of cloud droplets ($N_{\mathrm{d}}$), changing cloud reflectance and thus the shortwave radiation absorbed by the climate system, which exerts a radiative forcing on climate (radiative forcing by aerosol–cloud interactions or RFaci; ). While our knowledge of RFaci is uncertain , an even thornier issue is cloud adjustments to the $N_{\mathrm{d}}$ perturbation, where multiple processes acting at different scales from cloud droplets to planetary circulation result in a multiscale dynamics prediction problem that is impervious to any one silver bullet solution . Estimates of ACI adjustments are, therefore, based on multiple, and often conflicting, lines of evidence . Those lines of evidence are, broadly, modeling at the cloud process scale (large eddy simulation or LES), global modeling, and observations at different scales.

In the following, we focus on stratocumulus (Sc) clouds, which play a large role in the energy budget due to their high albedo and frequent occurrence. Our understanding of adjustments in Sc is that two effects compete: an anthropogenic increase in $N_{\mathrm{d}}$ suppresses precipitation , increasing cloud liquid water path ($\mathcal{L}$); but the $N_{\mathrm{d}}$ increase also promotes increasing turbulent entrainment of subsaturated air at cloud top , decreasing $\mathcal{L}$. These mechanisms are regime-dependent; precipitation suppression only plays a role in clouds that would have precipitated in the absence of the aerosol perturbation, and the entrainment mechanisms depend strongly on the turbulence generation mechanisms (for example, cloud-top radiative cooling). The regime dependence of the underlying processes leads to “process fingerprints” in $N_{\mathrm{d}}$–$\mathcal{L}$ space in LES for the very limited set of boundary conditions where LES is available. Similar bifurcation behavior appears in satellite observations, where mean $\mathcal{L}$ as a function of $N_{\mathrm{d}}$ first increases in precipitating clouds, then reaches a peak that roughly coincides with the transition to non-precipitating clouds, and then decreases again . There is evidence that this inverted-V relationship between $\mathcal{L}$ and $N_{\mathrm{d}}$ overestimates the strength of the causal effect of $N_{\mathrm{d}}$ on $\mathcal{L}$ , but qualitatively it is consistent with process understanding from LES. Integrated over all meteorological boundary conditions, the overall satellite correlation between $\mathcal{L}$ and $N_{\mathrm{d}}$ is negative. The satellite inverted V; satellite observations of natural laboratories , where the origin of the perturbation is evident ; and process-modeling lines of evidence lead to the assessment that the adjustment of $\mathcal{L}$ to anthropogenic aerosol is a reduction in $\mathcal{L}$ – that is, a positive contribution to the effective radiative forcing by ACI (ERFaci; ).

Global climate models – which, currently, means general circulation models (GCMs) run at a roughly 1° latitude–longitude spatial resolution – tell a different story. They would project an increase rather than a decrease in $\mathcal{L}$ when aerosols are increased from preindustrial (PI) to present-day (PD) concentrations . The GCM line of evidence is discounted in multiline assessments because it conflicts with the other lines and because those lines are assumed to provide more reliable information. This assumption rests on the representation of the relevant processes in GCMs. In these models, precipitation is initiated by a microphysical parameterization with an explicit dependence on $N_{\mathrm{d}}$ (or, largely equivalent, droplet size) so that the $\mathcal{L}$ increase by precipitation suppression is explicitly parameterized. Reduced $\mathcal{L}$ by enhanced evaporation, on the other hand, depends critically on meter-scale or smaller interactions between turbulence, radiation, and microphysics at the cloud edge. These interactions fall between several parameterizations and are therefore tricky to formulate in GCMs. As a perverse consequence, this causes us to fret that GCMs may be structurally incapable of representing turbulent entrainment scales, while we often mistakenly consider the many-orders-of-magnitude-smaller-scale precipitation processes a parametric problem (e.g., ).

In this work, we show that some GCMs from the Coupled Model Intercomparison Project 6 (CMIP6) era, unlike earlier model generations, are capable of producing inverted-V $N_{\mathrm{d}}$–$\mathcal{L}$ relationships in agreement with global observations and LES. Based on these PD correlations and on the $N_{\mathrm{d}}$ change between PI and PD (i.e., mimicking the information available to observation-based ERFaci estimates), these models predict a reduction in $\mathcal{L}$, which is consistent with assessments that use multiple lines of evidence. However, the causal effect of anthropogenic $N_{\mathrm{d}}$ changes on $\mathcal{L}$, as diagnosed by model experiments where all climatic boundary conditions apart from aerosols are held fixed, remains as in previous GCM generations: an anthropogenic $N_{\mathrm{d}}$ increase leads to an increase in average $\mathcal{L}$, consistent with a dominant role of the precipitation suppression mechanism parameterized in the model microphysics.

2 Data and methods

We use an ensemble of GCMs to perform fixed-sea-surface-temperature model experiments with PD and PI emissions and archive instantaneous aerosol and cloud information, with a sufficient frequency (3 h) to resolve the diurnal cycle and with a sufficient length (1–5 years with the large-scale winds nudged to PD meteorology) to draw statistically robust conclusions. The model ensembles used are the CMIP5-era AeroCom indirect effect experiment (AeroCom IND3) simulations on the one hand and four newer-version models prepared for CMIP6 on the other. The AeroCom models are described in and . The CMIP6-era models are the US Department of Energy's Energy Exascale Earth System Model (E3SM) Atmosphere Model version 1 (EAMv1; ), the NASA Goddard Institute for Space Studies (GISS) ModelE3 configuration Tun2, the Geophysical Fluid Dynamics Laboratory (GFDL) atmospheric model AM4.0 , and the Community Earth System Model version 2 Community Atmosphere Model version 6 (CESM2-CAM6; ). The CMIP6-era models were run for 1 year for the baseline experiment. E3SM was further run for 5 years for additional experiments that needed more data to perform stratification by confounding variables (see Sect. ). For E3SM, the Cloud Feedback Model Intercomparison Project (CFMIP) Observation Simulator Package (COSP) satellite simulator , mimicking the MODerate Resolution Imaging Spectroradiometer (MODIS) cloud retrievals , and a number of vertically resolved fields were archived over a limited area over the northeast Pacific (NEP) Sc region for further analysis of confounders; only this experiment uses COSP-simulated $N_{\mathrm{d}}$ and $\mathcal{L}$.

2.1 Cloud selection

From the model output, we select liquid clouds, defined by the absence of ice (ice water path $<{\mathrm{10}}^{-\mathrm{3}}$ kg m${}^{-\mathrm{2}}$ and ice cloud cover $=\mathrm{0}$) in the column. To mimic passive satellite analyses and to simplify the application of entrainment diagnostics in Part 2 of the series , we require near-overcast (liquid cloud cover $>\mathrm{0.9}$) conditions. For these liquid clouds, we calculate “in-cloud” cloud-top $N_{\mathrm{d}}$ and $\mathcal{L}$ by dividing the grid-mean $N_{\mathrm{d}}$ and $\mathcal{L}$ by the projected cloud cover. Only clouds over ocean are considered in this analysis. We refer to these clouds as “overcast clouds”.

In addition to these globally occurring overcast clouds, we also study smaller cloud subsets defined by the dynamical regime following . In this classification, the stratocumulus regime is based on vertical velocity ($\mathit{\omega }$) at 700 and 500 hPa (${\mathit{\omega }}_{\mathrm{700}}>\mathrm{10}$ hPa d${}^{-\mathrm{1}}$ and ${\mathit{\omega }}_{\mathrm{500}}>\mathrm{10}$ hPa d${}^{-\mathrm{1}}$) and lower tropospheric stability (LTS), which we define here as the difference in potential temperature ($\mathit{\theta }$) between 1000 and 700 hPa (${\mathit{\theta }}_{\mathrm{700}}-{\mathit{\theta }}_{\mathrm{1000}}>\mathrm{18.55}$ K). We further restrict the clouds to occur in grid boxes where these conditions are met at least 30 % of the time, a subjective choice that selects the subtropical Sc regions and rejects midlatitude Sc. The occurrence fraction ($f_{\text{Sc}}$) of these conditions is shown in Fig. . In addition to the requirements, all of the abovementioned warm-cloud criteria are applied. We refer to these clouds as Sc-regime clouds.

Figure 1

Occurrence fraction, $f_{\text{Sc}}$, of Sc conditions by the criteria in EAMv1, shown where $f_{\text{Sc}}>\mathrm{0.1}$. The $f_{\text{Sc}}>\mathrm{0.3}$ threshold used in the analysis selects the NEP, southeast Pacific, and southeast Atlantic Sc regions, while limiting the northward extent of the NEP region beyond the subtropics and rejecting all midlatitude Sc.

[Figure omitted. See PDF]

2.2 Analysis methods

From $\mathcal{L}$ and $N_{\mathrm{d}}$, we construct the conditional probability $P\left(\mathcal{L}\right|N_{\mathrm{d}})$, following . For ease of comparison among models and configurations, we collapse the two-dimensional $P\left(\mathcal{L}\right|N_{\mathrm{d}})$ into one dimension by calculating the geometric mean $\mathcal{L}$ in each $N_{\mathrm{d}}$ bin, also following .

For the MODIS simulator analysis in Sect. , we transform the simulated $\mathit{\tau }$ and droplet effective radius ($r_{\mathrm{e}}$) into $N_{\mathrm{d}}$ and $\mathcal{L}$ using a power-law relationship for adiabatic updrafts with constant $N_{\mathrm{d}}$ : $$\begin{array}{}\text{1}& {\displaystyle }{N}_{\mathrm{d}}={\displaystyle \frac{\sqrt{\mathrm{5}}}{\mathrm{2}\mathit{\pi }k\sqrt{{\mathit{\rho }}_{\mathrm{w}}Q}}}\sqrt{f_{\text{ad}}\mathrm{\Gamma }}{\mathit{\tau }}^{\mathrm{1}/\mathrm{2}}{r}_{\mathrm{e}}^{-\mathrm{5}/\mathrm{2}},\text{2}& {\displaystyle }\mathcal{L}={\displaystyle \frac{\mathrm{5}}{\mathrm{9}}}{\mathit{\rho }}_{\mathrm{w}}\mathit{\tau }{r}_{\mathrm{e}},\end{array}$$ where we take the ratio $k=(r_{\mathrm{v}}/r_{\mathrm{e}}{)}^{\mathrm{3}}$ between the volumetric mean radius $r_{\mathrm{v}}$ cubed and effective radius cubed to be 1; subadiabatic factor, $f_{\text{ad}}$, to be 1; scattering efficiency, $Q$ to be 2; and adiabatic condensation rate, $\mathrm{\Gamma }$, to be $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}$ kg m${}^{-\mathrm{4}}$. These assumptions minimize the complications involved in showing results that are mostly power-law behavior independent of these constant factors. (This does neglect important modifications that can arise if these factors are not, in fact, constant; .)

To analyze confounding by planetary boundary-layer (PBL) depth, $h$ (Sect. ), we identify the top of the Sc-like boundary layer by the first model level where temperature increases with height in Sc-regime overcast columns. This produces well-mixed profiles of the liquid water potential temperature, ${\mathit{\theta }}_{\mathrm{l}}$, and total water mixing ratio, $q_{\mathrm{w}}$. (Other definitions of PBL top – i.e., the model level of the greatest gradient in ${\mathit{\theta }}_l$ or $q_{\mathrm{w}}$ – yield very similar results.) As we see in Sect. , cloud and aerosol properties are remarkably stratified by PBL depth in E3SM; to keep the properties as distinct as possible as a function of PBL depth, we retain the native model vertical discretization instead of converting the hybrid pressure levels to pressure or geometric height.

Table  summarizes the emissions, cloud selection, and model run duration for each experiment.

Table 1

GCM experiments used in this analysis.

Experiment	Emissions	Cloud selection	Duration
Multimodel PD	Present-day	Global overcast	2010, nudged
Multimodel PI	Preindustrial	Global overcast	2011, nudged
E3SM $+$ COSP simulator	Present-day	NEP Sc regime	2010, nudged
E3SM precipitation	Present-day	Sc regime	2010–2014, nudged

3 Results and discussion

In Fig. , we show the behavior of the AeroCom IND3 (CMIP5-era) GCMs in $N_{\mathrm{d}}$–$\mathcal{L}$ space: with the exception of one model, $\mathcal{L}$ increases monotonically as a function of $N_{\mathrm{d}}$. In some models, the slope decreases at high $N_{\mathrm{d}}$, but only one model, the Hadley Center Global Environment Model version 3 (HadGEM), has quantitatively similar behavior to the inverted-V satellite $N_{\mathrm{d}}$–$\mathcal{L}$ plot. The behavior of these models (with the exception of HadGEM) is consistent with the interpretation that the predominant mechanism linking $\mathcal{L}$ and $N_{\mathrm{d}}$ is precipitation suppression.

Figure 2

AeroCom IND3 state-of-the-art models' $N_{\mathrm{d}}$–$\mathcal{L}$ relationship. The satellite inverted-V relationship is indicated by the solid line.

[Figure omitted. See PDF]

3.1

CMIP6-era models produce inverted-V $N_{\mathrm{d}}$–$\mathcal{L}$ relationships

A funny thing happened on the way to CMIP6: three of the four US CMIP6-era GCMs had an inverted V with a pronounced negative slope. The behavior of these models is contrasted with the AeroCom models' behavior in Fig. . The geographic distribution of the regression slope between $\log \mathcal{L}$ and $\log N_{\mathrm{d}}$ is predominantly negative in the models with an inverted V (Fig. ), as found in satellite retrievals.

Figure 3

(a) AeroCom IND3 state-of-the-art models' marginal probability distributions and $N_{\mathrm{d}}$–$\mathcal{L}$ conditional probability distributions compared to the (b) CMIP6-era state-of-the-art models' $N_{\mathrm{d}}$–$\mathcal{L}$ relationship. The satellite inverted-V relationship is indicated by the dashed line. Three of the four CMIP6 models examined are qualitatively similar to the satellite result in the sense that the $N_{\mathrm{d}}$–$\mathcal{L}$ correlation turns negative at moderate $N_{\mathrm{d}}$.

[Figure omitted. See PDF]

Figure 4

Geographic distribution of $\mathrm{d}\log \mathcal{L}/\mathrm{d}\log N_{\mathrm{d}}$ in the CMIP6-era models. Model output is aggregated to 5° $\times$ 5° latitude–longitude boxes before calculating linear regression slopes of $\log \mathcal{L}$ against $\log N_{\mathrm{d}}$.

[Figure omitted. See PDF]

One of these models (ModelE) was designed to better represent the entrainment behavior to which the negative slope is attributed in process-scale modeling. The other two (CAM6 and EAMv1), however, were not; if the negative slope is due to an entrainment ACI mechanism, it is an emergent behavior not explicitly parameterized into the turbulence scheme. It is doubly surprising that these models produce a negative slope considering that their closely related predecessor, CAM5.3-CLUBB-MG2, was part of the AeroCom ensemble and showed, at best, a slightly negative relationship between $N_{\mathrm{d}}$ and $\mathcal{L}$.

Unfortunately, the intersection between CMIP5 AeroCom models and CMIP6 models in this study is small (only CAM5 and CAM6). We hope that an updated AeroCom experiment will provide comparisons between the CMIP5 and CMIP6 versions of more models. We also hope perturbed physics ensembles of each of the inverted-V models will explore the effect of physics choices on the $N_{\mathrm{d}}$–$\mathcal{L}$ relationship.

3.2

The negative correlation between $N_{\mathrm{d}}$ and $\mathcal{L}$ does not predict the sign of PI to PD change in $\mathcal{L}$

The bulk of the $N_{\mathrm{d}}$ population lies in the part of the inverted V with a negative $N_{\mathrm{d}}$–$\mathcal{L}$ correlation. If we regarded this relationship as indicative of a causal influence of $N_{\mathrm{d}}$ on $\mathcal{L}$ – that is, that an increase in $N_{\mathrm{d}}$ causes $\mathcal{L}$ to decrease – then we would predict a decrease in $\mathcal{L}$ as $N_{\mathrm{d}}$ increases from its PI value to its PD value due to anthropogenic emissions.

We can compare the change in $\mathcal{L}$ predicted by the $N_{\mathrm{d}}$–$\mathcal{L}$ correlation in PD internal variability to the outcome of a model experiment designed to measure the causal effect of $N_{\mathrm{d}}$ on $\mathcal{L}$. This experiment fixes all climatic boundary conditions affecting cloud state (i.e., solar constant, greenhouse gases, and sea-surface temperature) with the exception of anthropogenic aerosols. The change in $\mathcal{L}$ in this experiment can therefore only be due to the anthropogenic aerosol emissions change. This model experiment shows that the causal effect of the $N_{\mathrm{d}}$ increase is to increase $\mathcal{L}$ on average, contradicting the prediction of a decrease in $\mathcal{L}$ based on PD internal variability (Fig. ). The correlation seen in PD internal variability in these models therefore cannot be causal. Plotting the correlations within PD and PI, as shown in Fig. , provides a glimpse at what is happening instead: a secular increase in $N_{\mathrm{d}}$ does not lead to a secular reduction in $\mathcal{L}$ by shifting the $\mathcal{L}$ population along the correlation line, as would be expected for a causal relationship. Instead, the correlation line shifts along with the secular shifts in $N_{\mathrm{d}}$ and $\mathcal{L}$ (mostly to the right given that the change in $N_{\mathrm{d}}$ is far greater than the change in $\mathcal{L}$) in a way that is not predicted by the correlation line itself.

Figure 5

PI–PD $\mathcal{L}$ change from the causal experiment (solid arrow) contrasted with the change predicted from the PD internal variability (dashed arrow). The mean $\log \mathcal{L}$ as a function of $\log N_{\mathrm{d}}$ from the PD model run (Fig. ) is used to predict PI mean $\log \mathcal{L}$ from the PI $\log N_{\mathrm{d}}$ distribution. Even though the PD $N_{\mathrm{d}}$–$\mathcal{L}$ correlation is negative, ${\mathcal{L}}_{\text{PD}}>{\mathcal{L}}_{\text{PI}}$.

[Figure omitted. See PDF]

Figure 6

PD (solid) and PI (dashed) $N_{\mathrm{d}}$ and $\mathcal{L}$ marginal probability distributions and $N_{\mathrm{d}}$–$\mathcal{L}$ conditional probability distributions in two GCMs with unrelated turbulence schemes. The $N_{\mathrm{d}}$–$\mathcal{L}$ relationship based on internal variability within one climate state is not universal across the states.

[Figure omitted. See PDF]

This contradiction raises three questions. First, what produces the noncausal negative $N_{\mathrm{d}}$–$\mathcal{L}$ correlation? We provide a few hypotheses in the following section. Second, considering that these models can replicate the observed PD correlation, what can we infer about the causality of the relationship in observations, where we are unable to conduct direct experimental tests of causality? We discuss this question in Sect. . Third, is any part of the negative relationship between $N_{\mathrm{d}}$ and $\mathcal{L}$ in the models causal? Any such causal mechanism would have to involve a direct or indirect $N_{\mathrm{d}}$ dependence in cloud-top entrainment. In ModelE, the turbulence scheme provides an explicit entrainment closure. have shown that the combination of the Cloud Layers Unified By Binormals (CLUBB; ) cloud and turbulence scheme and the Morrison–Gettelman microphysics scheme can reproduce entrainment-mediated enhanced evaporation at high $N_{\mathrm{d}}$ in single-column experiments. This behavior has not been documented in three-dimensional GCM experiments, but CAM6 and EAMv1 use related cloud–turbulence and cloud–microphysics schemes, so it is conceivable that $N_{\mathrm{d}}$-dependent entrainment mechanisms contribute to the $N_{\mathrm{d}}$–$\mathcal{L}$ relationship in these three models. A deeper investigation of this question merits a separate paper (Part 2 of this series, ).

3.3

Sources of covariability that produce noncausal $N_{\mathrm{d}}$–$\mathcal{L}$ relationships

Noncausal relationships between two variables often originate from a third (possibly unobserved) variable that exerts a causal relationship on the two variables being correlated. This third variable is termed a confounding variable . In its most striking form, confounding can lead to a sign reversal between causation and correlation – for example, in Simpson's paradox . Cloud properties respond strongly to the circulation at the scales of the Sc cellular organization (mesoscale) and greater. Thus, the mesoscale-to-synoptic-scale circulation is a natural place to look for confounding variables that lead to noncausal correlations between cloud properties.

3.3.1 Mesoscale cloud regimes

Mesoscale circulation manifests as cloud regimes (e.g., ). ACI mechanisms likely differ between cloud regimes (e.g., ). This could result in different $N_{\mathrm{d}}$–$\mathcal{L}$ slopes in open- or closed-cell Sc or shallow cumulus or, as the positively and negatively sloped legs of the inverted-V relationship perhaps show, in precipitating and non-precipitating cloud regimes. Due to GCMs' coarse resolution, it is doubtful that they can correctly represent these mesoscale cloud regimes, their ACI mechanisms, or their coupling to the circulation. Nevertheless – or perhaps precisely because we can probably discount cloud-scale causal links between $N_{\mathrm{d}}$ and $\mathcal{L}$ due to the mismatch with the GCM resolved scale – we can use GCMs to test whether the existence of cloud regimes is, on its own, a confounding mechanism for the $N_{\mathrm{d}}$–$\mathcal{L}$ relationship.

To assess whether regime-induced confounding effects may exist in the model $N_{\mathrm{d}}$–$\mathcal{L}$ relationship, we stratify the E3SM model clouds by surface rain rate, $R$. These bins of rain rate are our stand-ins for precipitation regimes. We focus on the surface rain rate because, unlike mesoscale morphological regime definitions (which are subgrid-scale in the GCM), the precipitating–non-precipitating regime delineation has a somewhat clear analog in the GCM. Because the model rain rate has a very long low tail, we do not attempt to define a binary non-precipitating versus precipitating categorization but rather divide the cloud sample into quantiles of rain rate. Specifically, we use sextiles, balancing the need for a meaningful range of rain rates with the need to maintain a large sample of clouds within each bin. The CloudSat precipitation detection sensitivity at the GCM spatial resolution ($R\approx \mathrm{0.01}$ mm d${}^{-\mathrm{1}}$; ) falls roughly into the third rain-rate bin; so, by this definition, half the bins approximately represent precipitating and half non-precipitating clouds.

Figure  shows the results. The model, perhaps unrealistically, produces clouds that generate surface-reaching rain at all droplet concentrations; however, in bins with higher rain rates, the $N_{\mathrm{d}}$ distribution peak is shifted lower, as might be expected from the negative-exponent power law that parameterizes the autoconversion of cloud water to rain, and as is expected from observations . At the same time, $\mathcal{L}$ is higher in bins with higher rain rates, again, as might be expected from the parameterized autoconversion and accretion. Superimposing the bin mean $N_{\mathrm{d}}$ and $\mathcal{L}$ for each rain-rate bin on the unbinned $N_{\mathrm{d}}$–$\mathcal{L}$ distribution, we find that the negative correlation among the bin means echoes the unbinned correlation. This is the case even though, in very classic fashion, the correlations among five out of the six $R$ bins are positive. Comparing a monovariate linear regression of the form $\log \mathcal{L}={\mathit{\alpha }}_{\mathrm{0}}+{\mathit{\alpha }}_{N_{\mathrm{d}}}\log N_{\mathrm{d}}$ (with ${\mathit{\alpha }}_{\mathrm{0}}$ and ${\mathit{\alpha }}_{N_{\mathrm{d}}}$ being the intercept and slope parameters) with a bivariate regression of the form $\log \mathcal{L}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{N_{\mathrm{d}}}\log N_{\mathrm{d}}+{\mathit{\beta }}_R\log R$ (with ${\mathit{\beta }}_{\mathrm{0}}$, ${\mathit{\beta }}_{N_{\mathrm{d}}}$, and ${\mathit{\beta }}_R$ the intercept and slope parameters) provides further confirmation of the Simpson's paradox-like behavior: while the monovariate regression slope is negative (${\mathit{\alpha }}_{N_{\mathrm{d}}}=-\mathrm{0.27}\pm \mathrm{0.001}$; fraction of explained variance is 3.3 %), the bivariate slope is positive (${\mathit{\beta }}_{N_{\mathrm{d}}}=+\mathrm{0.13}\pm \mathrm{0.001}$; fraction of explained variance is 53 %). Both regressions are performed over the range where the $N_{\mathrm{d}}$–$\mathcal{L}$ correlation is negative ($N_{\mathrm{d}}>\mathrm{20}$ cm${}^{-\mathrm{3}}$). The fraction of the sample where $R=\mathrm{0}$ is approximately ${\mathrm{10}}^{-\mathrm{4}}$; these data points are discarded in both regressions. The numbers quoted are ordinary least-squares regression slopes and their standard errors. Thus, the opposing influences of $N_{\mathrm{d}}$ and $\mathcal{L}$ on rain rate can, without any involvement of entrainment or evaporation mechanisms, generate a noncausal negative correlation between $N_{\mathrm{d}}$ and $\mathcal{L}$.

Figure 7

Precipitation-stratified $N_{\mathrm{d}}$ and $\mathcal{L}$ marginal probability distributions and $N_{\mathrm{d}}$–$\mathcal{L}$ conditional probability distributions (colored lines; surface rain-rate intervals given in brackets) in E3SM. The dashed black line shows the unstratified $N_{\mathrm{d}}$–$\mathcal{L}$ relationship. The solid black line connects the mean $(N_{\mathrm{d}},\mathcal{L})$ in each precipitation sextile (colored dots). Binning by precipitation intensity exposes a precipitation-mediated negative $N_{\mathrm{d}}$–$\mathcal{L}$ covariability, with a much steeper slope than the overall $N_{\mathrm{d}}$–$\mathcal{L}$ correlation, even though the $N_{\mathrm{d}}$–$\mathcal{L}$ correlation within all but the least-precipitating sextile is positive.

[Figure omitted. See PDF]

We note that the mechanism generating this noncausal correlation is unusual. The contrast with the causal precipitation suppression mechanism is clear: there, causation runs from $N_{\mathrm{d}}$ to autoconversion to $\mathcal{L}$. Here, causation runs jointly from both $N_{\mathrm{d}}$ and $\mathcal{L}$ to precipitation. How or whether the causal chain then returns from precipitation to $N_{\mathrm{d}}$ and $\mathcal{L}$, as in the classic confounding mechanism, is an open question.

We further note that precipitation already appears to have a qualitative effect on the model's $N_{\mathrm{d}}$–$\mathcal{L}$ relationship at rain rates far below the CloudSat sensitivity threshold: even in the second-lowest $R$ bin, the correlation between $N_{\mathrm{d}}$ and $\mathcal{L}$ is already positive. This suggests that the parameterized precipitation may exert such a strong influence on ACI even for clouds with a low precipitation rate that other ACI adjustment mechanisms, while they may in principle be represented in the model, could be so overwhelmed by the parameterized precipitation suppression that their effect is not discernible in the climate response.

3.3.2 Synoptic-scale air mass advection

At the synoptic to planetary scales, covariability between cloud and aerosol properties can lead to spurious correlations in ACI metrics . Synoptic-scale meteorological covariability can take the form of continental versus marine air mass advection. When an air mass originates over land, it typically has a higher temperature, lower relative humidity (contributing to lower $\mathcal{L}$), and higher aerosol concentration (contributing to higher $N_{\mathrm{d}}$) than when an air mass originates over ocean. This contrast between air masses creates an anticorrelation between $N_{\mathrm{d}}$ and $\mathcal{L}$ even in the absence of any causal effect of $N_{\mathrm{d}}$ on $\mathcal{L}$ . Additionally, sea-surface temperature is coldest and climatological subsidence strongest near the coast, resulting in shallow marine boundary layers. The model's conception of this synoptic-scale covariability in space can be seen in Fig. , with shallow boundary layers and high cloud condensation nuclei (CCN) concentrations near shore and deeper boundary layers with a low CCN farther offshore. A similar covariability exists at particular locations in time. Figure  illustrates the mechanism in the NEP Sc region: at any given location, PBL depth and CCN concentration are strongly linked via the synoptic-scale circulation. Presumably, the position of the anticyclonic subtropical subsidence governs both the PBL depth and whether continental or maritime air is advected.

Figure 8

Within the Sc regime, synoptic-scale meteorology results in strong spatial covariability of (a) PBL depth $h$ and (b) CCN concentration at 0.2 % supersaturation averaged over the depth of the PBL in E3SM. Both of these variables are functions of air mass continentality.

[Figure omitted. See PDF]

Figure 9

Within the Sc regime, synoptic-scale meteorology results in strong temporal covariability between CCN and PBL depth, $h$, in E3SM. This is exemplified by the NEP; the Southern California Bight is depicted in top right corner. The star indicates the grid point with the highest occurrence fraction of Sc conditions according to the criteria of .

[Figure omitted. See PDF]

To assess the synoptic meteorological confounding effect, we stratify the E3SM model clouds in the NEP Sc region by PBL depth. We choose PBL depth as the confounding variable because it appears to act as a proxy for air mass continentality in the model (Fig. ), without a direct parameterized relationship to either aerosols or cloud. PBL depth is nevertheless strongly correlated with both CCN concentration (temporal- and regional-mean vertical profiles are shown in Fig. ) and $\mathcal{L}$. For a fairly wide range of PBL depths (representing the central 90 % of the PBL depth distribution for Sc-regime cloud columns; 10–15 model levels, corresponding approximately to 750–1400 m), the relationship between mean $N_{\mathrm{d}}$ and mean $\mathcal{L}$, stratified by PBL depth, mimics the slope of the unstratified $N_{\mathrm{d}}$–$\mathcal{L}$ relationship quite closely (Fig. ). Based on this, it is plausible that synoptic-scale meteorological covariability contributes substantially to the overall negative $N_{\mathrm{d}}$–$\mathcal{L}$ correlation in the model.

Figure 10

Temporal- and regional-mean CCN concentration profiles in the NEP Sc region stratified by PBL depth in E3SM. Within the Sc regime, CCN is strongly sorted by PBL depth, illustrating the strong covariability between PBL thermodynamic structure and aerosol advection. The central 90 % of the PBL depth range (between 10 and 15 model levels, corresponding to approximately 750–1400 m) is shown to avoid outliers in the low-statistics PBL depth bins. The PBL depth is measured in units of model levels $k$, with $k$ increasing downward from the level of the PBL-capping inversion $k_{\text{pbl}}$ to the model level closest to the surface ($k=\mathrm{72}$ in EAMv1).

[Figure omitted. See PDF]

Figure 11

Within the Sc regime, PBL depth–CCN covariability leads to a negative $N_{\mathrm{d}}$–$\mathcal{L}$ correlation with a slope similar to the overall $N_{\mathrm{d}}$–$\mathcal{L}$ correlation in E3SM. The dashed black line shows the $N_{\mathrm{d}}$–$\mathcal{L}$ relationship not stratified by PBL depth. The solid black line connects the mean $(N_{\mathrm{d}},\mathcal{L})$ at each PBL depth (colored dots). The central 90 % of the PBL depth range (between 10 and 15 model levels, corresponding to approximately 750–1400 m) is shown to avoid outliers in the low-statistics PBL depth bins.

[Figure omitted. See PDF]

We note several caveats. The synoptic-scale covariability of aerosol advection and PBL depth is a feature of the general circulation and can therefore be expected to be modeled reasonably well in a GCM. However, the interaction of the advected aerosol with boundary-layer clouds depends on mixing between the free troposphere and the boundary layer, which is likely much less well represented in GCMs that have coarse vertical resolutions. Furthermore, cloud adjustments mediated by circulation changes are artificially suppressed as a side effect of the fixed sea-surface temperature and nudged winds that we employ to reduce internal-variability noise in the ERFaci estimates. Whether the synoptic-scale confounding signature in the model mimics the real climate is therefore uncertain. Further, the synoptic-scale covariability differs depending on the geographic particulars of each Sc basin; we have only analyzed the NEP Sc in detail. Finally, while the PBL depth-stratified negative $N_{\mathrm{d}}$–$\mathcal{L}$ relationship in the model is consistent with observational analyses (e.g., ), the model does not reproduce the weakening of the $N_{\mathrm{d}}$–$\mathcal{L}$ correlation within each PBL depth bin that is found in observations . We compare (a) a monovariate linear regression of the form $\log \mathcal{L}={\mathit{\alpha }}_{\mathrm{0}}+{\mathit{\alpha }}_{N_{\mathrm{d}}}\log N_{\mathrm{d}}$ (with ${\mathit{\alpha }}_{\mathrm{0}}$ and ${\mathit{\alpha }}_{N_{\mathrm{d}}}$ being the intercept and slope parameters) over the unstratified data, (b) a monovariate linear regression of the form $\log \mathcal{L}={\mathit{\alpha }}_{\mathrm{0}}+{\mathit{\alpha }}_{N_{\mathrm{d}}|h}\log N_{\mathrm{d}}$ (with ${\mathit{\alpha }}_{\mathrm{0}}$ and ${\mathit{\alpha }}_{N_{\mathrm{d}}|h}$ the intercept and slope parameters) stratified by PBL depth, and (c) a bivariate regression of the form $\log \mathcal{L}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{N_{\mathrm{d}}}\log N_{\mathrm{d}}+{\mathit{\beta }}_R\log h$ (with ${\mathit{\beta }}_{\mathrm{0}}$, ${\mathit{\beta }}_{N_{\mathrm{d}}}$, and ${\mathit{\beta }}_h$ the intercept and slope parameters). In all three cases, the slope parameters (${\mathit{\alpha }}_{N_{\mathrm{d}}}$, ${\mathit{\alpha }}_{N_{\mathrm{d}}|h}$, and ${\mathit{\beta }}_{N_{\mathrm{d}}}$) lie in a range from $-\mathrm{0.33}$ to $-\mathrm{0.30}$. The ordinary least-squares regressions are performed over the $N_{\mathrm{d}}$ range where the $N_{\mathrm{d}}$–$\mathcal{L}$ correlation is negative ($N_{\mathrm{d}}>\mathrm{20}$ cm${}^{-\mathrm{3}}$) and over the central 90 % PBL depth range (10–15 model levels, corresponding approximately to 750–1400 m).

3.3.3 Phase-space boundaries

Correlations between $N_{\mathrm{d}}$ and $\mathcal{L}$ can also arise simply because not all parts of the $N_{\mathrm{d}}$–$\mathcal{L}$ phase space are uniformly populated by clouds. This can be illustrated by applying the MODIS simulator to the model. The MODIS simulator provides the optical thickness, $\mathit{\tau }$, and droplet effective radius, $r_{\mathrm{e}}$, diagnosed consistently with the MODIS satellite cloud retrievals . Power-law adiabatic relationships $\mathcal{L}(r_{\mathrm{e}},\mathit{\tau })$ and $N_{\mathrm{d}}(r_{\mathrm{e}},\mathit{\tau })$ can be used to transform the MODIS output into $N_{\mathrm{d}}$–$\mathcal{L}$ space (e.g., ). In logarithmic coordinates, this is a linear transformation, yielding the correlation shown in Fig. . This S-curve correlation, like the model-native $N_{\mathrm{d}}$–$\mathcal{L}$ correlation, shows a steep rise in $\mathcal{L}$ at low $N_{\mathrm{d}}$ and a steep drop at moderate $N_{\mathrm{d}}$. It also shows another steep rise at high $N_{\mathrm{d}}$ that the model may hint at but does not exhibit clearly. Investigating the data before and after the coordinate transformation to $N_{\mathrm{d}}$–$\mathcal{L}$ space is instructive. In $\log \mathit{\tau }$–$\log r_{\mathrm{e}}$ space, the MODIS simulator output falls within a rectangle, bounded by the limits the model prescribes on its clouds. Upon transformation to $\log \mathcal{L}$–$\log N_{\mathrm{d}}$ space, the population is bounded by a parallelogram (see the isolines in Fig. ). These limits on the phase space strongly sculpt the behavior of the mean $\log \mathcal{L}$ as a function of $\log N_{\mathrm{d}}$ because the parts of phase space that are not populated do not contribute to the mean $\mathcal{L}$ as a function of $N_{\mathrm{d}}$.

Figure 12

Joint probability distribution of E3SM MODIS-simulated $\mathcal{L}$ and $N_{\mathrm{d}}$. Isolines of $r_{\mathrm{e}}$ and $\mathit{\tau }$, from which adiabatic $N_{\mathrm{d}}$ and $\mathcal{L}$ are retrieved, are overlaid. The mean $\log \mathcal{L}$ as a function of $\log N_{\mathrm{d}}$ is shown as a blue line. Because the model imposes a rectangular limit on $r_{\mathrm{e}}$ and $\mathit{\tau }$, the $N_{\mathrm{d}}$–$\mathcal{L}$ phase space has parallelogram-shaped boundaries. At least part of the rise, fall, and repeated rise in $\mathcal{L}$ as a function of $N_{\mathrm{d}}$ (blue line) is due to these phase-space boundaries cutting off the upper and lower parts of the $\mathcal{L}$ distribution.

[Figure omitted. See PDF]

3.4 Persistent disagreement with other lines of evidence

Before these results, it was only logical to discount the GCM evidence on the basis that it could not reproduce the observed $N_{\mathrm{d}}$–$\mathcal{L}$ relationship in PD internal variability. Now that some GCMs match the other lines of evidence in PD internal variability, what do we make of the fact that the disagreement on the sign of the causal climatic $\mathcal{L}$ adjustment to RFaci persists?

In observations, it is more difficult to establish causality than in the GCMs, where it is as simple as changing the aerosol emissions while fixing all other boundary conditions. The most reliable causal evidence in observations comes from observational natural laboratories, where the aerosol perturbation is known and an unperturbed control can be identified clearly . Such laboratories indicate unchanged or reduced $\mathcal{L}$ in the perturbed clouds . But, such laboratories are rare, and there is no rigorous extrapolation from these laboratories to the full diversity of cloud regimes found in the climate. The most representative observations – that is, the global satellite-retrieved inverted-V correlations – have the opposite problem: they are representative, but are the correlations causal? The correlation is more negative than the estimate of the causal interannual $\mathcal{L}$ response to $N_{\mathrm{d}}$ perturbations using an effusive volcano as a laboratory (; albeit for shallow Cu rather than Sc). argue that satellite $N_{\mathrm{d}}$–$\mathcal{L}$ correlations are negatively biased not only by covariability confounding but also by retrieval errors. applied a causal network approach to the temporal evolution in geostationary satellite data and found that the causal negative $N_{\mathrm{d}}$–$\mathcal{L}$ relationship is weaker than the $N_{\mathrm{d}}$–$\mathcal{L}$ correlation. Strong regional increasing and declining trends on multidecadal timescales in the satellite record may also contribute to disentangling covariability and causality .

In LES, as in GCMs, causality can be established by varying aerosol concentration, while keeping the other boundary conditions constant. This provides very clear evidence that precipitation suppression and entrainment feedbacks lead to process fingerprints of positive and negative $\mathcal{L}$ tendencies in $N_{\mathrm{d}}$–$\mathcal{L}$ space that translate into steady Sc states . But these LES experiments are expensive, so boundary conditions are carefully curated to a very small subset of the high-dimensional space of meteorological conditions present in the climate. We simply do not know whether the process fingerprints would be as unambiguous if a broader spectrum of boundary conditions were simulated or if the clouds were able to interact with larger scales in the multiscale climate problem instead of evolving to a steady state.

In summary, GCMs are still the odd ones out in their negative $\mathcal{L}$ adjustment component of ERFaci. The observational and LES modeling lines of evidence have clear confounding and representativeness problems. Are these problems severe enough to flip the sign of the adjustment? It seems unlikely, but our GCM results show that it is possible; addressing the representativeness and confounder questions in the other lines of evidence thus takes on a renewed urgency.

4 Conclusions

expressed the hope that global models, after a long stretch of playing the odd line of evidence out in assessments of global energy budget problems , might return to a more equal role in the balance and struggle between conflicting lines of evidence. One way in which the global model perspective shores up the strength of the multiline assessment by providing information not available from the other lines of evidence: being able to test causality and showing that PD internal variability may not even correctly predict the sign of the causal cloud water adjustment to the anthropogenic cloud droplet perturbation.

Causality (or, in this case, lack of causality) is easy to establish in a model experiment but very difficult in observations. Where the noncausal correlation originates is another question that models can, in principle, answer definitively by shutting off confounding model mechanisms in mechanism-denial experiments. In Part 2 of this series, we more fully use the power of models as hypothesis testers by performing perturbed-physics and mechanism-denial experiments. In this paper, we have restricted ourselves to slice-and-dice analyses that could, in principle, also be performed on observations. We hope that they will be performed on observations, especially if Lagrangian investigation of the cloud life cycle (e.g., ) and observational fingerprints of loss processes (e.g., ) can be included.

Whether lack of causality in the model system implies lack of causality in the real atmosphere is a question that models alone cannot address, so we do not yet know how worried we need to be about the sign difference between correlation and causation in the model $N_{\mathrm{d}}$–$\mathcal{L}$ relationship. When it comes to the non-GCM lines of evidence, one can quibble with the representativeness of the causal evidence and with causality in the representative evidence – at the very least, these model results are a flashing-red warning sign hanging over our interpretation of the $\mathcal{L}$ adjustment component of ERFaci.

Code and data availability

Model output for the AeroCom IND3 experiment is available (subject to acknowledgment of a data policy at https://aerocom.met.no/FAQ/data_access, last access: 26 June 2024) from ssh://aerocom-users.met.no/metno/aerocom-users-database/AEROCOM-PHASE-II-IND3 . The US CMS model output is available at https://doi.org/10.5281/zenodo.10449670 . The analysis code is available at https://doi.org/10.5281/zenodo.10449750 .

Author contributions

The idea for this research came from ASA, SEB, MWC, AMF, SD, EG, LRL, PLM, YM, JM, JQ, ACV, HW, KZ, and YZ.

ASA, AG, XL, JM, DN, DGP, PS, TT, and YZ performed model experiments.

ASA, AMF, SD, AG, EG, PLM, JM, JQ, FT, and YZ contributed to analysis.

The original manuscript draft was prepared by JM, and all co-authors provided comments.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Atmospheric Chemistry and Physics. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We thank the two reviewers, Anna Possner and Jianhao Zhang, for suggestions that have resulted in improvements in the analysis and presentation. We thank Christopher Bretherton, Susannah Burrows, Yao-Sheng Chen, Leo Donner, Graham Feingold, Colleen Kaul, Jan Kazil, Naser Mahfouz, Daniel McCoy, Isabel McCoy, Christina Sackmann, Rob Wood, Jianhao Zhang, and Xiaoli Zhou for comments and discussion. This work arises from the 2020 US Climate Modeling Summit held virtually and chaired by Gavin Schmidt and Susanne Bauer.

Financial support

Johannes Mülmenstädt was supported by the Office of Science, US Department of Energy (DOE) Biological and Environmental Research as part of the Regional and Global Model Analysis program area and used resources of the National Energy Research Scientific Computing Center (NERSC), a US DOE Office of Science user facility located at Lawrence Berkeley National Laboratory (contract no. DE-AC02-05CH11231). Andrew S. Ackerman, Ann M. Fridlind, Florian Tornow, and Susanne E. Bauer were supported by the NASA Modeling, Analysis, and Prediction Program; their computational resources were provided by the NASA Center for Climate Simulation (NCCS) at Goddard Space Flight Center. Edward Gryspeerdt was supported by a Royal Society University Research Fellowship (URF$\backslash$R1$\backslash$191602). David Neubauer is supported by the European Union’s Horizon 2020 research and innovation program (grant no. 821205). Toshihiko Takemura was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI grant (grant no. JP19H05669) and the Environment Research and Technology Development Fund S-20 (grant no. JPMEERF21S12010) of the Environmental Restoration and Conservation Agency provided by the Ministry of the Environment, Japan. Philip Stier was supported by the European Research Council project RECAP under the European Union’s Horizon 2020 research and innovation program (grant no. 724602) and the FORCeS project under the European Union’s Horizon 2020 research and innovation program (grant no. 821205). The Pacific Northwest National Laboratory (PNNL) is operated for DOE by Battelle Memorial Institute (contract no. DE-AC05-76RLO1830).

Review statement

This paper was edited by Matthew Lebsock and reviewed by Anna Possner and Jianhao Zhang.