Increased aerosols can reverse Twomey effect in water clouds through radiative pathway

Aerosols play important roles in modulations of cloud properties and hydrological cycle by decreasing the size of cloud droplets with the increase of aerosols under the condition of fixed liquid water path, which is known as the first aerosol indirect effect or Twomey-effect or microphysical effect. Using high-quality aerosol data from surface observations and statistically decoupling the influence of meteorological factors, we show that highly loaded aerosols can counter this microphysical effect through the radiative effect to result both the decrease and increase of cloud droplet size depending on liquid water path in water clouds. The radiative effect due to increased aerosols reduces the moisture content, but increases the atmospheric stability at higher altitudes, generating conditions favorable for cloud top entrainment and cloud droplet coalescence. Such radiatively driven cloud droplet coalescence process is relatively stronger in thicker clouds to counter relatively weaker microphysical effect, resulting the increase of cloud droplet size with the increase of aerosol loading; and vice-versa in thinner clouds. Overall, the study suggests the prevalence of both negative and positive relationships between cloud droplet size and aerosol loading in highly polluted regions.


Results
AERONET provides aerosol data generated from direct sun measurements 37 and almucantar measurements 38 . As direct sun measurement-based products have higher temporal resolution than almucantar measurement-based products, data counts are higher in the former than in the latter. Taking this advantage, aerosol data obtained from direct sun measurements are used to understand aerosol effects on cloud droplet effective radius (CER) and cloud optical thickness (COT). Aerosol and cloud data are binned for LWP with a bin width of 10 g/m 2 to determine ∂CER´/∂AOT´ and ∂COT´/∂AOT´ for each LWP bin after decoupling the effects of meteorological factors through multiple linear regression analysis 39 (see "Methods"). Here, x´ (x is CER, COT, or AOT) is the normalized anomaly (see "Methods"). Figure 1a,b show ∂CER´/∂AOT´ for CER corresponding to wavelength of 3.7 µm (MODIS band 20) and ∂COT´/∂AOT´ for COT corresponding to wavelength of 0.645 µm (MODIS band 1), respectively, for LWP bins having sample count greater than 30 for the regression analysis. In both of them, AERONET AOT is for 0.5 µm. The significance tests for regression models used to determine ∂ CER  Table 1. As suggested by Table 1, the regression models for different LWP bins are significant with more than 99.5% confidence level (α = 0.005), except for three cases for Gandhi College site: ∂ CER for 70 g/m 2 ≤ LWP < 80 g/m 2 (p-value = 0.01). Overall speaking, all results, except for ∂ CER ′ 3.7 /∂AOT ′ 0.5 for 0 g/m 2 ≤ LWP < 10 g/m 2 , are within 95% confidence interval. Similarly, as shown in Table 1, the r value ranges from 0.21 to 0.73 in different cases. In general, the r value increases with the increase of LWP. As the results shown in Fig. 1 are derived by considering two meteorological factors-lower tropospheric stability (LTS) and precipitable water content(PWC)-in the regression model (see "Methods"), a sensitivity study is performed by increasing these meteorological factors in the regression model. In the sensitivity study, ∂CER   Fig. 1. As 3.7 µm wavelength has a weaker cloud-depth penetration compared to shorter wavelengths (e.g., 1.6, 2.1 µm) used in cloud remote sensing, this wavelength provides information of droplet size for upper cloud layers. Figure 1a suggests that the response of increased AOT on CER (3.7 µm) is different for relatively thin (low LWP) and thick (high LWP) clouds: the thin (LWP < ~ 25 g/m 2 ) and thick (LWP > ~ 25 g/m 2 ) clouds clearly show the negative and positive relationships, respectively. Figure 1a further reveals that the positive relationship becomes stronger with the increase of LWP. On the other hand, as shown in Supplementary Fig. S2, CER corresponding to 1.6 µm, which represents the droplet size for cloud layers deeper than that for 3.7 µm, hardly shows such increase of the strength of positive relationship with the increase of LWP. The p-value and r value for the regression models corresponding to results of Supplementary Fig. S2 are given in Supplementary  Table S1. In general, aerosols near the cloud base influence the lower and/or middle cloud layers through direct aerosol-cloud interaction process, and this effect gradually expands towards the cloud top 10,12 . Such difference regarding LWP dependent strength for positive relationship between CER(3.7 µm) and CER(1.6 µm), as noted in Figs. 1a and S2, indicates that the aerosol-cloud interaction process initiated by aerosols below the cloud base can be opposed by the next process occurring near the cloud top (see "Discussion"). Further, by coinciding with Fig. 1a,b suggests the opposite response of increased aerosols on COT with respect to that on CER. For fixed LWP, such opposite response is not surprising as LWP can be approximated as k × COT × CER, where k is a constant term depending on the vertical inhomogeneity of cloud layers 40,41 . In other words, decreased (increased) CER can increase (decrease) cloud droplet concentration to increase (decrease) total cross section area and COT when LWP remains unchanged. Therefore, the positive relationship is associated with not only CER increment, but also COT decrement. Though a large volume of studies discussed CER increment with the increase of aerosols; however, to the best of our knowledge, this consequent response of increased aerosols on COT decrement has not been discussed in detail. The most plausible explanation for such COT decrement along with CER increment is cloud droplet coalescence process. This suggests that cloud droplet coalescence process can cause the positive relationship between CER and aerosol loading.    www.nature.com/scientificreports/ In order to get better insight into such cloud droplet coalescence process over our study regions, we analyzed the difference in CER between 3.7 µm and 1.6 µm (CER 3.7 -CER 1.6 ). It is because CER corresponding to 3.7 µm becomes smaller than that for 2.1 µm or 1.6 µm wavelength, if coalescence growth process dominates over condensation growth process, and vice-versa 42,43 . As shown in Fig. 2, the mean values of CER 3.7 -CER 1.6 are negative for all LWP bins, suggesting the dominance of coalescence growth process over both study regions. Figure 2 further shows that the CER 3.7 -CER 1.6 becomes more negative with the increase of LWP. This suggests that coalescence growth process becomes stronger in thicker clouds than in thinner clouds. It is because thicker clouds provide sufficient time and opportunity than thinner clouds for cloud droplets to collide. This behavior is important for the occurrences of negative and positive relationships in relatively thinner and thicker clouds, respectively, as noted in Fig. 1. It is because the microphysical effect, which is primarily initiated from the cloud base, is an indispensable part regardless of the presence or absence of the radiative effect that occurs when the microphysical effect becomes saturated (see "Discussion"). As aerosols near the cloud base can easily interact with the middle layers or even higher layers in thinner clouds than in thicker clouds, the microphysical effect can strongly modify cloud properties in thinner clouds than in thicker clouds. Since the coalescence growth process is weaker in thinner clouds, the strong microphysical effect that occurs before the saturation point is weakly opposed in thinner clouds, resulting the net effect dominated by the beforehand microphysical effect. On the other hand, stronger coalescence growth process can overcome weaker microphysical effect in thicker clouds, resulting the net effect dominated by coalescence growth process. This phenomenon plausibly describes the negative and positive relationships for relatively thinner and thicker clouds, respectively, over highly aerosol loaded regions, such as IGP. Overall speaking, prevalence of negative and positive relationships is the outcome of the competition between the microphysical and radiative effects.
Aerosol Index (AI), a product of AOT and Ångström exponent (AE), is also used as a proxy of aerosol number concentration (N tot ) while quantifying aerosol effects on cloud properties 5,13,14,28 , assuming that AI can represent N tot better than AOD 5 . AE is estimated either by paring AOTs of two wavelengths or through a least square fit of AOTs observed at discrete wavelengths. Therefore, depending on the choice of wavelengths, the calculated AE can differ 44,45 to affect AI as well, and thereby quantification of aerosol led impacts on cloud properties. Supplementary  Figure S3 confirms that quantified aerosol effects on cloud properties are subject to change depending on the choice of wavelengths for AE calculation. It is because AE is related to aerosol size distribution 46 ; and fine-mode and coarse-mode aerosols exhibiting larger sensitivities to the longer and shorter wavelengths, respectively, are differently weighted depending on wavelengths in AE calculation 47 . Figure S3 shows that such wavelength choice for AE (and AI) calculation can affect the quantified values for both thinner and thicker clouds; however, relatively, thinner clouds may be affected more prominently than thicker clouds. Among different wavelength combinations, the one with relatively wider spectral range (0.44-1.02 µm) shows results closer to Fig. 1a, which use AOT as a proxy of aerosol concentration. This suggests that AI estimated from AE of relatively wider spectral range may better represent aerosol loading, as such wider spectral range can more reasonably weight fine-and coarse-mode contributions in aerosol size distribution than narrow spectral range. Taking an opportunity of aerosol size distribution data (dV/dlnr) available in almucantar measurements 38 , we further tested if AI is indeed a good proxy for N tot . For this purpose, the calculated N tot values from aerosol size distributions (see "Methods") are compared with both AOT (0.5 µm) and AI separately, as shown in Supplementary Fig. S4. Figure S4 shows that the varying AE can largely deteriorate the relationship between AI and N tot , suggesting that AI may not be a better choice to consider a proxy of N tot over regions of varying AE (or aerosol type). On the other hand, the improved agreement for AOT and N tot relationship suggests www.nature.com/scientificreports/ AOT, rather than AI, better represents N tot over such regions. Therefore, AOT can better quantify aerosol led impacts on cloud properties over the study regions of this study. The aerosol led impacts on cloud properties are generally quantified in the form of aerosol indirect effect (AIE) as dlnCLD/dlnALD 13,14,21 , where CLD is CER or COT, and ALD is AOT or AI. We further followed this approach to quantify AIEs for each LWP bin. The quantified values corresponding to CER of 3.7 µm and COT of 0.645 µm are shown in Fig. 3a,b, respectively. Both Fig. 3a,b show qualitative agreement with Fig. 1a,b, respectively, by suggesting the negative and positive relationships for relatively thinner (low LWP) and thicker (high LWP) clouds, respectively. However, Fig. 3 does not show the increase in the strength of positive relationship with the increase of LWP very clearly, such as those observed in Fig. 1. This could be due to the influence of meteorological factors on cloud properties 6,28 , as such influences are left untouched in this approach. Despite it, Fig. 3 supports the robustness for negative and positive relationships for relatively thinner and thicker clouds, respectively, over our study regions, as shown in Fig. 1. At the same time, by highlighting the influences of meteorological factors on AIE, Fig. 3 suggests it important to disentangle these meteorological influences to better understand aerosol led impacts on cloud properties.

Discussion
There are two pathways for aerosols to influence cloud properties: microphysical and radiative 48,49 . In the former, aerosols directly interact with clouds by acting as CCN. But, in the latter, aerosols play an important role by absorbing 21,48 and/or scattering 50 the solar radiation and then altering the meteorological conditions, e.g., atmospheric stability and moisture content. The microphysical effect increases logarithmically with the increase of AOT until reaching a saturation point 48 . Though it depends on the meteorological condition, the microphysical effect generally saturates when AOT at 0.5 µm increases up to ~ 0.3 or more 21,48 . After saturation, the radiative effect gradually overrides the former microphysical effect 48 . Thus, the radiative pathway can play a critical role to modulate cloud properties over highly polluted regions with AOT (0.5 µm) exceeding ~ 0.3. Supplementary  Fig. S5 shows the frequency distributions of AOT (0.5 µm), corresponding to the direct sun measurement, and single scattering albedo (SSA) at 0.44 µm, corresponding to the almucantar measurement, for both study regions. The peak value of AOT (0.5 µm) falls above 0.3 over both study regions, suggesting the prevalence of high aerosol loadings over both study regions. Furthermore, the peak value of SSA (0.44 µm) is less than 0.95 over both study regions, suggesting that these highly loaded aerosols are quite light absorbing as well. Overall, aerosols of both study regions are highly capable to trap (absorb) solar radiation within the atmosphere with potential of modulating cloud properties through the radiative pathway mentioned above. Supplementary Fig. S6 shows the scatterplot for aerosol loading and the amount of solar energy trapped within the atmosphere, i.e., atmospheric forcing (ATM), which is the difference in aerosol radiative forcing between the top of the atmosphere (TOA) and the surface, over our study regions. These calculated ATMs and collocated meteorological data (see "Data") are used to illustrate how aerosol led atmospheric heating via solar radiation absorption can modulate the vertical distributions of atmospheric stability and moisture content. For this purpose, the calculated ATMs are sorted in an ascending order and binned into three segments of equal sample numbers. Then, the potential temperature (θ) gradient (-dθ/dp) and specific humidity (q) profiles of pressure (p) coordinated atmospheric layers corresponding to data of the first segment ("low atmospheric heating" case) and the third segment ("high atmospheric heating" case) of sorted ATMs are averaged. The mean values of ATM corresponding to the first and third segments are 26.9 ± 8.7 Wm −2 (28.4 ± 8.9 Wm −2 ) and 82.5 ± 20.2 Wm −2 (84.6 ± 20.3 Wm −2 ), respectively, for Kanpur (Gandhi College). Note that ATM can be converted into atmospheric heating rate as k' × ATM/ΔP, where k' is a constant term and ΔP is the pressure difference between the surface and TOA 51,52 . Along with the mean values of cloud top pressure (CTP), the vertical profiles of q (blue lines) and − dθ/dp (red lines) for "low atmospheric heating" case (solid lines) and "high atmospheric heating" case (dashed lines) are shown in Fig. 4. The positive value of − dθ/ dp suggests a stable atmospheric layer, and vice-versa. Figure 4 suggests that the increased atmospheric heating is associated with the decrease of moisture content within the atmosphere. Since both study regions have high  Figure 4 further suggests higher values of -dθ/dp at altitudes higher than cloud top heights in "high atmospheric heating" case compared to "low atmospheric heating" case. This suggests that atmospheric stability at higher altitudes can be enhanced by increased aerosols via radiative pathway. Such increased stability at such higher altitudes plays an important role to affect cloud properties near the cloud top, as revealed from Fig. 5. Figure 5 shows the correlations of LTS, the difference in θ between 700 and 1000 hPa, with CER (3.7 µm) and CER (1.6 µm). Over both study regions, LTS correlated strongly with CER (3.7 µm) than with CER (1.6 µm), indicating the important role of such increased stability on modulating the properties of clouds near the cloud top. Further, as stable air mass has a tendency to sink, enhanced atmospheric stability at such higher altitudes can amplify the entrainment mixing process near the cloud top 53 . The entrainment mixing process can also play an important role to decrease the cloud top height. As a result, Fig. 4 suggests the decrease of cloud top height in "high atmospheric heating" case compared to "low atmospheric heating" case. It is important to note that the entrained air causes both dilution and evaporation of cloud droplets, helping to broaden cloud droplet spectra [54][55][56] . As larger cloud droplets have higher terminal velocities, they fall faster and collide with smaller cloud droplets. The aggregated droplets fall even faster to collide and coalescence with smaller droplets in their path. This process becomes more efficient for wider distribution of cloud droplet size, as terminal velocity depends on cloud droplet size. Therefore, cloud droplet spectra broadened by entrainment mixing process generates a better condition for cloud droplet coalescence process. Thus, increased aerosols can contribute in cloud droplet coalescence process by increasing atmospheric stability at higher altitudes over both study regions. A modelling study has also shown the increase (decrease) of CER (cloud droplet concentration) through entrainment process in polluted clouds 53 to support the results of this study.

Data
AERONET. Level   . Vertical profiles of specific humidity (q) and potential temperature gradient (-dθ/dp) along with CTP for "low atmospheric heating" and "high atmospheric heating" cases for (a) Kanpur (KP) and (b) Gandhi College (GC) sites.  28 by further helping to accumulate more reasonable data of aerosols and clouds for aerosol-cloud interaction study, as aerosols beneath the cloud is difficult to observe during cloudy sky condition. The hourly metrological data of each pressure level were converted into time resolution of 1-min through cubic spline interpolation. Then, metrological data of ± 30-min centering the time ahead of cloud observation time by 3-h were averaged. We used air temperature (T), pressure (P), and specific humidity (q) data in this study. The potential temperature (θ) was calculated as where P 0 is pressure at 1000 hPa, R (= 287.05 Jkg −1 K −1 ) is the gas constant of air, and C p (= 1004 Jkg −1 K −1 ) is the specific heat capacity at constant pressure. The potential temperature gradient for any atmospheric layer i was calculated as The lower tropospheric stability (LTS) was calculated as Similarly, lifting condensation level (LCL), the height at which the relative humidity of air parcel becomes saturated, is calculated from vertical profiles of P and q using Metpy Python.
Calculation of aerosol radiative effect. In order to calculate the amount of solar flux trapped within the atmosphere due to aerosol absorption, we used AERONET observed spectral values of AOT, SSA, ASY, and surface reflectance at 0.44, 0.67, 0.87, and 1.02 µm wavelengths, and PWC. Using those data in a Santa Barba (3) LTS = θ 700hPa − θ 1000hPa . Multiple linear regression analysis. In order to quantify aerosol effects on cloud property (COT or CER) by decoupling the effect of meteorological parameters, a multiple linear regression method was applied by treating a cloud property as an independent variable and AOT and meteorological parameters as dependent variables. We formulated their relationship as where, and In these equations, CLD is a cloud property (COT or CER) and M i is any specific meteorological parameter (e.g., PWC). Further, bar and σ represent the mean and standard deviation values, respectively. In Eq. (10), a and b i are constant terms. We derive a by considering two meteorological factors, LST and PWC. The constant term a is ∂CER´/∂AOT´ (∂COT´/∂AOT´) when CLD is CER (COT).

Code availability
Data analyses codes are available upon request to the corresponding author.