Constraining the aerosol inﬂuence on cloud liquid water path

. The impact of aerosols on cloud properties is one of the largest uncertainties in the anthropogenic radiative forcing of the climate. Signiﬁcant progress has been made in constraining this forcing using observations, but uncertainty remains, particularly in the magnitude of cloud rapid adjustments to aerosol perturbations. Cloud liquid water path (LWP) is the leading control on liquid-cloud albedo, making it important to observationally constrain the aerosol impact on LWP. Previous modelling and observational studies have shown that multiple processes play a role in determining the LWP re- 5 sponse to aerosol perturbations, but that the aerosol effect can be difﬁcult to isolate. Following previous studies using mediating variables, this work investigates use of the relationship between cloud droplet number concentration (N d ) and LWP for constraining the role of aerosols. Using joint probability histograms to account for the non-linear relationship, this work ﬁnds a relationship that is broadly consistent with previous studies. There is signiﬁcant geographical variation in the relationship, partly due to role of meteorological factors (particularly relative humidity). However, the N d -LWP relationship is negative in 10 the majority of regions, suggesting that aerosol induced LWP reductions could offset a signiﬁcant fraction of the instantaneous radiative forcing from aerosol-cloud interactions (RFaci). the N d -LWP the causal N d impact on LWP due to the role of confounding factors. The weaker LWP reduction implied by these “natural experiments” means that this work provides an upper bound to the radiative forcing from 15 aerosol-induced changes in the LWP.

radiative forcing due to changes in cloud properties vary significantly between different global climate models (Zelinka et al., 2014;Heyn et al., 2017), highlighting the need for observational constraints on the impact of aerosol on cloud properties.
Unlike greenhouse gases, aerosol properties vary strongly in space and time. This co-variation of aerosol and cloud properties in the present day atmosphere has been used to infer the impact of aerosols on cloud properties (e.g. Sekiguchi et al., 2003;Kaufman et al., 2005;Koren et al., 2005). Such observed relationships have been used to estimate the instantaneous radiative 5 forcing (RFaci) from a change in N d (e.g. Quaas et al., 2008;Stevens et al., 2017;McCoy et al., 2017;Gryspeerdt et al., 2017) and of the aerosol induced change in CF (Chen et al., 2014;Goren and Rosenfeld, 2014;Gryspeerdt et al., 2016;Christensen et al., 2017). As the leading order term for determining cloud albedo (Engström et al., 2015), it is also vital to constrain aerosol effects on the in-cloud liquid water path (LWP), separate from changes in the CF. Existing studies show a mixed picture; while some model (Quaas et al., 2009;Koren et al., 2014;Seifert et al., 2015;Grosvenor et al., 2017;Neubauer et al., 2017) and 10 observational studies (Gryspeerdt et al., 2014b;McCoy et al., 2018) suggest an increase in LWP with increasing aerosol, other studies (Wang et al., 2003;Small et al., 2009;Chen et al., 2014;Michibata et al., 2016;Christensen et al., 2017;Sato et al., 2018) find a reduction in LWP as aerosol increases. Some studies find both an increase and a decrease in LWP, depending on the meteorological conditions (Han et al., 2002;Ackerman et al., 2004;Bretherton et al., 2007;Xue et al., 2008;Toll et al., 2017;Bender et al., 2018), while other studies suggest a very weak LWP response to aerosol (Wang et al., 2012;Malavelle 15 et al., 2017). The main aim of this work is to reconcile these previous studies and develop a constraint on the aerosol impact on LWP.

Isolating an aerosol effect
The key difficulty in interpreting observed aerosol-cloud relationships is separating the causal impact of aerosols (the change in LWP caused by an aerosol perturbation) from the confounding role of local meteorology (e.g. Quaas et al., 2010) and retrieval 20 errors (e.g. Várnai and Marshak, 2009). Relative humidity in particular has been shown to obscure the causal relationship between aerosol optical depth (AOD) and CF (Quaas et al., 2010;Chand et al., 2012;Grandey et al., 2013). As many cloud properties are correlated to CF, the factors that obscure the aerosol-CF relationship can also confound other aerosol-cloud relationships, even those involving "intrinsic" cloud properties (Chen et al., 2014), such as cloud top pressure (Gryspeerdt et al., 2014a), and LWP Neubauer et al., 2017). Recent work (Gryspeerdt et al., 2016) has suggested 25 that the use of a mediating variable such as N d can be used to account for the confounding influence of relative humidity.
Following from this, the potential of the N d -LWP relationship to constrain the aerosol impact on LWP is investigated in this work.
Similar to the aerosol-LWP relationship, where both potential aerosol effects and confounders can influence the strength of the relationship, several effects may influence the observed N d -LWP relationship. 30 E1 Aerosol effects An increased aerosol concentration is likely to increase N d . This increase in N d may affect cloud processes and in turn modify the LWP. There are several hypothesised pathways for a causal effect of aerosol on LWP, varying in relative strength with the local meteorological conditions and aerosol environment: (a) Precipitation suppression (Albrecht, 1989) -an increased N d at initially unchanged LWP implies reduced cloud droplet sizes, suppressing the formation of precipitation. This reduction in the cloud water loss to precipitation could subsequently increase cloud depth (Pincus and Baker, 1994) and thus LWP. While it has been demonstrated that a reduction in droplet size suppresses precipitation (Suzuki et al., 2013), it is not clear how strongly this impacts LWP.

5
(b) The sedimentation-entrainment feedback (Ackerman et al., 2004;Bretherton et al., 2007) -the reduction in droplet radius from increased N d reduces the sedimentation flux in stratiform clouds, concentrating liquid water in the entrainment zone at the cloud top and increasing cloud-top evaporative and radiative cooling, increasing the entrainment rate. This increases the evaporative cooling in a positive feedback that depends on the above-cloud relative humidity, with drier air above cloud tops implying a larger LWP decrease. Negative N d -LWP relationships 10 in recent observational studies were suggested to have been due to this effect (Chen et al., 2014;Michibata et al., 2016;Sato et al., 2018).
(c) Evaporation-entrainment feedbacks (Wang et al., 2003;Xue and Feingold, 2006;Jiang et al., 2006;Small et al., 2009;Dagan et al., 2017) -smaller droplets have a faster evaporation timescale, enhancing the cooling and hence the negative buoyancy at the edge of cumulus clouds. This intensifies the horizontal buoyancy gradient, increasing 15 entrainment and hence evaporation, reducing the LWP with an expected similar meteorological dependency to E1b.
Aircraft observations have found increased horizontal buoyancy gradients and reductions in cloud liquid water content (LWC) in polluted clouds (Small et al., 2009).
(d) Warm cloud invigoration (Koren et al., 2014) -when N d is low, a lack of droplet surface area slows the cloud liquid water content (LWC) growth, increasing the local supersaturation. In this N d limited state, increasing the N d in 20 polluted clouds increases the LWC and so the latent heat release, allowing the cloud to achieve a larger vertical extent, which may increase the LWP.
E2 Retrieval errors The MODIS LWP and N d both depend on the retrieved cloud top droplet effective radius (r e ) and cloud optical depth (τ c ) and involve assumptions of varying validity (e.g., Grosvenor et al., 2018).
(a) Random errors in the retrieval of cloud properties (τ c , r e ) becoming correlated errors in LWP and N d . Using N d and 25 LWP calculated using the adiabatic assumption, random errors in τ c will generate a positive N d -LWP sensitivity , while errors in r e will generate a negative sensitivity ( d ln L d ln N d = −0.4), see appendix A for details.
(b) Sub-adiabatic clouds. Both the LWP and the N d retrieval make assumptions about the adiabaticity of clouds. Variations in the adiabaticity (Merk et al., 2015), even across a single cloud can therefore generate a positive N d -LWP

30
(c) Other systematic retrieval errors. Systematic biases in r e and τ c (particularly in broken cloud regions) may also affect the N d -LWP relationship. Other possibilities include variations in the vertical distribution of cloud water, assumptions about the droplet size spectrum, a dependence on satellite and solar zenith angle (Eastman and Wood, 2016;Grosvenor and Wood, 2014) and non-linearities in the retrieval (Zhang and Platnick, 2011 covariation problematic (Pearl, 1994).
(a) Precipitation preferentially occurs at large LWP. Precipitation scavenging of aerosol can reduce the amount of aerosol available for future activation to cloud droplets, reducing N d . Conversely, if an increased N d decreases the precipitation rate, this could result in a further increase in the N d through a reduction in wet scavenging and an increase in the available aerosol (a positive feedback).

10
(b) The impact of entrainment on the retrieved N d . The retrieved N d depends on the r e and the impact of entrainment on r e depends on the mixing type. Extreme inhomogeneous mixing (Baker et al., 1980) leads to a reduction in N d and LWP, but no immediate change in the droplet size distribution and hence no change in the r e or the retrieved N d . In contrast, homogeneous mixing (Warner, 1973)  These effects are depicted in Fig. 1. To constrain the causal aerosol influence on LWP, the impact of E1 has to be identified and isolated from that of E2-4. This would allow the aerosol impact on LWP to be constrained using the N d -LWP relationship. 25 It is necessary to understand the role of these different processes on the N d -LWP relationship in order to determine the impact of aerosols on the LWP. Using a variety of different satellite retrievals along with reanalysis data, the N d -LWP relationship is investigated globally and the impact of meteorology is explored. To understand the role of feedbacks (E3) and additional confounders (E4), natural experiments are used to examine the N d -LWP relationship in regions where there is a strong aerosol perturbation. Finally, the observed relationship is converted to a radiative forcing, allowing it to be compared to 30 other observational studies and to be used for further analysis of the aerosol impact on clouds and the climate.

Methods
This work is based on observational data from the Aqua satellite, specifically the moderate resolution imaging spectroradiometer (MODIS), the advanced microwave scanning radiometer-EOS (AMSR-E) and the clouds and the Earth's radiant energy system (CERES) instruments for a three year period (2007-2009 inclusive).
N d is retrieved using the level 2, collection 6, MODIS cloud property dataset (MYD06_L2) at a 1 km by 1 km resolution, 5 making use of the adiabatic assumption Quaas et al., 2006). Following the work of Grosvenor and Wood (2014) and Bennartz and Rausch (2017), the N d is filtered to include only liquid, single layer clouds with a top warmer than 268 K at 1 km resolution. In addition, pixels with an optical depth smaller than 4 or an effective radius less than 4 µm are excluded due to the uncertainty of these retrievals (Sourdeval et al., 2016). Pixels with a 5 km cloud fraction less than 0.9 are excluded to remove pixels close to cloud edges, and only pixels with a solar zenith angle of less than 65 • and a sensor 10 zenith angle of less than 41.4 • are used to reduce the impact of known biases (Grosvenor and Wood, 2014;Eastman and Wood, 2016;Grosvenor et al., 2018). Finally, only pixels with an inhomogeneity index (Cloud_Mask_SPI) of less than 30 are used to account for biases in the effective radius (r e ) in inhomogeneous scenes (Zhang and Platnick, 2011). Trials using a more stringent upper limit of 10 show little difference to the results presented here (not shown). The N d is gridded to a 1 • by 1 • resolution and finally, the condensation rate temperature correction from Gryspeerdt et al. (2016) is applied. 15 The MODIS LWP is gridded to a 1 • by 1 • resolution from MYD06_L2, selecting only liquid, single layer clouds with tops warmer than 268 K. The extra filtering applied to the N d is not applied to the LWP at the pixel resolution as the LWP is less sensitive to r e biases and this filtering would significantly bias the LWP against AMSR-E by selecting primarily high LWP scenes. However, only 1 • by 1 • gridboxes with a N d retrieval are retained for this analysis, resulting in an implicit filtering by satellite and solar zenith angles. 20 As both the MODIS LWP and N d rely on the adiabatic assumption and the same retrieved cloud properties, there is a significant potential for errors in these properties due to failures of the adiabatic assumption (Merk et al., 2015) and consequent correlated errors generating a N d -LWP relationship (E2b). The N d retrieval is better able to deal with non-adiabatic clouds than the effective radius retrieval alone (Painemal and Zuidema, 2011). For the majority of this work, the LWP is determined using V6 of the AMSR-E Ocean product (Wentz and Meissner, 2004), a passive microwave product that does not depend on 25 the adiabatic assumption. Clear-sky bias corrections are applied following Lebsock and Su (2014) at the pixel level. As the windspeed and sea surface temperature retrievals are unreliable in precipitating scenes, they are interpolated to precipitating locations by fitting a cubic mesh (Jones et al., 2001). To determine the in-cloud LWP, the AMSR-E LWP is divided by the MODIS cloud retrieval cloud fraction (CF) at the AMSR-E pixel level (14 km), with pixels having a CF of less than 10% being excluded due to the large uncertainty in the resulting in-cloud LWP. Finally, the AMSR-E data is gridded from the sensor 5 footprint of 14 km to a 1 • by 1 • resolution.
As a linear sensitivity ( d ln L d ln N d ) is not able to fully describe the non-linear relationship between N d and LWP, a piecewise relationship of the form (Eq.1) is used. L p and N p d are the LWP and N d values at the intersection between the two parts of the curve, while m l and m h are the gradients of the fit for the low and high N d portions of the curve. This curve is fit to the N d -LWP joint probability histogram (P(L|N d )), using the Levenberg-Marquardt algorithm in log-space (Jones et al., 2001). By fitting to 10 the joint probability histogram, each N d bin is given equal weight, rather than the weighting by the present day N d probability distribution implicit in the standard linear regression. This leads to a clearer picture of the overall form of the relationship, as the shape of the relationship does not change as the N d distribution changes (as might happen due to anthropogenic aerosol emissions; Gryspeerdt et al., 2017). Note that this method using "snapshots" of cloud fields, restricts the analysis to inferring information about cloud development, rather than studying their evolution directly (e.g. Matsui et al., 2006;Meskhidze et al., 15 2009; Gryspeerdt et al., 2014b).
To convert a change in LWP to a change in top of atmosphere radiation, data from the CERES 1 degree daily Single Scanner 20 Footprint, Edition 4 dataset is used (Wielicki et al., 1996). The all-sky albedo from CERES (α) is histogrammed as a function of the CF (f c ), LWP and N d , creating a single, global, joint probability histogram (P (α|f c , L, N d )). Given the retrieved cloud properties for a location (f c , LWP and N d ), this histogram produces a distribution of consistent values of the all-sky albedo (P(α)). This can be used to calculate the mean oceanic albedo to within 1% in the tropics, with an RMS error in the tropics of 1%, increasing to around 5% near the poles. These variations are primarily due to differences in the mean solar zenith angle 25 between the MODIS and CERES datasets, such that they have a small effect when determining the albedo sensitivities in this work.
Following Eq. 2, the N d -LWP and N d -f c relationships can be used to determine a change in scene/all-sky albedo as a function of an N d change. The relationships are treated as conditional probabilities (P(L|N d )= P (L,N d ) P (N d ) ), following Gryspeerdt et al. (2016). When combined with the N d sensitivity to aerosol (τ a ) changes P(N d |τ a ), this allows the scene albedo as a function of 30 aerosol (P(α|τ a )) to be calculated for a given scene of liquid clouds (Eq.3), where the circumflex indicates that a variable has been set to a certain value (the causal relationship). Note that this is different from the observed relationship P(α|τ a ), due to the confounding effects of local meteorology (Pearl, 1994;Gryspeerdt et al., 2016). It also makes the assumption that the observed conditional probabilities represent the causal relationship (i.e. P(L|N d )=P(L|N d ), representing only E1), an assumption that will be investigated in this work.
The albedo sensitivity to aerosol through modifications of each of the components of the albedo (N d , L, f c ) can be determined by replacing probabilities conditioned on N d with unconditional probabilities. For example, the sensitivity due only to N d variations (the Twomey effect; Twomey, 1974) can be determined by removing any dependence of CF and LWP on N d The change in planetary albedo is then determined by multiplying each gridbox by 1-f ice c (the ice cloud fraction), making the implicit assumption that there is no change in the ice cloud albedo or f ice c .

10
This is converted to a radiative forcing by multiplying by the incoming solar flux and anthropogenic aerosol fractions from Bellouin et al. (2013) and Kinne (2019).
To avoid uncertainties associated with the aerosol anthropogenic fraction inherent in estimates of the aerosol radiative forcing, the ERF due to LWP changes is not reported directly, only as a fraction of the RFaci calculated using the same dataset (Bellouin and Quaas, in p). The value for the forcing due to LWP changes can be re-constructed using an appropriate estimate 15 of the RFaci if required (e.g. Quaas et al., 2008;Stevens et al., 2016;McCoy et al., 2017;Gryspeerdt et al., 2017). A similar negative relationship is observed when using the AMSR-E LWP, both the all-sky LWP (Fig 2c) and the in-cloud LWP (Fig 2e). The relationship in Fig. 2c, using the all-sky LWP, is much weaker than the in-cloud LWP in Fig. 2e, which is 25 the most strongly negative linear sensitivity of the three relationships in Fig. 2. A strong positive relationship remains in the East China Sea.
The N d -LWP joint histograms shown in the right hand column of Fig. 2 show that the N d -LWP relationship is highly nonlinear at a global scale. All of the histograms show an increase in the LWP with increasing N d at low N d , followed by a decrease in the LWP at high N d . Despite global variations in N d and LWP retrieval biases (e.g. Grosvenor and Wood, 2014) and in N d , 30 this non-linearity is not obvious in the global plots of the linear sensitivity. However, a similar variation in the sensitivity   Xue and Feingold (2006) and Xue et al. (2008), where a high N d can result in a LWP reduction in clouds where precipitation does not reach the surface.
The differences between the fits of Eq.1 to the MODIS (Fig 2b) and the AMSR-E (Fig 2f) histograms demonstrate how a simple linear regression for calculating a sensitivity does not capture the strength or nature of the relationship. The AMSR-E relationship in Fig. 2f has a slightly weaker negative relationship at high N d (m h ) than that found using MODIS data (Fig. 2b

Regional relationships
Due to the difficulty of visualising joint histograms globally, and the sparse nature of the histograms in some regions making fitting Eq.1 error prone, a clustering approach is used to select regions with similar microphysics. A k-means clustering method (Anderberg, 1973) is used on the N d -LWP joint probability histograms representing each 1 • by 1 • gridbox. The algorithm is 20 modified to deal with missing data (k-POD; Chi et al., 2016), resulting in two distinct clusters over ocean with each gridbox being assigned to a single cluster (Fig.3). The clustering algorithm fills in missing data in the histograms with data interpolated from the clusters. This may reduce the number of clusters, but it suffices for demonstrating global variation. The first cluster  to N d (Fig. 3c) is similar to the results of Malavelle et al. (2017), who suggest a weak overall response of LWP to N d variations in a region where cluster two dominates. However, it means that the mid-latitude response may be a poor constraint on the response of the subtropical stratocumulus to N d perturbations, an issue that is of particular importance given the large role of the stratocumulus decks for the global aerosol forcing (Gryspeerdt and Stier, 2012).

5
While the overall form of the relationship remains the same, there is some variation in the N d -LWP joint histogram as a function of the meteorological state (Fig. 4). Following previous studies (Chen et al., 2014;Michibata et al., 2016), the data are separated by low troposphere stability (LTS) and relative humidity at 750 hPa (RH 750 ; approximately cloud top). Although the saturation deficit is more closely related to evaporation rates, we use RH 750 for consistency with previous work.
The response to LTS variations is small, occurring primarily in the part of N d -LWP space where precipitation is expected 10 (Figs.4c,f). The weak response to LTS is different from previous studies, which have shown a similar sized response to LTS and RH changes (Chen et al., 2014). A comparison between Figs. 4a,b shows that this variation in the linear sensitivity is partly due to variations in the N d distribution. At high LTS (Fig. 4b), the mean N d is larger than that found at low LTS (Fig. 4a), resulting in a more negative linear sensitivity. However, the high-N d sensitivity from the fitted relationship (m h ) is very similar at both high and low LTS. The difference in the precipitating region sensitivity (m l ) may be due to variations in the precipitation processes or regime dependent retrieval errors for shallow cumulus (low LTS) and stratocumulus clouds (high LTS). However, the low frequency of occurrence of these low-N d conditions (the histograms under each joint histogram in Fig. 4)  These changes in m h as a function of RH 750 and LTS fit the conclusions of previous studies (Ackerman et al., 2004;Chen et al., 2014;Michibata et al., 2016); increased entrainment at higher N d results in a reduction of the LWP, with a stronger decrease at lower cloud top humidities. Results using the saturation deficit are similar, but with an increased magnitude ((see S.I.)). The resulting decrease in LWP with increasing N d would reduce cloud albedo, offsetting the RFaci (also due to an increase in N d ) and reducing the overall ERFaci. In situations where there is a loop or feedback in the causal graph (e.g. Fig. 1), an experiment is required to determine the 20 strength of the causal relationship. Although the capability to artificially alter N d over a large spatial and temporal scale does not exist, large aerosol perturbations are able to alter the CCN environment and hence N d independently of any feedbacks or confounders (E2-4; Fig. 1). The N d -LWP relationship produced by these "natural experiments" would therefore be expected to be closer to the causal impact of aerosol on LWP than the relationship determined in Sec. 4.
(2011), the Kilauea volcano on the island of Hawai'i is used as an exogeneous aerosol perturbation. Previous work has shown a stronger linear AOD-N d sensitivity downwind of the Hawai'i than in surrounding regions, demonstrating the strong impact the SO 2 from Kilauea has on N d (Gryspeerdt and Stier, 2012). There is significant variability in the SO 2 emitted from the volcano.
Comparing a year with strong SO 2 emissions (2008) with a low emissions year (2007) shows that the variation in aerosol index 30 (AI; AOD times Ångström exponent; Nakajima et al., 2001) downwind from the volcano comes primarily from the variation in aerosol (Fig. 5a), rather than in meteorological conditions. Despite the strong negative N d -LWP relationship observed in sub-tropical regions (Fig. 3b), there is no change in the LWP ( Fig. 5b) in the region with a strong change in AI (region A). This lack of a LWP response to volcanic emissions is similar to the results of Malavelle et al. (2017), but is within the area covered by the more sensitive cluster (Fig. 3). The weak LWP response to aerosol variations suggests that the strong negative N d -LWP relationship (Figs. 2, 3) is unlikely to describe the impact of N d variations on LWP. and upwind (B; Fig. 5f) of the volcano, with a strongly negative m h and negative linear sensitivity. However, in the high aerosol environment of 2008 (Fig. 5c), this negative relationship becomes much weaker in the volcanic plume (m h =-0.15), whilst little change is observed upwind of the island (Fig. 5e). There is a difference in the LTS between 2007 and 2008 of around 1 K properties cannot explain the changes in region A. This means that the inter-annual difference in region A can be attributed primarily to aerosol variations (E1).
In the absence of feedbacks (E3), additional confounders (E4) and meteorological variations, the N d -LWP relationship should be insensitive to the cause of the N d variations. Given the similarity in the meteorological conditions between the years, the difference in the N d -LWP relationship in region A therefore suggests that the relationship is modified by feedbacks (E3) 5 or additional confounders (E4). Due to the high volcanic emissions, the 2008 N d -LWP relationship in region A is known to be strongly controlled by aerosol variations (E1) and has a reduced impact of other processes (E2-4), such that it is likely closer to the causal N d -LWP relationship. This indicates a considerably weaker role for N d than determined in Sec. 4. With an m h of -0.15, the in-plume results are much closer to the results from LES simulations (Ackerman et al., 2004;Bretherton et al., 2007;Xue et al., 2008, m h < −0.2) and in-situ observations of shiptracks, where decreases in LWP have been observed in particularly 10 polluted conditions (Ackerman et al., 2000;Noone et al., 2000;Christensen and Stephens, 2011;Goren and Rosenfeld, 2014).
The consequently weaker LWP response to aerosol is in better agreement with the weak LWP changes observed in Fig. 5b and Malavelle et al. (2017). The Kilauea volcano affects primarily shallow cumulus clouds (Oreopoulos et al., 2014), which exert a weak control on the ERFaci from LWP changes due to their low liquid CF. The processes responsible for a reduction in LWP (E1c) may be different 15 from those controlling stratocumulus clouds (E1b). Shipping provides another source of exogeneous aerosol perturbations , generating shiptracks that are primarily concentrated in the high CF stratocumulus regions. Using a database of shiptracks from Christensen et al. (2014), the relationship between the in-shiptrack N d and LWP increase in the shiptrack compared to the control region around the track (dLWP) indicates how the LWP responds to N d perturbations. As the N d values always increase from the control region to the inside the shiptrack, dLWP shares a sign with the gradient of the N d -LWP relationship. Note that due to the required spatial resolution, the LWP for these shiptracks is retrieved using MODIS, rather than AMSR-E.
For low control values of the LWP (Fig. 6), increases in LWP (positive values of dLWP) are seen at lower in-shiptrack values 5 of N d , but as the shiptrack N d gets higher, the dLWP reduces to close to zero, with a negative dLWP for the most polluted cases.
When the control LWP is high, dLWP is consistently weakly negative, although this likely is due to regression to the mean effects (the mean control LWP is 82 gm −2 ). This suggests that the LWP becomes insensitive to further aerosol/N d perturbations once the LWP reaches a sufficient magnitude, consistent with an aerosol suppression of precipitation (E1a). These small dLWP values at high N d are consistent with the Kilauea results, suggesting a weak LWP response at high N d . If the LWP response in 10 shiptracks followed the relationships from Sec. 4, a strong negative dLWP should be visible at high N d , in contrast to the weak negative response actually observed (Fig. 6).
By selecting situations where aerosol is known to be responsible for N d variations (so-called "natural experiments"), the impact of feedbacks (E3) and additional covariations (E4) can be reduced (although not completely removed). In these situations, the N d variations are driven by exogeneous aerosol perturbations, such that the LWP variations are a response to (rather than 15 a driver or indicator of) the change in N d (E1 only). This means that the N d -LWP relationship during these "natural experiments" provides better information on the LWP response to N d variations, such that the strong negative N d -LWP relationships observed in Sec. 4 likely overestimate the decrease in LWP in response to aerosol perturbations. While the satellite-derived relationships may therefore be unsuitable as a direct estimate on the aerosol impact on LWP, they could be used as a lower bound on the LWP change (an upper bound on the radiative forcing) from aerosol-induced LWP decreases.

The implied ERFaci
The planetary albedo sensitivities to aerosol perturbations are shown in Fig.7 following Eq. 2. Due to the difficulty of visualising joint histograms globally, linear sensitivities are determined from the joint histograms (P(α|τ a )) by weighting by the present day aerosol distribution (see Gryspeerdt et al., 2016). The first three subplots show the albedo sensitivity through modifying the N d (constant CF and LWP; Fig. 7a), CF (constant LWP; Fig. 7b) and AMSR-E LWP (constant CF; Fig. 7c). Both changes 25 in N d and CF increase the scene albedo, which results in a negative radiative forcing. They have somewhat different spatial patterns, with the albedo sensitivity to N d changes being concentrated in the centres of the stratocumulus decks due to the high liquid cloud fraction. The sensitivity to CF changes is highest at the edges of the stratocumulus decks, where the greatest potential for modifying the cloud fraction exists, as found in previous studies (Gryspeerdt et al., 2016;Christensen et al., 2017;Andersen et al., 2017). 30 The sensitivity to LWP changes is also strongly dependent on the liquid CF and so is strongest in the centres of the stratocumulus decks (Fig. 7c). As a reduction in LWP with increasing N d is observed in these regions (Fig. 3) albedo sensitivities in Fig. 7 with the anthropogenic aerosol fraction from Bellouin et al. (2013) implies a positive radiative forcing from LWP changes that offsets 62% of the RFaci calculated using the same data, resulting in a weakening of the RFaci.
The offset is similar (59%) when using the anthropogenic fraction from Kinne (2019). This is likely the upper bound on the fraction of the RFaci offset by LWP reductions, following the results of Sec. 5 and supported by the weaker offsetting in regions with larger aerosol perturbations (e.g. the East China sea, the tropical and north Atlantic). Despite the reduced albedo 5 sensitivity due to the LWP reduction, the overall albedo sensitivity to aerosols is still positive (Fig. 7d), resulting in a negative ERFaci from liquid clouds due to the strong implied forcing from the N d -CF relationship (approximately a 200% increase above the RFaci).
There remains considerable uncertainty in the magnitudes of these effects. The albedo change is only calculated over ocean.
Observational studies suggest the N d change and RFaci over land are small, but it is possible that the LWP adjustments could have a very different character and relationship to the RFaci over land. The variation in the N d -LWP relationship in the Kilauea volcanic plume (Fig. 5) and the response of the LWP in shiptracks (Fig. 6), suggest that the LWP change determined in Fig. 7 is overly strong. This would then place a 62% offset of the RFaci as the upper bound on the radiative forcing from LWP changes 5 (larger offsets are unlikely). This is consistent with previous work, where an increase in cloud albedo is found in response to a change in aerosol (Lebsock et al., 2008;Chen et al., 2014;Christensen et al., 2016), such that a LWP reduction cannot completely offset the RFaci.

Discussion
This work demonstrates that a non-linear relationship exists between N d and LWP (Fig. 2). These results are in agreement with 10 previous studies, with an increase in LWP with N d at low N d from precipitation suppression (E1a), but a decrease at high N d due to increased cloud top or lateral entrainment (E1b, c). The similarity in the relationship when using different measures of LWP suggests that this relationship is not primarily due to LWP retrieval errors (E2). There are global variations in the N d -LWP relationship and significant changes accompany variations in meteorological factors, particularly RH 750 (Fig.4). The observed N d -LWP relationship implies a reduction in LWP with increasing aerosol and N d , resulting in a positive radiative forcing that 15 offsets around 60% of the RFaci.
The analysis in Sec. 5 suggests that the negative N d -LWP relationship observed over much of the world may be overestimated, resulting in too strong a corresponding positive radiative forcing due to aerosol induced LWP adjustments. A precipitation feedback (E3a) would produce a positive N d -LWP relationship and so is unlikely to be responsible. An entrainment-based feedback on the N d (E3b) or an additional confounder (E4) could be responsible for the negative N d -LWP relationship.

20
The albedo sensitivity to aerosol via LWP changes is particularly strong in the stratocumulus regions (Fig. 7), due to the high liquid cloud fraction. This implies an important role for the sedimentation-entrainment feedback (E1b). With the entrainment of dry environmental air at the cloud top, the assumptions in the N d retrieval of a linearly increasing liquid water content and vertically constant N d no longer hold as the cloud is no longer adiabatic, such that the cloud top N d is no longer representative of the cloud base N d . A reduction in the cloud top r e by homogeneous mixing during entrainment would produce an increase in However, although some studies have found evidence of homogeneous mixing in stratocumulus cloud (Breon and Doutriaux-Boucher, 2005;Yum et al., 2015), many studies have found that inhomogeneous mixing dominates, particularly at cloud 30 top (Nicholls, 1987;Pawlowska et al., 2000;Gerber et al., 2005;Burnet and Brenguier, 2007;Yum et al., 2015). While inhomogeneous mixing reduces the N d , in extreme cases it does not result in a r e change, so may not be detected by satellite.
As such, some proportion of homogeneous mixing is required for E3b to generate a negative N d -LWP relationship in satellite data. A discrepancy between satellite retrieved and in-situ N d as a function of humidity or entrainment rate might be one indicator of this process. Further investigation into the mixing and behaviour of these retrievals at cloud top is necessary to establish the impact of E3b on the N d retrievals and the N d -LWP relationship.
An additional, unknown confounder (E4) is also a possible explanation for the results in Sec. 5. This effect would have to act on both N d and LWP together -A process that only affects one would not generate the systematic bias required. Even if 5 such an unknown, additional confounding process exists, the conclusion drawn from Sec. 5 would still hold -that the implied aerosol impact on LWP in Fig. 7 is likely too strong.
By using 1 • by 1 • average values, this work ignores the impact of sub-grid variability of the N d and LWP retrievals (Zhang et al., 2018). Preliminary work indicates that this may modify the relationship, with the strength of the relationship changing when it is determined at smaller spatial and temporal scales. If the interpretation of the results from natural experiments is 10 followed, it implies that these small scale N d -LWP relationships are strongly influenced by E2-4, due to the lack of aerosol variation to drive the N d variation necessary to highlight the impact of E1. The cause of this scale dependence will be investigated in future studies.
Although volcanic emissions (Fig. 5) and shiptracks are exogeneous sources of aerosol, the datasets linked to these sources are limited. They occur in relatively restricted locations on the globe and there are a small number of the high N d retrievals 15 required to populate the N d -LWP histogram (Fig.6). While the shiptrack dataset is concentrated in stratocumulus region , it is still possible that the effect on shallow cumulus clouds could be large enough to overcome the relatively small CF in this regime which has previously been shown to restrict the contribution of shallow cumulus clouds to the RFaci (Gryspeerdt and Stier, 2012). Given the importance of the N d to this work, an improved understanding of the behaviour of the N d retrieval through a comparison with in-situ data is particularly important. Future studies are planned to expand this dataset 20 of exogeneous aerosol perturbations in marine clouds such that a more representative global study of this type can be performed. Process-resolving simulations of these cases and a comparison to the global results are necessary to fully understand the behaviour of the satellite retrievals and how accurately they can represent the aerosol-N d -LWP system to better constrain the aerosol impact on LWP. which has been previously suggested to be due to the droplet size impact on entrainment (E1b/c, Fig.2). This non-linearity of the N d -LWP relationship restricts the ability of linear regressions to characterising the relationship. The reduction in LWP with increasing N d is only slightly stronger when using MODIS LWP compared to the in-cloud LWP from AMSR-E, suggesting that although correlated errors in the MODIS LWP and N d can play a role (E2), they do not dominate the magnitude of the N d -LWP relationship.
By clustering the N d -LWP joint histograms, it is shown that the primary variation in the histograms comes from variations 5 in the LWP behaviour at high N d (Fig. 3). In the subtropical subsidence regions, there is a clear LWP reduction with increasing N d , whilst in other regions, LWP remains constant or even increases with LWP even at high N d . The global relationship is dominated by the subtropical relationship due to the high liquid CF and higher N d variation in these regions, but the regional variations in the N d -LWP relationship make it difficult to use the results from one region to constrain others.
Part of this variability come from regional differences in meteorological conditions. Significant variations in the N d -LWP  The observed N d -LWP relationship suggests that LWP adjustments could offset up to 60% of the RFaci/Twomey effect ( Fig. 7), as a positive radiative forcing. This represents an upper bound on the positive radiative forcing expected from a LWP reduction. The results from natural experiments suggest that the LWP response is likely weaker than this (Figs. 5, 6), as the causal N d -LWP relationship is obscured by feedbacks (E3) and additional confounders (E4) in many cases. Further work is required to bound the LWP response, but these results suggest that the overall ERFaci is likely to be negative, supported by 25 previous studies that have found a complete offset of the RFaci is unlikely (Chen et al., 2014).
Although it has been demonstrated in this work that the N d -LWP relationship has a substantial impact on the ERFaci, it is clear that significant uncertainties remain. The satellite retrieved N d -LWP relationship has several features that are similar to the relationship predicted by high resolution models (Ackerman et al., 2004;Sato et al., 2018), but it is not clear the extent to which these relationships represent the causal relationship and so can be used to constrain aerosol-cloud interactions. A wider 30 study of the effect of aerosols on LWP due to exogenous aerosol perturbations in a variety of cloud regimes would provide one avenue for progress, as would finding a suitable mediating variable within the N d -LWP relationship.
Author contributions. All of the authors participated in the design of the study. EG performed the analysis and wrote the paper. MC provided the shiptrack dataset. All of the authors assisted in the interpretation of the results and commented on the paper.

Appendix A: Expected sensitivities
If the LWP and N d are calculated from MODIS data using the adiabatic assumption (Wood, 2006;Quaas et al., 2006), they take the form where 0 < f ad ≤ 1 is the adiabatic factor (f ad =1 is completely adiabatic) and c(T) is the temperature correction to the 10 condensation rate from Gryspeerdt et al. (2016). By similar logic, the sensitivity expected at a constant r e from variations in τ c is d ln L d ln N d re = 2 (A8) 20 Note that the cause of these variations is not specified. A variation in r e due to retrieval errors or N d variations would produce the same effect. As both the LWP and N d relate to the adiabatic factor in the same way as the optical depth, the expected sensitivity from adiabatic factor variation is also 2. Andersen, H., Cermak, J., Fuchs, J., Knutti, R., and Lohmann, U.: Understanding the drivers of marine liquid-water cloud occurrence and