Untangling causality in midlatitude aerosol–cloud adjustments

Aerosol–cloud interactions represent the leading uncertainty in our ability to infer climate sensitivity from the observational record. The forcing from changes in cloud albedo driven by increases in cloud droplet number (N_d) (the first indirect effect) is confidently negative and has narrowed its probable range in the last decade, but the sign and strength of forcing associated with changes in cloud macrophysics in response to aerosol (aerosol–cloud adjustments) remain uncertain. This uncertainty reflects our inability to accurately quantify variability not associated with a causal link flowing from the cloud microphysical state to the cloud macrophysical state. Once variability associated with meteorology has been removed, covariance between the liquid water path (LWP) averaged across cloudy and clear regions (here characterizing the macrophysical state) and N_d (characterizing the microphysical) is the sum of two causal pathways linking N_d to LWP: N_d altering LWP (adjustments) and precipitation scavenging aerosol and thus depleting N_d. Only the former term is relevant to constraining adjustments, but disentangling these terms in observations is challenging. We hypothesize that the diversity of constraints on aerosol–cloud adjustments in the literature may be partly due to not explicitly characterizing covariance flowing from cloud to aerosol and aerosol to cloud. Here, we restrict our analysis to the regime of extratropical clouds outside of low-pressure centers associated with cyclonic activity. Observations from MAC-LWP (Multisensor Advanced Climatology of Liquid Water Path) and MODIS are compared to simulations in the Met Office Unified Model (UM) GA7.1 (the atmosphere model of HadGEM3-GC3.1 and UKESM1). The meteorological predictors of LWP are found to be similar between the model and observations. There is also agreement with previous literature on cloud-controlling factors finding that increasing stability, moisture, and sensible heat flux enhance LWP, while increasing subsidence and sea surface temperature decrease it. A simulation where cloud microphysics are insensitive to changes in N_d is used to characterize covariance between N_d and LWP that is induced by factors other than aerosol–cloud adjustments. By removing variability associated with meteorology and scavenging, we infer the sensitivity of LWP to changes in N_d. Application of this technique to UM GA7.1 simulations reproduces the true model adjustment strength. Observational constraints developed using simulated covariability not induced by adjustments and observed covariability between N_d and LWP predict a 25 %–30 % overestimate by the UM GA7.1 in LWP change and a 30 %–35 % overestimate in associated radiative forcing.

1 Introduction

Uncertainty in the radiative forcing due to aerosol–cloud interactions is the leading uncertainty limiting our ability to accurately diagnose the Earth's climate sensitivity from the observational record (Forster, 2016). The best estimate of the radiative forcing due to aerosol–cloud interactions (also called the first indirect effect; Twomey, 1977) has narrowed to −1.2 to −0.34 W m⁻² in a recent survey of forcing from aerosol–cloud interactions (Storelvmo, 2017; Bellouin et al., 2020), but the sign and strength of the forcing due to changes in cloud macrophysical properties in response to aerosol (aerosol–cloud adjustments) remain uncertain (Bellouin et al., 2020). This uncertainty reflects the difficulty in disentangling the many factors that determine cloud macrophysical properties. Unlike cloud droplet number concentration (N_d), which is primarily driven by the availability of suitable aerosol, cloud macrophysical properties are primarily determined by the state of the atmosphere but may be modulated by N_d (Stevens and Feingold, 2009).

Here, we focus on liquid clouds. Within liquid clouds two main processes are hypothesized to alter cloud liquid content in response to changes in N_d. As N_d is increased the following may occur: (1) a suppression of precipitation by enhanced N_d (Albrecht, 1989; Pincus and Baker, 1994) and (2) strengthened entrainment of dry air (Ackerman et al., 2004; Bretherton et al., 2007; Wang et al., 2003; Xue and Feingold, 2006). These processes drive the liquid content of cloud in opposite directions following a perturbation in microphysics, further complicating the interpretation of covariability between cloud macrophysical and microphysical properties. Because these processes exist at a timescale and length scale far below those resolved in climate models, they must be parameterized, resulting in substantial uncertainty. Because of this, constraining the cloud macrophysical response to changes in cloud microphysics using observations is an essential step towards constraining aerosol forcing. Aerosol–cloud adjustments are likely to contribute a relatively small fraction of overall variability in cloud macrophysical properties (Bender et al., 2019) – making the data volume available from remote-sensing observations particularly relevant to providing an observational constraint by detecting the small signal from aerosol–cloud adjustments hidden in meteorologically induced cloud variability. The literature has produced numerous careful analyses of the observational record, but these analyses have produced divergent estimates of the aggregate effect of increased N_d on liquid clouds. These impacts range from increased N_d increasing liquid content in clouds (Chen et al., 2014; Rosenfeld et al., 2019) to almost no response in in-cloud and area-averaged liquid water path (LWP; Toll et al., 2017; Malavelle et al., 2017) and to decreasing liquid content (Gryspeerdt et al., 2019; Sato et al., 2018; Toll et al., 2019). In contrast, the response in cloud cover to N_d has tended to be consistently inferred as producing a forcing equal to or larger than the first indirect effect (Gryspeerdt et al., 2016; Christensen et al., 2017; Andersen et al., 2017). A simulation of deep, open-cellular boundary layers has shown an increase in the cloud fraction and a decrease in in-cloud LWP, ultimately resulting in an increase in reflected shortwave in the simulation (Possner et al., 2018). Given the wide range of potential observational constraints it remains difficult to offer advice on the cloud microphysical parameterizations that are the most realistic.

However, the causal link flowing from N_d to clouds is not the only mechanism relating cloud macrophysics to cloud microphysics. As shown in Wood et al. (2012), the spatial pattern in N_d observed in nature is primarily determined by precipitation scavenging. That is to say, N_d may both alter cloud macrophysical properties and be altered by them, and interpreting observations as a constraint on only microphysical to macrophysical causality is erroneous, as discussed in Gryspeerdt et al. (2019). To be able to constrain cloud responses to aerosol-driven changes in N_d we must be able to characterize the effects of cloud macrophysical properties on N_d (by precipitation scavenging) and covariability between N_d and LWP induced by confounding factors such as the relative location of aerosol sources and climatological cloud cover or meteorological modulation of aerosol and clouds.

Bender et al. (2019) proposed utilizing the observed covariance between N_d and LWP as an aggregate measure of aerosol–cloud-radiation behavior across climate models and reanalyses. For a model to be able to realistically reproduce aerosol–cloud behavior one precondition is that covariance flowing from clouds to aerosol and from aerosol to clouds as well as the non-causal covariance induced by air mass history must result in total covariance consistent with observations. Similarly, Gryspeerdt et al. (2019) examined covariance between N_d and in-cloud LWP but proposed to separate variability induced by adjustments from other sources of covariability by using clustering on global data. This study follows these earlier studies and attempts to partition covariance between LWP and N_d related to causality flowing from N_d to LWP from other covariance. We utilize empirical analysis of observations and model output to try to disentangle variability driven by meteorological “cloud-controlling factors” (Stevens and Brenguier, 2009) and variability related to changes in cloud microphysics. Within this framework we focus on the midlatitudes. Previous work has performed similar analysis on cyclonic midlatitude systems, finding a distinct increase in LWP at a fixed precipitation rate (McCoy et al., 2018b). These large synoptic systems account for roughly half of the midlatitudes (Bodas-Salcedo et al., 2014). Here, we turn our attention to the remaining interstitial regions between cyclones. This regime tends to be less cloudy, have an overall lower albedo, and have a lower cloud optical depth (an example is shown in Fig. 1a). Thus, while these regions do not host dramatic fronts and cloud shields, changes in their liquid cloud due to anthropogenic aerosol may contribute strongly to the overall aerosol forcing because cloud areal coverage and cloud optical depth are further below saturation.

Figure 1 (a) An example of a cyclone off the coast of North America from MODIS Aqua. The cyclone center has been identified south of the Aleut Peninsula based on MERRA-2 sea-level pressure (SLP). A line shows the edge of the area considered to lie within the cyclone following Field and Wood (2007). (b) The fraction of data outside of cyclones averaged between 2003 and 2015. The mean within the 30–60^∘ N region is 42 %.

2 Methods

2.1 Region and variable selection

The empirical analysis on both observational and simulated data in this work follows pioneering studies investigating meteorological controls on subtropical (Myers and Norris, 2013, 2015; Qu et al., 2015; Klein et al., 2017; Seethala et al., 2015) and midlatitude (Miyamoto et al., 2018) marine boundary-layer cloudiness. The selection of variables used in this study is based on the findings of these earlier studies. The predictors examined here are sea surface temperature (SST), large-scale subsidence (ω) at 550 hPa, the strength of the inversion at the top of the boundary layer (estimated inversion strength – EIS – as defined in Wood and Bretherton, 2006), the total water vapor path (WVP), and sensible heat flux (SHF). A summary is given in Table 1. A priori we expect that increasing SSTs should lead to decreasing cloudiness (Myers and Norris, 2015; Qu et al., 2015; McCoy et al., 2017; Bretherton and Blossey, 2014). Increasing subsidence should decrease the liquid water path (Myers and Norris, 2013). Increasing inversion strength should lead to an increasing cloud fraction (Wood and Bretherton, 2006; Klein and Hartmann, 1993; Qu et al., 2015; Myers and Norris, 2015; Seethala et al., 2015). As shown in Miyamoto et al. (2018) from examination of midlatitude ocean fronts, increased SHF should increase cloudiness. Finally, while less extensively studied, an overall increase in WVP seems likely to lead to increased cloudiness as the atmospheric column moisture increases, holding all else equal (McCoy et al., 2019).

Table 1 List of observed variables.

The central goal of this study is to develop an empirical constraint on aerosol–cloud adjustments. However, direct measurements of the aerosol state to constrain aerosol–cloud adjustments have been found to be sensitive to errors in pristine conditions (Ma et al., 2018). To develop a constraint on aerosol–cloud adjustments we examine the cloud microphysical state in addition to the meteorological parameters listed above. Cloud droplet number concentration is the state variable relating cloud macrophysics to aerosol concentration (Wood, 2012; Grosvenor et al., 2018) and is the predictor used in this study to characterize aerosol–cloud adjustments.

We hasten to note that the list of predictors described above is not intended to be complete (that is to say, we do not anticipate that we will be able to explain 100 % of the variability in clouds with these predictors). Our goal is to explain enough variability with these predictors to be able to infer statistically robust relationships between cloud macrophysical and microphysical properties in the variance unexplained by meteorological variability. The predictive ability of these inferred relationships will be tested within the context of a global climate model (GCM).

Having discussed predictors, we will now discuss predictands. Outside of cyclones cloud cover is primarily liquid, boundary-layer cloud, which substantially affects shortwave radiation but has little effect on longwave radiation (Hartmann and Short, 1980). Thus, we focus on the liquid water path (LWP) averaged over cloudy and clear regions. This variable allows observations and model output to be compared without needing to simulate output from a passive spectroradiometer (Bodas-Salcedo et al., 2011) or considering cloud overlap, as would be required for in-cloud LWP.

Analysis is carried out in the Northern Hemisphere (here 30–60^∘ N) over oceans. This choice has been made for several reasons: first, there is a large variability in N_d across the northern midlatitudes that is not strongly driven by the seasonal cycle. In the Southern Ocean cloud condensation nuclei (CCN) variability is primarily due to biogenic sources, and their variability is strongly driven by the seasonal cycle, making analysis in the context of meteorological controls difficult to interpret (Ayers and Gras, 1991; McCoy et al., 2015a; Charlson et al., 1987). Second, Southern Ocean N_d is still poorly represented in the Met Office Unified Model GA7.1 (discussed in Sect. 2.4), with both too low a mean value and generally too small a seasonal cycle – in contrast to the Northern Hemisphere, which compares well with observations (Mulcahy et al., 2018).

2.2 Outside of cyclone compositing

The passage of synoptic systems is the central mode of variability in midlatitude cloudiness. McCoy et al. (2018b) focused on cyclone systems and showed that once meteorological variability was accounted for, the effects of changes in cloud microphysics driven by changes in N_d on liquid content are revealed. This work utilized the Field and Wood (2007) cyclone compositing algorithm, which uses sea-level pressure (SLP) to identify cyclone centers. Here we utilize the same compositing approach but focus on all data that are 2000 km away from a cyclone center, following Bodas-Salcedo et al. (2014). This separates the midlatitudes into times when a cyclone is nearby and when there is no cyclone nearby. Once cyclone centers are identified all data within 2000 km of an identified cyclone center are masked out of the data set. All data over land are also masked because microwave observations are unavailable (see next section).

A schematic representation of the frequency of occurrence of out-of-cyclone states is shown in Fig. 1b and is 42 % in the 30–60^∘ N region focused on in this study. All data are regridded to a common 1^∘ × 1^∘ spatial resolution before analysis. All data are daily mean, except for observed N_d, which is only available for the MODIS Aqua overpass; however the diurnal cycle of N_d is thought to be relatively slight (Dong et al., 2014), and it is unlikely that this retrieval limitation significantly impacts our results. Our choice of resolution is motivated by the natural resolution of the Multisensor Advanced Climatology of Liquid Water Path (MAC-LWP) data set. Performing the analysis with a coarser resolution might change the regression coefficients derived in Sect. 2.5, as has been reported for analysis of the first indirect effect (McComiskey and Feingold, 2012). However, in Sect. 3.4 we validate the adjustments inferred from this resolution against the true adjustments in Met Office Unified Model (UM) GA7.1 (the atmosphere model of HadGEM3-GC3.1 and UKESM1) and find agreement, which demonstrates that our choice of resolution does not affect our results.

2.3 Observations

The cloud liquid water path (LWP) is calculated based on the aggregated observations from multiple satellite microwave sensors following the methodology of the multisensor advanced climatology of LWP (MAC-LWP; Elsaesser et al., 2017) at a daily mean resolution. Microwave radiometers are sensitive to total liquid. This means that the retrieval must calculate a partitioning between precipitating and non-precipitating liquid. Here, the liquid water path in clouds is calculated following Elsaesser et al. (2017). Previous analysis has shown that the midlatitude partitioning of rain and cloud in MAC-LWP compares favorably with convection-permitting simulations (McCoy et al., 2018b). The MAC-LWP data set estimates that 7 % of total liquid path is rain in the NH midlatitudes outside of cyclones. Because of the low fraction of rainwater in the midlatitudes outside of cyclones it is unlikely that the need to partition the observations into rain and cloud water substantially affects observational constraints as calculated in Sect. 3.3. The liquid water path is defined as the average of cloudy and clear regions and is insensitive to overlying ice cloud (unlike passive optical retrievals). This definition of LWP is consistent with the definition from GCMs. Microwave retrievals are only available over the ocean.

Observations of N_d are calculated based on MODIS observations of cloud optical depth and droplet effective radius. Level 2 swath data (joint product) from MODIS collection 5.1 (King et al., 2003) are filtered by removing pixels with solar zenith angles greater than 65^∘ to eliminate problematic retrievals at a pixel level (Grosvenor and Wood, 2014). The daily mean N_d at 1^∘ × 1^∘ resolution is calculated from the filtered level 2 swath data and only low (cloud tops below 3.2 km), liquid clouds were used to calculate N_d. Only 1^∘ × 1^∘ regions where the liquid cloud fraction exceeds 80 % are considered valid (Bennartz et al., 2011), and the N_d is calculated using the effective radius from the 3.7 µm MODIS channel. This data set is evaluated in Grosvenor et al. (2018) and in McCoy et al. (2018a), where it shows consistency with measurements from aircraft.

We feel that this combination of radiometers in constraining adjustments is particularly advantageous because there is no shared information between the microwave retrieval of LWP and the shortwave and near-infrared retrieval of N_d. This is in contrast to the common practice of utilizing MODIS retrievals of cloud droplet effective radius and optical depth to calculate both the in-cloud LWP and N_d, making it conceptually difficult to cleanly infer changes in LWP due to changes in N_d.

The meteorological state of the atmosphere is characterized using reanalysis from MERRA-2 (Molod et al., 2015) and microwave observations of WVP from MAC-LWP. The period during which all of these observations are available stretches from 2003 to 2015. A list of variables and data sources is given in Table 1.

2.4 Simulations

Simulations are carried out in atmosphere-only simulations in the UM version 10.8 at N96 horizontal resolution. The version of the atmosphere model is GA7.1 coupled to the UKCA (UK Chemistry and Aerosols) model as described in (Mulcahy et al., 2018). The GA7.1 atmosphere model is the physical basis for the atmosphere model included in the HadGEM3-GC3.1 and UKESM1 climate model submission to CMIP6 (Walters et al., 2019). All simulations are run for 2 years and 3 months, starting in September 2013. Simulations are nudged to European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis winds above model level 16 (approximately 1.8 km). Paired simulations are conducted with emissions of aerosol and precursor gases set to preindustrial (henceforth labeled PI) and to present-day (labeled PD) levels. Anthropogenic emissions use the CMIP6 inventory (Eyring et al., 2016), natural emissions over land use the MEGAN-MACC inventory (Sindelarova et al., 2014), and natural emissions over the ocean follow the POET inventory (Granier et al., 2005). These simulations are conducted to calculate the response of cloud to anthropogenic aerosol. Additional sensitivity studies were carried out, starting with nudging at level 27 (approximately 5 km) and without nudging. In the latter case simulations of the PI and PD were run for 10 years and 3 months, starting in September 2003. The results of this analysis are shown in the Supplement and do not differ substantially from the results in the main text.

In GA7.1 aerosol number concentration has the potential to affect cloud microphysics in three places: first, the conversion of cloud water to rain, where enhanced N_d suppresses the conversion of cloud water to rainwater. This is parameterized using the scheme in Khairoutdinov and Kogan (2000), bias corrected as described in Boutle et al. (2014). The conversion rate of cloud to rain is parameterized as

$$\begin{array}{}\text{(1)}& P_{\text{Cloud}\to \text{Rain}}\propto q_{\text{Cloud}}^{\mathrm{2.47}}\cdot N_{\mathrm{d}}^{-\mathrm{1.79}},\end{array}$$

where q_Cloud is the cloud liquid water mixing ratio. Second, it can affect cloud microphysics in the gravitational settling of cloud droplets, where enhanced N_d suppresses the settling of cloud droplets out of the cloud. The flux of water settling out of a cloud layer is given as

$$\begin{array}{}\text{(2)}& P_{\text{Settle}}\propto q_{\text{Cloud}}^{\mathrm{5}/\mathrm{3}}{N}_{\mathrm{d}}^{-\mathrm{2}/\mathrm{3}}.\end{array}$$

Third, it can affect cloud albedo through changes in N_d (Twomey, 1977).

However, clouds may affect aerosol concentrations via wet scavenging. UKCA allows scavenging to be configured in two ways: the default configuration or using the scavenging coefficients in the ECHAM5-HAM model (Stier et al., 2005). To evaluate some of the uncertainty related to scavenging we perform simulations with both sets of coefficients. A list of model simulations conducted is given in Table 2.

Table 2 List of simulations: the Met Office model designation; a short description; the scavenging coefficients used; the value of N_d used in the radiation, settling, and autoconversion parameterizations; and aerosol emissions period used (present day – PD – or preindustrial – PI). A tilde in N_d for a given parameterization indicates that N_d was not set to a constant value.

2.5 Analyzing covariance in the system

Our goal in this paper is to characterize how N_d drives LWP. However, to reveal robust correlations between N_d and LWP, meteorological variability must be accounted for. To achieve this we follow previous empirical analysis of cloud-controlling factors and utilize multiple linear regression to characterize the dependence of various cloud properties as predictands on the predictors described above (Myers and Norris, 2015). As with all empirical analyses of observations, correlation between predictors and the predictand does not necessarily mean that they are causally linked. The meteorological predictors selected here have had mechanisms hypothesized to link them to cloud variability. Another issue with this and similar analyses of the response of cloud to cloud-controlling factors is that thermodynamic and dynamical predictors within the Earth's atmosphere are correlated (McCoy et al., 2017; Myers and Norris, 2015, 2013). As the variance shared by predictors grows, the uncertainty in the coefficients calculated by multiple linear regression increases. While this is an issue, it is inescapable in trying to disentangle the effects of different cloud-controlling factors on cloud properties. Examples of meteorological variables that are correlated but have opposing effects on cloud cover include subsidence and stability and stability and SST (Myers and Norris, 2013; Qu et al., 2015; Wall et al., 2017).

To reduce predictor covariance we bin our 1^∘ × 1^∘ daily mean data from model output and observations into the space of SST and WVP. All output over oceans is considered in the analysis. We choose these predictors to bin our data before performing regression analysis because they represent two of the most basic meteorological state variables, because they covary strongly, and because LWP varies significantly as a function of WVP and SST in both observations and simulations (Fig. 2); 9 bins are used to span WVP (0–40 kg m⁻²), and 10 are used for SST (275–305 K). We propose that a multiple-linear-regression model of the form,

$$\begin{array}{ll}{\displaystyle }\text{LWP}=& {\displaystyle }{a}_{\text{ccn}}\text{CCN}+a_{\mathrm{2}}{\mathit{\omega }}_{\mathrm{550}}+a_{\mathrm{3}}\text{EIS}+a_{\mathrm{4}}\text{WVP}+a_{\mathrm{5}}\text{SST}\\ \text{(3)}& {\displaystyle }& {\displaystyle }+a_{\mathrm{6}}\text{SHF}+a_{\mathrm{7}},\end{array}$$

be trained in each bin of SST and WVP. All the predictors within the regression model are environmental drivers that are external to the clouds being influenced. However, robust remote-sensing observations of CCN are lacking and ln(N_d) is used in its place. As discussed above, N_d is predominantly a function of CCN availability. Thus, a multiple-linear-regression model,

$$\begin{array}{ll}{\displaystyle }\text{LWP}=& {\displaystyle }{a}_{\mathrm{1}}\ln \left(N_{\mathrm{d}}\right)+a_{\mathrm{2}}{\mathit{\omega }}_{\mathrm{550}}+a_{\mathrm{3}}\text{EIS}+a_{\mathrm{4}}\text{WVP}+a_{\mathrm{5}}\text{SST}\\ \text{(4)}& {\displaystyle }& {\displaystyle }+a_{\mathrm{6}}\text{SHF}+a_{\mathrm{7}},\end{array}$$

is trained in each bin of WVP and SST. Bins that contain fewer than 3000 data points are excluded. (Units of LWP are in g m⁻², N_d is in cm⁻³, subsidence is in Pa s⁻¹, EIS is in K, WVP is in kg m⁻², SST is in K, and SHF is in W m⁻².) A 99 % confidence interval is used to determine if a coefficient is significant (different from zero). The initiation of precipitation in boundary-layer cloud can substantially impact LWP within a given cloud with a timescale longer than the daily sampling timescale utilized here (Berner et al., 2013). However, the spatial scale of our analysis (1^∘ × 1^∘) samples clouds at many stages of their life cycle. This means that we do not need to consider rain to be a predictor of LWP and may think of it as scaling with cloud LWP.

Figure 2 MAC-LWP-observed (a) and UM-simulated (b) LWP in the space of SST and WVP in the region 30–60^∘ N. White lines show the distribution of data in SST–WVP space. The contours of number of observations in each bin are shown in white.

Training the regression model in Eq. (4) produces a measure of the covariance between predictors and predictands. In particular, it produces a measure of the covariability between ln (N_d) and LWP (a₁). However, this covariance does not show causality. N_d may enhance LWP by, for instance, suppressing rain (Albrecht, 1989), and it is this relationship where causality runs from N_d to LWP that is of interest in understanding adjustments. However, the observed covariability between N_d and LWP is also the product of causality running from LWP to N_d due to precipitation removing aerosol and reducing N_d (Wood et al., 2012). The overall covariance between N_d and LWP is the product of these two causal pathways as well as correlation induced by environmental factors driving both clouds and aerosol and by geographic and seasonal distribution of sources of aerosol relative to cloud. Figure 3 shows a schematic of these causal flows.

Figure 3 A schematic illustration of the mechanisms linking cloud, aerosol, and meteorology considered in this study. Observed quantities are shown in green (in this case N_d and LWP). Adjustments (shown in red) are the causal link that flow from N_d to LWP, and inferring them is the focus of this study. In the “scavenging-only” simulation the N_d seen by the cloud physics is set to a constant value.

How do we constrain adjustments based on examination of modeled and observed N_d and LWP if we cannot disentangle causality? The solution we propose is to use a simplified version of the model to characterize the sensitivity of N_d to precipitation and other factors and use this sensitivity to interpret the full model. To conceptualize why this process can estimate the sensitivity of LWP to N_d we can write an equation describing small changes in LWP in the space of small changes in N_d. First, we consider that N_d acting on cloud microphysics ($N_{{\mathrm{d}}_{\text{Cloud}}}$) and N_d when precipitation scavenging is the only causal factor linking cloudiness and N_d ($N_{{\mathrm{d}}_{\text{Scav}}}$) as separate entities. In the latter case air mass history and sinks unrelated to precipitation can affect N_d, in addition to the effect of precipitation scavenging, but precipitation scavenging is the only causal link between cloud and aerosol (see Fig. 3). We consider changes in LWP relative to $N_{{\mathrm{d}}_{\text{Cloud}}}$ and $N_{{\mathrm{d}}_{\text{Scav}}}$. That is to say we consider $\text{LWP}=fn\left(\ln \right(N_{{\mathrm{d}}_{\text{Cloud}}}),\ln (N_{{\mathrm{d}}_{\text{Scav}}}\left)\right)$, where the logs are added for consistency with Eq. (4) and all other terms in Eq. (1) are held fixed and the first term is expanded.

Figure 4 A simple conceptual model demonstrating how variance due to precipitation driving of N_d and due to N_d driving changes in LWP through cloud microphysics affects total covariance between N_d and LWP in observations and models. N_d where scavenging is the only causal link between cloud and aerosol (N_dScav) and N_d affecting the cloud microphysics (N_dCloud) are considered to be distinct entities. The contour map in (a) illustrates a potential topological LWP space in terms of N_dCloud and N_dScav. The case where N_dCloud is set to a constant value is shown with a blue line, and the case where N_dCloud is increased independently is shown with a red line. In a causal sense this implies N_dCloud driving LWP (shown in b) and N_dScav being driven by LWP (shown in c), where N_dScav is negatively correlated with LWP because $\mathrm{d}{N}_{\text{dScav}}/\mathrm{d}R<\mathrm{0}$ and $\mathrm{d}\text{LWP}/\mathrm{d}R>\mathrm{0}$, where R is rain rate. The observed variability in LWP in terms of N_d is shown in (d) for the cases N_dCloud=N_d and N_dScav=constant (red), N_dScav=N_d and N_dCloud=constant (blue), and $N_{\text{dCloud}}=N_{\text{dScav}}=N_{\mathrm{d}}$ (the real world; dashed line).

Conceptually, LWP is increased with increasing $N_{{\mathrm{d}}_{\text{Cloud}}}$ (Fig. 4b), and $N_{{\mathrm{d}}_{\text{Scav}}}$ is reduced by increasing rain rates, which increase with LWP (Fig. 4c). In the topographic space describing LWP in terms of $N_{{\mathrm{d}}_{\text{Cloud}}}$ and $N_{{\mathrm{d}}_{\text{Scav}}}$, reduced $N_{{\mathrm{d}}_{\text{Scav}}}$ corresponds to increased precipitation rates and, in turn, LWP. That is to say,

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}{N}_{\text{dScav}}}{\mathrm{d}R}}<\mathrm{0},{\displaystyle \frac{\mathrm{d}R}{\mathrm{d}\text{LWP}}}>\mathrm{0},\end{array}$$

where R is the rain rate, so the reciprocal of the product of these terms gives

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}\text{LWP}}{\mathrm{d}{N}_{\text{dScav}}}}<\mathrm{0}.\end{array}$$

Thus, a small change in LWP in terms of changes in these two variables alone is

$$\begin{array}{ll}{\displaystyle }\mathrm{d}\text{LWP}=& {\displaystyle }{\displaystyle \frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Scav}}}\right)}}\mathrm{d}\ln \left(N_{{\mathrm{d}}_{\text{Scav}}}\right)\\ \text{(7)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Cloud}}}\right)}}\mathrm{d}\ln \left(N_{{\mathrm{d}}_{\text{Cloud}}}\right).\end{array}$$

We assume that scavenging is unaffected by changes in $N_{{\mathrm{d}}_{\text{Cloud}}}$ affecting rain rates, so there are no higher-order terms, although this may not be the case in reality (Wood et al., 2012), and this feedback is found to substantially enhance the strength of adjustments in some models (Jing and Suzuki, 2018). A visualization of Eq. (7) is shown in Fig. 4a. The best fit of LWP to N_d gives the sum of the partial derivatives of the two terms

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{d}\text{LWP}}{\mathrm{d}\ln \left(N_{\mathrm{d}}\right)}}=& {\displaystyle }{\left({\displaystyle \frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Scav}}}\right)}}\right)}_{N_{\text{dCloud}}}{\displaystyle \frac{\mathrm{d}\ln \left(N_{{\mathrm{d}}_{\text{Scav}}}\right)}{\mathrm{d}\ln \left(N_{\mathrm{d}}\right)}}\\ \text{(8)}& {\displaystyle }& {\displaystyle }+{\left({\displaystyle \frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Cloud}}}\right)}}\right)}_{N_{\text{dScav}}}{\displaystyle \frac{\mathrm{d}\ln \left(N_{{\mathrm{d}}_{\text{Cloud}}}\right)}{\mathrm{d}\ln \left(N_{\mathrm{d}}\right)}},\end{array}$$

where subscripts on the partial derivatives indicate that they are evaluated while holding that variable constant. In the control simulation, $N_{\mathrm{d}}=N_{{\mathrm{d}}_{\text{Cloud}}}=N_{{\mathrm{d}}_{\text{Scav}}}$ because N_dCloud and N_dScav are not distinct entities. That is to say, the cloud microphysics and radiation see that the N_d calculated by the nucleation scheme and aerosol experiences the effects of scavenging and meteorology. In this particular case the derivatives of $N_{{\mathrm{d}}_{\text{Scav}}}$ and $N_{{\mathrm{d}}_{\text{Cloud}}}$ each go to unity because the N_d seen by the cloud and affected by scavenging is equal, and we are left with

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\mathrm{d}\text{LWP}}{\mathrm{d}\ln \left(N_{\mathrm{d}}\right)}}={\left({\displaystyle \frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Scav}}}\right)}}\right)}_{N_{{\mathrm{d}}_{\text{Cloud}}}}+{\left({\displaystyle \frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Cloud}}}\right)}}\right)}_{N_{{\mathrm{d}}_{\text{Scav}}}}.\end{array}$$

We always sample along the dashed line in Fig. 4a in the observational record or in a model where clouds are sensitive to N_d and aerosol is scavenged, giving the curve in Fig. 4d. This curve combines the effects of causality flowing from N_d to LWP and from LWP to precipitation to aerosol and onto N_d. The term $\frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Cloud}}}\right)}$ is the key to constraining adjustments because it describes the sensitivity of LWP to changes in N_d (in this analogy the slope of LWP along the red line in Fig. 4a). While we cannot disentangle causality in observations or in a model where N_d affects cloud microphysics, we propose that in a model configuration where cloud microphysics (e.g., radiation, settling, and autoconversion in UM GA7.1) are unaffected by N_d we may estimate the term $\frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Scav}}}\right)}$. That is to say, if $N_{{\mathrm{d}}_{\text{Cloud}}}=\text{constant}$ and $N_{\mathrm{d}}=N_{{\mathrm{d}}_{\text{Scav}}}$, then

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{d\text{LWP}}{d\ln \left(N_{\mathrm{d}}\right)}}={\left({\displaystyle \frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Scav}}}\right)}}\right)}_{N_{{\mathrm{d}}_{\text{Cloud}}}},\end{array}$$

and the best fit of LWP to N_d in the scavenging-only simulation is just a measure of the scavenging of aerosol by precipitation and the non-causal covariance between N_d and LWP induced by other confounding factors.

Thus, the sensitivity $\frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Cloud}}}\right)}$ can be estimated as the difference between the regression of LWP on N_d trained in a simulation where clouds are insensitive to N_d ($N_{{\mathrm{d}}_{\text{Cloud}}}=\text{constant}$; scavenging only) and the regression in a control simulation ($N_{\mathrm{d}}=N_{{\mathrm{d}}_{\text{Cloud}}}=N_{{\mathrm{d}}_{\text{Scav}}}$) because

$$\begin{array}{}\text{(11)}& {\left({\displaystyle \frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Cloud}}}\right)}}\right)}_{N_{{\mathrm{d}}_{\text{Scav}}}}={\displaystyle \frac{\mathrm{d}\text{LWP}}{\mathrm{d}\ln \left(N_{\mathrm{d}}\right)}}-{\left({\displaystyle \frac{\partial \text{LWP}}{\partial \ln \left(N_{{\mathrm{d}}_{\text{Scav}}}\right)}}\right)}_{N_{{\mathrm{d}}_{\text{Cloud}}}}.\end{array}$$

We will use this correction for variability induced by factors besides adjustments throughout the paper to infer the effects of N_d on cloud properties in model configurations and the observations. We note that this methodology is subject to several limitations. First, as noted above, it assumes that feedbacks between adjustments and scavenging are small. This assumption will break down as adjustments become stronger. However, the response in cloud to transient changes in aerosol does not appear to permit very strong adjustments (Malavelle et al., 2017; Toll et al., 2017, 2019). Second, it is assumed here that the relationship between LWP and ln (N_d) is linear. Empirically this assumption appears to hold in the models and observations (as shown in the following section), and inclusion of higher-order polynomials of ln (N_d) in Eq. (3) did not alter our results.

To reiterate, non-causal factors also lead to correlation between N_d and LWP because air mass history is important for both clouds and aerosol (Mauger and Norris, 2007). In the discussion of the conceptual model presented above we focused on the impact of changes in N_d on LWP and vice versa. Spurious correlation due to air mass history will affect both terms on the right-hand side of Eq. (11) when they are calculated using regression. Thus, the correction for scavenging-induced variability in Eq. (11) also corrects for variability induced by air mass history simultaneously affecting N_d and LWP.

3 Results

3.1 Variance in LWP as a function of N_d

As described in the Methods section, data from observations and the control simulation in the UM between 30 and 60^∘ N over oceans and 2000 km from cyclone centers at a daily mean resolution are binned by WVP and SST, and the regression model shown in Eq. (4) is trained in each bin. If we plot the quantity,

$$\begin{array}{}\text{(12)}& \text{LWP}-a_{\mathrm{2}}{\mathit{\omega }}_{\mathrm{500}}-a_{\mathrm{3}}\text{EIS}-a_{\mathrm{4}}\text{WVP}-a_{\mathrm{5}}\text{SST}-a_{\mathrm{6}}\text{SHF}-a_{\mathrm{7}},\end{array}$$

as a function of N_d we can see that the coefficient relating LWP to N_d is greater than zero at 99 % confidence, once variability associated with meteorological predictors is removed (Fig. 5). The residual variance in LWP covaries with N_d similarly between the observations and UM simulations where observations exist. However, one key difference is the span of N_d values in the observations and simulations. Observed N_d almost never falls below 30 cm⁻³ (the grey line in Fig. 5), while the N_d calculated by the model occasionally falls below this value (∼ 10 % of modeled ln (N_d)). This is not surprising because the model suffers from no retrieval limitations and can always measure an N_d, even when cloud is extremely tenuous. Ultimately the goal of this study is to utilize the slope of LWP with respect to N_d to infer adjustment strength. However, using the slope from the observations and the models will produce quite different values due to the strong increase in LWP with N_d below 30 cm⁻³. In order to compare observations and models fairly we exclude data for $N_{\mathrm{d}}<\mathrm{30}{\mathrm{cm}}^{-\mathrm{3}}$.

Figure 5 Residual variance in LWP after removing variability associated with other predictors (Eq. 12) plotted as a function of N_d in the UM GA7.1 and observations (Obs.). Each subplot corresponds to a bin of WVP and SST. Due to retrieval limitations, values of N_d below 30 cm⁻³ (grey line) are almost never observed but occur in model output.

For $N_{\mathrm{d}}>\mathrm{30}{\mathrm{cm}}^{-\mathrm{3}}$ residual LWP decreases with increasing N_d in both observations and models. This agrees with the notion that precipitation variability drives N_d (Wood et al., 2012) and agrees with previous analysis (Gryspeerdt et al., 2019). Clearly, this does not mean that the model reduces cloudiness in response to aerosol. The cloud microphysics within the UM GA7.1 must increase LWP in response to a change in N_d, and analysis of the model response to changes from preindustrial (PI) to present-day (PD) aerosol confirmed a negative forcing from adjustments (Mulcahy et al., 2018). As discussed in the Methods section, the difference in slopes between the scavenging-only and control simulation is a proxy for the response in cloud microphysics to N_d. Before analyzing this difference, we will briefly discuss the covariance between LWP and meteorological predictors in the observations and the control simulations.

3.2 Covariances between LWP and meteorology

Before returning to discussing aerosol–cloud adjustments, we discuss the covariances between observed LWP and the other predictors in Eq. (4). This is shown in Fig. 6. The explained variance (R²) by the regression model exceeds 40 % except at low WVP and SST. The variance explained across all bins of SST and WVP is 46 % of daily mean 1^∘ × 1^∘ variability. That is to say, nearly half of the day-to-day variability in LWP across the midlatitudes away from cyclones can be explained as a simple linear combination of five variables.

Figure 6 The regression coefficients relating observed LWP to predictors in the phase space of SST and WVP. Regression coefficients not significant at 99 % confidence are marked with a dot. White contours show the distribution of observations in the WVP–SST phase space. The first plot shows the explained variance in each bin. The explained variance within the entire data set is noted in the title along with the number of observations. The remaining plots show the coefficients relating each predictor to LWP. In the title the weighted mean of the absolute value of the standardized coefficient ($\left|\partial \mathit{\sigma }\text{LWP}/\partial \mathit{\sigma }x\right|)$ is shown in brackets to give an estimate of contribution of each predictor to the variance. Regression coefficients are as described in Eq. (4).

Several of the predictors have uniform effects on LWP: LWP decreases with subsidence (consistent with Myers and Norris, 2013), increases with WVP, decreases with increasing SST (consistent with Qu et al., 2015, and references therein), and increases with sensible heat flux (consistent with Miyamoto et al., 2018). We note that at very low SST it appears that LWP may increase with SST. It is possible that this feature is related to ice-to-liquid transitions (McCoy et al., 2015b, 2016; Tsushima et al., 2006; Senior and Mitchell, 1993; Tan et al., 2016), but this region only accounts for a small fraction of the overall data volume. The only meteorological predictor that correlates positively and negatively with LWP across the SST–WVP phase space is EIS. At higher WVP and lower SST, EIS and LWP are negatively correlated, while at low WVP and high SST, they are positively correlated. The latter effect is consistent with increasing EIS increasing the cloud fraction (Wood and Bretherton, 2006). One possibility is that because clouds are closer to 100 % in areal coverage at higher WVP and lower SST, increasing stability suppresses the thickening of cloud in this regime, while the clouds cannot expand horizontally.

Examination of the control simulation shows very similar covariances between meteorological predictors and LWP (Fig. 7). Notable differences are the stronger positive covariance between SST and LWP at low SST and uniform negative covariance between EIS and LWP across the SST–WVP phase space. However, the relationship between LWP and meteorology is strikingly similar to that between observations and the UM. The explained variance by predictors tends to be higher in the UM GA7.1 (64 %), but it is hard to say how much of this difference in explained variance is due to simplification of the real atmosphere and cloud physics by the model and how much is due to observational error.

Figure 7 As in Fig. 6 but showing the covariance within the UM.

In this framework we also characterize the covariance between ln (N_d) and LWP. The correlation between ln (N_d) and LWP is primarily negative in the observations, with positive correlations only occurring at high WVP and low SST. The correlation between LWP and ln (N_d) in the model is almost uniformly negative (in the data set restricted to $N_{\mathrm{d}}>\mathrm{30}{\mathrm{cm}}^{-\mathrm{3}}$). Standardization of the predictor strength estimates the effect of a standard deviation change in each predictor in standard deviations of LWP. In both the observations and models the impact of a standard deviation in meteorological predictors dwarfs the effect of a standard deviation in ln (N_d), with contributions from subsidence and SHF dominating variability. That is to say, the relative contributions of variability in N_d are quite small in comparison to variability that is simply due to the weather.

3.3 Inferring the effect of N_d on LWP

As we have seen, the correlation between LWP and N_d is mostly negative. The slope relating ln (N_d) to LWP in the context of Eq. (4) is reproduced in Fig. 8a and b. However, if we manually set the N_d seen by autoconversion, radiation, and settling to 75 cm⁻³ in the UM GA7.1 (the approximate mean value in the study region) the negative covariance strengthens substantially (Fig. 8c). This model version is termed “scavenging only” because scavenging is the only causal link between cloud and aerosol (Fig. 3). As discussed in the Methods section, the effect of changes in N_d on LWP may be approximated by the difference between Fig. 8b and c. This is shown in Fig. 8e. This agrees with our expectation based on our knowledge of how the UM GA7.1 treats liquid cloud processes. We know that the effect of increasing N_d on autoconversion and settling is to inhibit the conversion of cloud to rain and reduce cloud droplet fall speed.

Figure 8 The observed covariance between N_d and LWP and the inferred effect of N_d on LWP. Notations in the bottom right of plots indicate the expected direction of causality. The covariance between N_d and LWP is shown in (a) the observational record (as in Fig. 6) and (b) the control simulation of the UM GA7.1 (as in Fig. 7); (c) a simulation where the N_d seen by the autoconversion, settling, and radiation is set to 75 cm⁻³ and only scavenging links N_d and LWP. The inferred strength of adjustments in observations (d) calculated as the difference in slope between (a) and (c) as in Eq. (11), and in (e) the UM is calculated as the difference in slope between (b) and (c). Inferred and true adjustment strength for the UM is compared in (f). The inferred adjustment strength in (e) is compared to the sensitivity of LWP to N_d in a set of simulations where N_d is fixed in the autoconversion, settling, and radiation to a range of values, forcing causality to flow from N_d to LWP. The 1 : 1 line is shown with dashes. Different symbols indicate whether the N_d seen by the autoconversion, settling, or radiation was varied while others were held constant. Each symbol corresponds to a bin in WVP and SST, as in (e).

To evaluate the accuracy of this method in inferring the sensitivity of LWP to N_d we conduct simulations where causality is forced to flow from the N_d to LWP. This is done by setting the N_d seen by the autoconversion, settling, and radiation to constant values. Two simulations are conducted for each mechanism, with N_d set to 100 and 300 cm⁻³. The N_d seen by the other mechanisms is held constant at 30 cm⁻³. A further three simulations are conducted, in which N_d seen by all three mechanisms is set to 30, 100, and 300 cm⁻³ (Table 2). This yields a set of simulations describing the effects on LWP of an increase in N_d in the autoconversion, settling, or radiation as well as via all three mechanisms simultaneously. Analysis of the effects of increased N_d seen by the radiation code showed negligible changes in LWP and is not shown here.

As in Fig. 7, the fit in Eq. (4) is trained in each bin of WVP and SST across each set of three simulations. The N_d seen by the model parameterization is used in the multiple linear regression (Eq. 4). As expected, given the functional form of these parameterizations (Eqs. 1 and 2) the effect of increasing N_d in either autoconversion or settling results in an increase in LWP, as does the case where N_d in the settling, autoconversion, and radiation are all varied together. The effect of a perturbation in the N_d seen by autoconversion, radiation, and settling on LWP is nearly identical to the slope inferred by subtracting the slope in the scavenging-only run from the control run following Eq. (11) (Fig. 8f). This comparison between the true model sensitivity of LWP to changes in N_d supports the validity of the multiple-linear-regression-based inference approach pursued here.

Autoconversion and settling are likely to interact with each other in a non-linear manner. As one is suppressed as a sink of cloud liquid, the other will grow (for example, if the conversion of cloud to rain is rendered inefficient, then the gravitational settling of cloud droplets will become more efficient as the liquid content grows). Analysis of crosstalk between these terms is not the focus of this article, but sensitivity of LWP to changes in each mechanism is shown in Fig. 8f. The effect of perturbing settling N_d on LWP is larger than the effect of perturbing the autoconversion N_d, singling it out as an important lever in controlling adjustment strength for future study. This result is consistent with the large process rate for settling relative to autoconversion found in weakly precipitating low-cloud regimes in CAM5 (Gettelman et al., 2013). It is also possible that interactions between settling and N_d may decrease the sensitivity of warm rain to autoconversion representation that is diagnosed in other models (Michibata et al., 2019; Jing et al., 2019). Overall, the results of this study are insensitive to whether settling or autoconversion has a stronger effect on aerosol–cloud adjustments because we assume that feedbacks between scavenging and adjustments are small.

Given that the scavenging-only simulation can be used to infer the adjustment strength in the UM, we also attempt to use the scavenging-only simulation to infer the strength of adjustments in the real world (Fig. 8d). This assumes that the statistical relationship between LWP and rain rates in the UM GA7.1 and the scavenging of CCN-relevant aerosol by precipitation in the UM GA7.1 is realistic.

Based on the covariance between N_d and LWP in the scavenging-only simulations in the UM GA7.1 and observations we offer an estimate of adjustment strength. The adjustment strength inferred from observations is stronger than the UM for low SST and high WVP. Unlike the UM, weak negative covariance between LWP and N_d exists for high SST and low WVP. This may indicate that the scavenging of aerosol by the UM is too efficient, or some other confounder of the relationship between LWP and N_d is poorly represented in the model, or it may point to the thinning of cloud via enhanced entrainment (Ackerman et al., 2004), which is not represented in the UM. Overall, the effect of increasing N_d in the UM is an increase in LWP.

We have discussed how to infer adjustment strength from the total covariance between N_d and LWP, which mixes covariance induced by causality flowing from cloud to aerosol and aerosol to cloud. In the next section we show how this inferred adjustment strength from the simple regression model is able to predict LWP adjustments in response to anthropogenic aerosols within the full GCM simulation.

3.4 Reproducing simulated adjustments between the preindustrial and present-day aerosol emission

We have produced an estimate of the sensitivity of LWP to changes in N_d in the UM GA7.1 model (Fig. 8e). Does this simple regression model hold skill in reproducing the actual change in LWP within the model between the present-day (PD) and preindustrial (PI) aerosol emissions? Simulations are carried out by setting aerosol emission to preindustrial values. This is done for simulations using the GA7.1 scavenging parameterization. A second set of simulations is conducted where the ECHAM-HAM5 scavenging coefficients are used to estimate the sensitivity to the representation of scavenging. PI, PD, and scavenging-only simulations are conducted for each scavenging configuration.

The change in LWP between the PI and PD (ΔLWP_PD−PI) simulated by the UM is 2.9 g m⁻², averaged across the 30–60^∘ N region. Using the UM-simulated change in N_d (ΔN_dPD−PI) and the slope relating ln (N_d) to LWP diagnosed from present-day variability (Fig. 8) we calculate ΔLWP_PD−PI. If the covariance between ln (N_d) and LWP observed in the control run (Fig. 8b) of the model and ΔN_dPD−PI is used to predict ΔLWP_PD−PI, the negative correlation introduced by precipitation scavenging leads to a negative ΔLWP_PD−PI (Fig. 9a; in disagreement with the UM-simulated ΔLWP_PD−PI). If the covariance in the scavenging-only simulation (causality flows from clouds to precipitation to aerosol to N_d) is used (Fig. 8c), the regression-model-predicted decrease in LWP between PI and PD aerosol levels doubles. However, if we subtract the sensitivity in scavenging-only simulation from the control simulation to yield the sensitivity of LWP to changes in N_d (Fig. 8e), the ΔLWP_PD−PI predicted by the regression model and the UM-simulated ΔN_dPD−PI agree well with the UM-simulated ΔLWP_PD−PI (Fig. 9). This result supports the utility of the inferred sensitivity of LWP to N_d in examining adjustments due to PD aerosol.

Figure 9 The change in LWP between the simulations with PI and PD aerosol emissions (ΔLWP_PD−PI) inferred by covariance between N_d and LWP in the PD and as simulated by the UM. The inferred ΔLWP_PD−PI is calculated using the UM-simulated change in N_d. Panel (a) shows the inferred ΔLWP_PD−PI based on variance in different simulations. The direction of causality is indicated for each case. From left to right: ΔLWP_PD−PI inferred from covariance in the control simulation (N_d⇔LWP, causality between N_d and LWP goes both directions; Fig. 8b), inferred from the scavenging-only simulation (N_d⇐LWP, causality from LWP to N_d; Fig. 8c), inferred by correcting the total covariance using covariance in the scavenging-only simulation (N_d⇒LWP, causality goes from N_d to LWP Fig. 8d – this is the adjustment strength in UM GA7.1 inferred by the method presented in this paper) and when the covariance in observations is combined with the scavenging-only simulation (listed as Obs. constraint). The ΔLWP_PD−PI simulated by the GCM is shown on the right in black. This is the true aerosol–cloud adjustment in the GCM that is compared to the value being inferred in the UM GA7.1, as shown by the red bars. ΔLWP_PD−PI is provided for simulations using the UM GA7.1 scavenging coefficients and ECHAM-HAM5 scavenging coefficients. Panel (b) as in (a) but resolved in latitude and only showing the case when the UM GA7.1 default scavenging is used.

What does the covariance between LWP and ln (N_d) in the real world tell us about the model? Use of the observational estimate of adjustment strength (Fig. 8d) calculates a ΔLWP_PD−PI that is two-thirds of the ΔLWP_PD−PI inferred from the control run and ΔN_dPD−PI (Fig. 9). This is because, while the inferred sensitivity of LWP to N_d from observations is quite strong in some regions, it is weakly negative in the most commonly occurring WVP–SST regimes (Fig. 8d). That is to say, if the efficacy of precipitation in removing CCN-relevant aerosol in the real atmosphere is near to the efficacy in the UM, the relationship between LWP and rain rates in the real atmosphere is near to the relationship in the UM, and the covariability between N_d and LWP induced by other confounders of the N_d–LWP relationship is realistic, then the adjustments (ΔLWP_PD−PI) simulated by the UM are not unreasonable, albeit a little large. This result is consistent with constraints provided by simulation of the Holuhraun eruption (Malavelle et al., 2017), which showed that the anomaly predicted by HadGEM3 was within the observed range. However, the observationally inferred sensitivity of LWP to changes in N_d is stronger in low-SST and high-WVP regimes than it is in the UM GA7.1 but much weaker, or even slightly negative, in low-WVP and high-SST regimes. This results in a strong latitudinal gradient in ΔLWP_PD−PI between 30 and 60^∘ N inferred from observations (Fig. 9b).

The predicted ΔLWP_PD−PI based on the control and scavenging-only simulations (Fig. 8e) underestimates the UM-simulated ΔLWP_PD−PI by around 10 %. This may result from disregarding feedbacks between precipitation suppression and N_d (Jing and Suzuki, 2018), may simply be due to shortcomings in the simple linear model relating ln (N_d) to LWP, or may relate to meteorological differences between the PI and PD simulations not accounted for by meteorological predictors or compositing. However, the fractional reduction in inferred ΔLWP_PD−PI when observations are used to constrain the sensitivity of LWP to N_d suggests that ΔLWP_PD−PI should be around 60 %–70 % of the value predicted by the UM GA7.1 averaged across the NH. ΔLWP_PD−PI inferred from observations is extremely dependent on latitude, with almost no change in LWP over lower latitudes and warmer SSTs. We hasten to note that the change in N_d between the PI and PD is uncertain and its prediction by UM GA7.1 may be biased, making the absolute values of the constrained ΔLWP_PD−PI less relevant than the fractional overestimation in UM-simulated ΔLWP_PD−PI. This analysis is repeated using simulations that are only nudged above level 27 of the model (Fig. S1) and not nudged at all (Fig. S2), with minimal differences from the results presented in the main text.

We have focused on changes in LWP in this work because it is a variable that we have good observations of, it can be compared between models and observations in a straightforward way, and it clearly links to adjustments. However, this variable is not the key variable in discussing radiative forcing and climate sensitivity. How does the radiative forcing from adjustments scale with LWP? Determining the precise effects of adjustments on the shortwave cloud radiative effect is difficult, but if we assume that the perturbations in LWP induced by adjustments are similar to the perturbations that are driven by meteorology, we can offer a simple estimate to inform our understanding of the modeled forcing from adjustments.

Examination of the relationship between LWP and albedo within the UM GA7.1 shows a rapid climb in albedo for low LWP, followed by saturation at higher values, as expected from saturation of the cloud fraction and optical depth (Fig. 10). This curve can be fit by a second-degree polynomial, meaning that the sensitivity of albedo to changes in LWP is a function of the LWP. That is to say, regions where LWP is already high are going to have a smaller increase in albedo for a unit change in LWP. The change in LWP between PI and PD is much smaller than the range of the variation in LWP in the climate mean state, so we use the average of monthly mean LWP between PI and PD to calculate the sensitivity of albedo to changes in LWP for the NH midlatitudes. This sensitivity (dα∕dLWP) is multiplied by the change in LWP to yield change in albedo. This change in albedo is scaled by the downwelling shortwave to give the change in reflected shortwave.

Figure 10 UM GA7.1 daily mean albedo as a function of LWP over ocean between 30 and 60^∘ N and outside of cyclones in equal quantiles of LWP (black circles). A second-order polynomial fit is shown using a red line.

The change in reflected shortwave predicted from LWP changes in the UM GA7.1 is 1.9–2.0 W m⁻², depending on the scavenging parametrization used. The change in shortwave inferred from the regression model of LWP trained in the control run and corrected by the scavenging-only run is 1.5 W m⁻². If the observed sensitivity of LWP to N_d is used to constrain ΔLWP, the predicted change in reflected shortwave is approximately 1.0 W m⁻². Thus, we estimate that GA7.1 overpredicts the change in reflected shortwave due to adjustments in response to a given change in N_d by around 50 %. This estimate is subject to the caveat that changes in LWP due to adjustments may not affect albedo in the same way as suggested by examining the total variability. For example, adjustments might only increase the liquid content of the thickest clouds and have a relatively slight impact on albedo. However, it is unclear how to provide a more complex calculation than the one presented here.

4 Discussion

As the possible range for the radiative forcing from the first indirect effect has narrowed, aerosol–cloud adjustments have become an increasingly central source of uncertainty in aerosol–cloud radiative forcing (Bellouin et al., 2020). Here we focus on the northern midlatitudes, where the majority of anthropogenic aerosol is emitted (Myhre et al., 2013). Previous work has examined aerosol–cloud adjustments in midlatitude cyclones, showing that cyclone liquid water path (LWP) increases with cloud droplet number concentration (N_d; McCoy et al., 2018b). This work focuses on the remainder of cases in the midlatitudes, when there is no cyclone center within 2000 km (roughly 42 % of the time between 30 and 60^∘ N; see Fig. 1).

Untangling the effect of cloud microphysics on cloud macrophysics from the total variability in cloud macrophysics is challenging. In the observational record we can only characterize the covariance between predictors and predictands. In interrogating the observed covariances between cloud properties and different meteorological predictors we find many of the relationships that have been described in the literature documenting cloud-controlling factors (Myers and Norris, 2015). Once meteorological variability is accounted for, statistically significant relationships between cloud microphysics (N_d) and cloud macrophysics (LWP) appear. In keeping with previous studies (Gryspeerdt et al., 2019), LWP and N_d are found to be negatively correlated. However, this is clearly not consistent with the UM GA7.1 GCM's actual response in LWP to changes in aerosol emissions between PI and PD. This negative correlation is due to a combination of N_d being driven by scavenging, as shown in Wood et al. (2012), and spurious correlation between N_d and LWP driven by external variables affecting both terms (Mauger and Norris, 2007), overwhelming any positive covariance driven by aerosol–cloud adjustments.

While we cannot disentangle the covariance due to adjustments and covariance due to scavenging and other confounders in the observational record, we can create simulations in which causality is forced to flow from clouds to N_d (scavenging only). Using this measure of non-adjustment-induced variability in N_d, we can infer the effect of changing N_d on LWP. This inferred sensitivity of LWP to N_d agrees well with simulations where N_d is manually varied in the cloud microphysics (Fig. 8f). The inferred sensitivity of LWP to N_d, combined with the UM-predicted change in N_d, reproduces the UM-predicted change in LWP between the PI and PD (ΔLWP_PD−PI; Fig. 9). These two tests support this method's relevance to understanding aerosol–cloud adjustments.

The analysis presented here, combined with the enhancement in LWP in cyclone systems shown in McCoy et al. (2018b), points toward an overall increase in LWP across the NH midlatitudes in response to anthropogenic aerosol. Ultimately, while the regimes examined here and in McCoy et al. (2018b) are very different, the detection of a change in LWP in response to changes in N_d rests on the ability of the technique to account for non-adjustment-induced variability in N_d. In McCoy et al. (2018b) this was done by stratifying the data set by the cyclone precipitation rate, which is well-constrained by the large-scale environment (Field and Wood, 2007); this in turn stratifies the data set by the scavenging of aerosol.

Assuming that scavenging in the UM GA7.1 model is realistic, that the relationship between LWP and precipitation is reasonable, and that the non-causal covariance produced by other factors is replicated in the model, we evaluate the strength of adjustments based on observations. This reveals that the present version of the UM GA7.1 overestimates the sensitivity of LWP to changes in N_d by approximately 50 % outside of midlatitude cyclones. Calculation of the implied change in shortwave shows a similar overestimate in forcing due to adjustments by the UM. Observations also imply that aerosol–cloud adjustments in the UM GA7.1 occur in the wrong regime. Adjustments in the UM favor warmer SSTs and lower WVP, while the observations favor colder SSTs and higher WVP. This difference in regime may be due to early precipitation onset in the UM GA7.1.

The present study and McCoy et al. (2018b) project a negative forcing due to adjustments in the midlatitudes. Extension of this study to the subtropics and tropics is reserved for a future analysis. While extrapolating the present results to the tropics is difficult, it seems reasonable to suppose that the pattern of adjustments in the subtropics is relatively similar to the inter-cyclone regime analyzed here and to expect a negative forcing in this region due to adjustments, although this may be weaker than the midlatitude adjustments due to decreasing LWP for increasing N_d over higher SSTs and lower WVP over some of the subtropics (Fig. 8d). One important limitation of this study is that it assumes weak feedbacks between adjustments and scavenging. If some regimes have very strong feedback between precipitation and adjustments, then our analysis will infer adjustments that are too weak.

Our result is in contradiction with previous empirical constraint studies that have postulated that changes in N_d greatly enhance LWP (Rosenfeld et al., 2019), have little effect (Malavelle et al., 2017; Toll et al., 2017), or reduce LWP (Sato et al., 2018; Gryspeerdt et al., 2019). We suggest that this diversity within the literature is because a range of constraints may be arrived at, depending on the degree to which precipitation scavenging and meteorological driving of aerosol and cloud occludes aerosol–cloud adjustments and the steps that are taken in the analysis to account for scavenging-induced and meteorologically induced covariability. Based on the analysis presented here we believe that positive, zero, or extremely strongly negative radiative forcings due to aerosol–cloud adjustments in the midlatitudes are not supported by the observations.