POLDER-3 aerosol retrievals

In this work, we use the POLDER-3 aerosol product retrieved by the SRON aerosol retrieval algorithm16,17 (processed for the year 2006) previously used for computing the direct radiative effect of aerosols18. Retrievals are based on Collection-3 level-1 POLDER-3 medium resolution data (18 \(\times\) 18 km\({}^{2}\)). We restrict ourselves to retrievals over ocean because of limited quality of over-land retrievals, especially for cases with low aerosol loading. Furthermore, we restrict the study to \(6{0}^{\circ }\)S \(<\ {\rm{latitude}}\ <6{0}^{\circ }\)N. POLDER-3 achieves global coverage in 1.5 days.

Cloud screening has been performed based on goodness-of-fit between forward model and measurement17. This means aerosol retrievals in (partly) cloudy scenes and directly next to clouds are excluded33. The aerosol products have been gridded on a 1\({}^{\circ }\times {1}^{\circ }\) latitude-longitude grid. The aerosol retrieval algorithm is based on a description of aerosols by a fine and a coarse mode (indicated by superscripts \(f\) and \(c\), respectively), where each mode is described by a log-normal function. Aerosol properties in the state vector are for both modes the effective radius \({r}_{{\rm{eff}}}\), effective variance \({v}_{{\rm{eff}}}\), the real and imaginary part of the refractive index \({m}_{{\rm{r}}}\) and \({m}_{{\rm{i}}}\), and the aerosol column number concentration \({N}_{a}\). The fraction of spheres34 \({f}_{{\rm{sp}}}\) of the coarse mode is included as a fit parameter in the retrieval state vector. In this study, we use the AOD, \({N}_{{\rm{a}}}\), \({r}_{{\rm{eff}}}\), and \({v}_{{\rm{eff}}}\) of the fine and coarse mode, and \({f}_{{\rm{sp}}}\).

\({N}_{{\rm{ccn}}}\) is computed from the log-normal bi-modal size distribution for the retrieved \({N}_{a}\), \({r}_{{\rm{eff}}}\) and \({v}_{{\rm{eff}}}\) of the fine and coarse mode as the column number of particles with radius \(> \, {r}_{\mathrm{lim}}=0.15\) µm. To investigate the capability of POLDER-3 to provide \({N}_{{\rm{ccn}}}\), we created 1000 synthetic POLDER-3 measurements with varying aerosol properties in the following range (superscripts \(f\) and \(c\) indicate fine and coarse mode, respectively): 0.02–0.3 µm for \({r}_{{\rm{eff}}}^{f}\), 0.65–3.5 µm for \({r}_{{\rm{eff}}}^{c}\), 0.1–0.3 for \({v}_{{\rm{eff}}}^{f}\), 0.4–0.6 for \({v}_{{\rm{eff}}}^{c}\), 1.33–1.60 for \({m}_{r}^{f}\) and \({m}_{r}^{c}\), 10\({}^{-8}\)−0.1 for \({m}_{i}^{f}\), 10\({}^{-8}\)–0.02 for \({m}_{i}^{c}\), 0.005–0.7 for AOD\({}^{f}\) and AOD\({}^{c}\), and 0–1 for \({f}_{{\rm{sp}}}\). We put a random error on the synthetic measurements with a standard deviation of 2% for radiance and 0.015 for degree of polarization, which is representative for POLDER-3 measurements over ocean35. From the synthetic experiment we conclude that the uncertainty on individual retrievals of \({r}_{{\rm{eff}}}^{f}\) is 0.034 µm (bias 0.016 µm) and on \({f}_{{\rm{sphere}}}\) 0.25 (bias 0.13). For \({N}_{{\rm{ccn}}}\) we find that for 71% of the data the difference between retrieved and true \({N}_{{\rm{ccn}}}\) is smaller than \(0.20\cdot {N}_{{\rm{ccn}}}+4\,\times1{0}^{6}\) cm\({}^{-2}\). We use this as an uncertainty estimate of \({N}_{{\rm{ccn}}}\).

Comparing POLDER \({N}_{{\rm{ccn}}}^{{\rm{pol}}}\) with corresponding values \({N}_{{\rm{ccn}}}^{{\rm{aer}}}\) computed from the aerosol size distribution of ground-based aerosol robotic network (AERONET) measurements, we find (Supplementary Fig. 1) \({R}^{2}=0.58\), a bias of 8.0 × 10\({}^{6}\) cm\({}^{-2}\), a mean absolute difference of 1.95 × 10\({}^{7}\) cm\({}^{-2}\), and a root mean square difference of 3.90 × 10\({}^{7}\) cm\({}^{-2}\). Totally, 51% of the POLDER-AERONET differences are smaller than the error bound found from the synthetic experiment (\(0.20\cdot {N}_{{\rm{ccn}}}+4\times 1{0}^{6}\) cm\({}^{-2}\)), which suggest that this is a reasonable error estimate given that the differences are also affected by errors in AERONET data. Most important, there is no significant trend of the relative difference \((({N}_{{\rm{ccn}}}^{{\rm{pol}}}-{N}_{{\rm{ccn}}}^{{\rm{aer}}})/{N}_{{\rm{ccn}}}^{{\rm{pol}}})\) with \({N}_{{\rm{ccn}}}^{{\rm{pol}}}\) (\({R}^{2}=0.01\)), which would affect the \({N}_{{\rm{ccn}}}-{N}_{{\rm{d}}}\) relationships.

MODIS cloud retrievals

We use the MODIS Collection-6 retrievals of cloud effective radius (CER) and cloud optical thickness (COT)19 to compute the cloud droplet number concentration \({N}_{{\rm{d}}}\) using the adiabatic approximation11,36, and aggregate the data at a 1\({}^{\circ }\times {1}^{\circ }\) horizontal and daily temporal resolution. Here, we consider only points for which CER \(> \) 4 µm and COT \(> \) 4, cloud fraction \(> \) 0.9 (at 5 km resolution), and with a sub-pixel inhomogeneity index (cloud_mask_SPI) \(<\) 30, as these are the ranges where meaningful CER and COT retrievals can be performed36. Further, we only consider liquid clouds by selecting 1\({}^{\circ }\times {1}^{\circ }\) grid cells with COT_ice = 0.

Deriving susceptibility

To derive the susceptibility from POLDER-3 and MODIS data, we use all grid cells for the year 2006 for which we have both a valid POLDER-3 \({N}_{{\rm{ccn}}}\) value and a valid MODIS \({N}_{{\rm{d}}}\) value. After selecting the data points in the region under consideration, we define a number of bins for \({N}_{{\rm{ccn}}}\) where each bin has an equal number of points. For deriving susceptibility we use 100 bins. The \({N}_{{\rm{ccn}}}\) value attributed to a certain bin is the median of all \({N}_{{\rm{ccn}}}\) values in that bin. The \({N}_{{\rm{d}}}\) value attributed to this bin is the median of all \({N}_{{\rm{d}}}\) values collocated with the \({N}_{{\rm{ccn}}}\) retrievals in that bin. This procedure gives us 100 paired values of \({N}_{{\rm{ccn}}}\) and \({N}_{{\rm{d}}}\) for the region under consideration. So, for determining the global (ocean) value for susceptibility we aggregate all valid POLDER3-MODIS data pairs into 100 bins. Having the same number of measurements in each bin ensures that each bin has the same statistical representativity. The susceptibility is determined by fitting a linear regression through the the binned points of \({N}_{{\rm{d}}}\) versus \({N}_{{\rm{ccn}}}\). The derived susceptibilities do not change significantly if a different number of bins is chosen (e.g., for 20 or 1000 bins the global susceptibilities agree within 0.01 with those for 100 bins).

An important assumption we make here is that the aerosol properties are uniformly distributed over a grid cell and that relative variations in the column integrated aerosol concentration are representative for the aerosol that interacts with the cloud at cloud base. Under this assumption, it is unnecessary to perform retrievals below clouds or directly next to clouds, which is not possible with current retrieval algorithms33. These assumptions may lead to an underestimation of susceptibility because if the retrieved aerosol is not representative for the aerosol that interacts with clouds, variations in retrieved aerosol properties would have no effect on \({N}_{{\rm{d}}}\). Aspects related to spatial sampling by the satellite are not expected to impact the derived susceptibility12.

Simulating the effect of measurement errors

To investigate the effect of measurement errors on \({N}_{{\rm{ccn}}}\), we use a simulator that models \({N}_{{\rm{d}}}\) as function of \({N}_{{\rm{ccn}}}\). For the ensemble of \({N}_{{\rm{ccn}}}\) we use the values measured by POLDER and we compute the corresponding \({N}_{{\rm{d}}}\) assuming \({N}_{{\rm{d}}}=C\ {N}_{{\rm{ccn}}}^{S}\), with S = 0.66 and C = 0.001. When we apply our procedure to determine susceptibility (see above) to the simulated data set we find exactly S = 0.66, as expected. When we put 50% error (1\(\sigma\)) error on \({N}_{{\rm{d}}}\) we still find S very close (within 0.01) to 0.66, showing that derived susceptibilities are not very sensitive to random fluctuations on \({N}_{{\rm{d}}}\). Next, we put random errors on \({N}_{{\rm{ccn}}}\) of the form \({{\rm{err}}}_{{\rm{ccn}}}\) = rel\({}_{{\rm{err}}}\)\({N}_{{\rm{ccn}}}\) + abs\({}_{{\rm{err}}}\), where rel\({}_{{\rm{err}}}\) and abs\({}_{{\rm{err}}}\) denote the relative- and absolute error on \({N}_{{\rm{ccn}}}\), respectively. We choose a range (1\(\sigma\)) for rel\({}_{{\rm{err}}}\) of 0–0.5 and abs\({}_{{\rm{err}}}\) of 0–5 × \(1{0}^{6}{\rm{cm}}^{-2}\). From the simulated data with error we determine the susceptibility \({S}_{{\rm{full}}}\) using the whole range of \({N}_{{\rm{ccn}}}\) and \({S}_{{\rm{opt}}}\) using only values \(> \) 10\({}^{7}\) cm\({}^{-2}\), leaving out 15% of the smallest values which are most heavily affected by the absolute term in the measurement uncertainty.

Supplementary Fig. 2 shows \({S}_{{\rm{full}}}\) (left panel) and \({S}_{{\rm{opt}}}\) (right panel) as function of abs\({}_{{\rm{err}}}\) for different values of rel\({}_{{\rm{err}}}\). It can be seen that \({S}_{{\rm{full}}}\) is being under-estimated in the presence of measurement errors, as is expected from earlier results based on AOD12. The values for \({S}_{{\rm{opt}}}\) are much closer to the true value of 0.66, although under-estimation is still possible (up to 0.18 for the range of \({N}_{{\rm{ccn}}}\) errors shown). The simulation results confirm that \({S}_{{\rm{opt}}}\) is a better estimate of susceptibility than \({S}_{{\rm{full}}}\). The expected measurement uncertainty of \(0.20\times{N}_{{\rm{ccn}}}+4\times1{0}^{6}\) cm\({}^{-2}\) indicates an underestimation in \({S}_{{\rm{opt}}}\) of 0.06. Choosing a larger cut-off value for \({N}_{{\rm{ccn}}}\) than 10\({}^{7}\) cm\({}^{-2}\) does not reduce the uncertainty range.

Radiative forcing calculation

We use five different global aerosol climate models (ECHAM6-HAM, CAM5.3, CAM5.3-CLUBB, CAM5.3-CLUBB-MG2, and SPRINTARS) to compute the increase \(\Delta {N}_{{\rm{ccn}}}\) (using the column CCN at 0.3% supersaturation) between PI and PD from a pair of nudged simulations that are the same except that the PI simulation uses pre-industrial, and the PD simulation, PD aerosol emissions. All these models, which were also used in the study by Gryspeerdt et al.8, have participated in the AEROCOM intercomparisons in current or previous model versions37,38. From \(\Delta {N}_{{\rm{ccn}}}\), for the different models, we compute \(\Delta {N}_{{\rm{d}}}\) using the values for \(S\) as derived from POLDER-3 and MODIS. From \(\Delta {N}_{{\rm{d}}}\) we compute the change in cloud albedo using the Twomey formula1 and the RF\({}_{{\rm{aci}}}\) using Eq. (3) of Gryspeerdt et al.8

$${{\rm{RF}}}_{{\rm{aci}}}=-{F}^{\downarrow }\ {f}_{{\rm{liq}}}{\alpha }_{{\rm{cld}}}(1-{\alpha }_{{\rm{cld}}})\,\frac{1}{3}\frac{\Delta {N}_{{\rm{d}}}}{{N}_{{\rm{d}}}},$$ (1)

where \({F}^{\downarrow }\) is the daily-mean incoming solar radiation flux at each grid-point for each day, \({f}_{{\rm{liq}}}\) the fractional cover by liquid–water clouds, and \({\alpha }_{{\rm{cld}}}\) the cloud albedo taken from CERES. Since we only derived values for susceptibility over ocean, the procedure above gives an estimate RF\({}_{{\rm{aci}},{\rm{ocean}}}\) = −0.76 W m\({}^{{\rm{-2}}}\) over the ocean and a range between −0.68 and −0.85 W m\({}^{{\rm{-2}}}\) from the range in \(\Delta {N}_{{\rm{ccn}}}\) from the different models. The procedure for computing RF\({}_{{\rm{aci}},{\rm{aod}}}\) and RF\({}_{{\rm{aci}},{\rm{ai}}}\) is the same as for \({N}_{{\rm{ccn}}}\) except that they are based on the PD–PI increase in AOD and AI, respectively.

To get an estimate of the land contribution to RF\({}_{{\rm{aci}}}\) we looked at the ratio RF\({}_{{\rm{aci}}}\)/RF\({}_{{\rm{aci}},{\rm{ocean}}}\) in 13 different aerosol climate models23. This gives us a range of values for RF\({}_{{\rm{aci}}}\)/RF\({}_{{\rm{aci}},{\rm{ocean}}}\) between 1.12 and 2.24, and a mean value of 1.5. Further, the simulator results indicate that for the expected error on \({N}_{{\rm{ccn}}}\) (\(0.20\;\times{N}_{{\rm{ccn}}}+4\;\times1{0}^{6}{\rm{cm}}^{-2}\)) the global value of \({S}_{{\rm{ccn}}}^{{\rm{opt}}}\) is underestimated by 10%, but given that this error estimate itself is quite uncertain also smaller and larger underestimations cannot be ruled out. Therefore, we model this uncertainty with a normal distribution of the RF\({}_{{\rm{aci}}}\) scaling factor with mean 1.05 and standard deviation 0.05, assuming a linear dependence of RF\({}_{{\rm{aci}}}\) on susceptibility. Combining this distribution with the different values for RF\({}_{{\rm{aci}},{\rm{ocean}}}\) and the different values of the ratio RF\({}_{{\rm{aci}}}\)/RF\({}_{{\rm{aci}},{\rm{ocean}}}\), we get a histogram of possible values for RF\({}_{{\rm{aci}}}\) We take the median of this this distribution, −1.14 W m\({}^{{\rm{-2}}}\), as our best RF\({}_{{\rm{aci}}}\) estimate, and define an uncertainty range using the 5 and 95 percentile values, respectively, which yield a range between −0.84 and −1.72 W m\({}^{{\rm{-2}}}\). If we ignore the uncertainty on \({S}_{{\rm{ccn}}}^{{\rm{opt}}}\) we obtain a range for RF\({}_{{\rm{aci}}}\) between −0.80 and −1.60 W m\({}^{{\rm{-2}}}\) and if we assume an error on \({S}_{{\rm{ccn}}}^{{\rm{opt}}}\) that is twice as large (mean 1.10 and standard deviation 0.10) we obtain a range between −0.86 and −1.84 W m\({}^{{\rm{-2}}}\). If we would only take into account the uncertainty in the ratio RF\({}_{{\rm{aci}}}\)/RF\({}_{{\rm{aci}},{\rm{ocean}}}\), the resulting RF\({}_{{\rm{aci}}}\) range would be from −0.85 to −1.70 W m\({}^{{\rm{-2}}}\). So, by far the largest part of the total uncertainty range can be explained by the uncertainty in the ratio RF\({}_{{\rm{aci}}}\)/RF\({}_{{\rm{aci}},{\rm{ocean}}}\) and the uncertainty caused by the uncertainty on \({S}_{{\rm{ccn}}}^{{\rm{opt}}}\) is small compared to the total uncertainty.

When computing the uncertainty ranges for RF\({}_{{\rm{aci}},{\rm{aod}}}\) (based on \({S}_{{\rm{aod}}}^{{\rm{full}}}\)) and RF\({}_{{\rm{aci}},{\rm{ai}}}\) (based on \({S}_{{\rm{ai}}}^{{\rm{full}}}\)) we do not take the uncertainty in S into account, in order to make these estimate comparable to previous studies.