Recent trend in the global distribution of aerosol direct radiative forcing from satellite measurements

Global distribution of aerosol direct radiative forcing (DRF) is estimated using Clouds and Earth's Radiant Energy System (CERES) synoptic (SYN) 1° datasets. During 2001–2017, a statistically significant change of global DRFs is revealed with a general decreasing trend (i.e., a reduced cooling effect) at the top of the atmosphere (DRFTOA ~ 0.017 W⋅m−2⋅year−1) and at the surface (DRFSFC ~ 0.033 W⋅m−2⋅year−1) with rapid change over the land compared to the global ocean. South Asia and Africa/Middle East regions depict significant increasing trend of atmospheric warming by 0.025 and 0.002 W·m−2⋅year−1 whereas, the rest of the regions show a decline. These regional variations significantly modulate the global mean DRF (−5.36 ± 0.04 W·m−2 at the TOA and − 9.64 ± 0.07 W·m−2 at the surface during the study period). The observed DRF trends are coincident with the change in the underlying aerosol properties, for example, aerosol optical depth, Ångström exponent and partly due to the increasing columnar burden of SO2 over some of the regions. This indicates that increasing industrialization and urbanization have caused prominent change in the DRF during recent decades.


| INTRODUCTION
Assessing global distribution of aerosol direct radiative forcing (DRF) and their long-term trends are vital to improve the state of understanding about the significance of aerosol climate forcing as well as the effectiveness of emission control policies (IPCC, 2013;Zhao et al., 2017;Yang et al., 2018;Aas et al., 2019). Study of aerosol radiation interactions are also essential for understanding the changes of the photolysis rates which influences other (viz., O 3 , NO 2 , SO 2 ) air pollutants (e.g., Li et al., 2014;Xing et al., 2015. In this regard, extensive measurements and radiative transfer simulations are used to estimate DRFs. However, the global estimation of DRFs has always been challenging especially due to inhomogeneity in both land surface and aerosol properties, posing most uncertain component for the estimation of cumulative radiative forcing (IPCC, 2013;Myhre et al., 2013;Murphy and Ravishankara, 2018). These uncertainties result from the uncertainty in (a) mixing state (including size and density) of aerosols (Seinfeld and Pandis, 2006;Kahn and Gaitley, 2015;Samset et al., 2018), (b) vertical profile of aerosols (Chung et al., 2005), (c) consideration of inaccurate aerosol species contributing to the total forcing in the models (Forster et al., 2007;Murphy, 2013;Myhre et al., 2013;Stevens, 2015;Paulot et al., 2018), their subsequent feedback to climate systems (Stocker et al., 2013); in addition to measurement errors and errors in models itself (IPCC, 2013). Lack of appropriate quantification of anthropogenic fraction of aerosol composites is also accountable to the large fraction of total uncertainty (Su et al., 2018;Myhre et al., 2013). For example, best estimated radiative forcing for the change in the net aerosol radiation interaction between 1,750 and 2005 is associated with a large uncertainty range of −0.85 to +0.15 WÁm −2 for the corresponding mean value of −0.35 WÁm −2 (IPCC, 2013).
The global average aerosol forcing due to both natural and anthropogenic aerosols was found to be −4.3 WÁm −2 at the top of the atmosphere (TOA), + 5.5 WÁm −2 in the atmosphere (ATM) and − 9.7 WÁm −2 at the surface (SFC) ( -2009( , Chung, 2012. Paulot et al. (2018) have reported an increase in the DRFs over Western Europe and eastern America by 0.7-1 and 0.9-1.4 WÁm −2 decade −1 and decrease over India by −1 to −1.6 WÁm −2 decade −1 during the period 2001-2015. This is in line with the fact that even though concentration of aerosols is seen to be decreasing in eastern America, central South America, Europe and northeast of Oceania, it is increasing in the rest of the world (Mao et al., 2014). Regionally, increase in aerosol burden and changes in their emission pattern is mostly attributed to anthropogenic activities (Dey and Girolamo, 2011;Babu et al., 2013;Moorthy et al., 2013;Kahn and Gaitley, 2015;Nair et al., 2016;Srivastava, 2017;Zhao et al., 2017). This is modulated further by the prevalent changes in meteorological parameters effectively enhancing the aerosol concentrations in the atmosphere (Yang et al., 2016). Thus, the change in DRFs is highly dependent on the aerosol types in addition to columnar abundance. For example, black carbon (BC) imparts larger positive effect (Bond et al., 2013;Zarzycki and Bond, 2010;Chung, 2012;Gogoi et al., 2017). Similarly, SO 4 and OC cause larger negative effect (Collins et al., 2002;Myhre et al., 2013;Paulot et al., 2018). Hence, trend analysis of DRFs itself is vital to estimate the change in radiative impact of aerosols. However, the global estimation of DRFs at different levels of the atmosphere that is, SFC, ATM and TOA is complex.
In the present study, the global distribution of DRFs and trends are uniquely studied using Clouds and Earth's Radiant Energy System (CERES; Wielicki et al., 1996) Synoptic (SYN) climate quality daily data sets. CER-ES_SYN1deg solar and longwave downwelling surface fluxes were evaluated against measurements at 85 land and ocean sites around the globe over the period from March 2000 to December 2007, revealing a monthly mean bias of 2.2 WÁm −2 (1.1%) for the shortwave and −4.1 WÁm −2 (−1.2%) for the longwave irradiances (Rutan et al., 2015). Moreover, CERES Science Team has applied the table of scaling factors for both Terra and Aqua to the observed TOA SW fluxes for the tuning of the computed fluxes. Thus, the synoptic data quality used in the present study is highly accurate and estimation of DRF over both regional and global scale at distinct atmospheric levels (i.e., SFC, ATM and TOA) adheres to improved estimation due to applications of updated CERES data sets. Based on this, the estimation of DRF and its trend in distinct geographical regions of the globe are themselves unique. The methodology adopted in the present study also serves as a primary proxy to accurately estimate global change in DRFs which is otherwise difficult to obtain from the limited in-situ measurements and highcomputation modeling. Another important aspect in the present study is the long-term data base (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017) chosen to examine the regional changes in aerosol loading in recent years over different parts of the globe. Thus, the estimation of DRF trends in the global scale in our study gives insight in to the identification of regional impact of aerosols on the Earth's energy budget associating distinct regional changes (increase/decrease).

| DATA AND METHODOLOGY
Daily global radiation data (CERES_SYN1deg-Day_Terra-Aqua-MODIS_Ed3A; Doelling et al., 2013) measured by CERES operating on-board Terra and Aqua satellite have been used to estimate the clear sky DRF at the TOA, ATM and SFC. The surface net radiations are estimated using the radiative transfer model with the atmospheric and surface properties obtained from MODIS, 3-hr geostationary (GEO) data and meteorological assimilation data from Goddard Earth Observation System (GEOS), while constraining both solar reflected and earth-emitted radiation for the TOA for all-sky and clear-sky conditions. Here, the radiative transfer calculation is made using the NASA Langley Research Center-Modified (LaRC) Fu-Liou radiative transfer model (Fu and Liou, 1993;Randles et al., 2013;Rose et al., 2013). The model is based on Fu-Liou scheme Liou, 1992, 1993). This code is a modified version of the original Fu-Liou code to improve treatment of Rayleigh scattering and updated shortwave (SW) and longwave (LW) gaseous absorption (Kato et al., 1999;Kratz and Rose, 1999). These computations use MODIS and geostationary cloud properties along with the atmospheric profiles of pressure, specific humidity and air temperature provided by NASA Global Modeling and Assimilation Office (GMAO; Reichle et al., 2009) reanalyzes as inputs. Furthermore, the absorption by water vapor, carbon dioxide, methane and oxygen are treated using the methodology followed by Kato et al. (1999). The computations are constrained by the observed CERES TOA fluxes (Rutan et al., 2015;Wang et al., 2017). The model is further constrained by the aerosol optical depth and scattering properties for the vertical column from the Model for Atmospheric Transport and Chemistry (MATCH; Gidhagen et al., 2005). The MATCH simulated aerosol properties are constrained by the observations from MODIS (Collins et al., 2001;Levy et al., 2013). The four modes of fluxes are used: pristine (clear, no-aerosols) and clear sky (clear, with aerosols) conditions at SFC and TOA. Here, 'clear' indicates 'cloud free aerosol skies'. The aerosol forcing at TOA (DRF TOA ) and surface (DRF SFC ) is estimated by differencing net clear sky from net pristine sky fluxes and in the atmosphere (DRF ATM ) by subtracting the two DRFs as follows: Using long-term values of DRF at TOA, ATM and SFC, trends have been estimated. The Mann-Kendall test (Kendall, 1975, Li et al., 2014Mann, 1945) has been used to identify significance of the changes at the 95% confidence level.
The MODIS on-board Terra and Aqua satellites provides aerosol retrievals over both land and ocean that are suitable for studying regional distribution of aerosol sources and trends. In this study, MODIS derived Aerosol optical depth (AOD, τ λ ), and Angstrom exponent (α) are used to understand aerosol loading and information on aerosol size distribution in the atmosphere. While AOD is the integrated extinction coefficient over columns of unit cross section measured at a particular wavelength (λ), values of given information on aerosol size distribution in the atmosphere. Following Angstrom (1929), α is estimated as: where, τ λ1 and τ λ2 are AOD at wavelength λ 1 and λ 2 , respectively.
α is inversely related with the particle size, that is, α ≤ 1 indicates size distribution dominated by coarse mode aerosols and α ≥ 1 by fine mode aerosols.
In addition, present study uses columnar SO 2 retrieved from OMI onboard Aqua satellite (Fioletov et al., 2013;Krotkov et al., 2016) for the period 2004-2017. OMI is the result of a partnership between NASA and the Dutch and Finnish meteorological institutes and space agencies flying on the NASA EOS Aura satellite Schoeberl et al., 2006). The OMI SO 2 product uses spectral measurements between 310.5 and 340.0 nm in the UV-2 (Li et al., 2013).  Figure S1). Collins et al. (2002) have reported that the surface insolation (under cloud free skies) over the Indian sub-continent is reduced by 40 WÁm −2 during spring (1999) due to the impact of natural and anthropogenic aerosols, accompanied by warming (25 WÁm −2 ) in the atmosphere. In the present study, the values of DRF ATM over the highly polluted Indo-Gangetic Plains (IGP) is as high as~30 WÁm −2 during 2017, along with a reduction in the surface reaching flux by about 45 WÁm −2 . East China shows higher surface dimming in the beginning of the study period which reduces during the later years. DRF TOA over Africa is comparatively low (~−12 WÁm −2 ) during both the years, but the DRF SFC and DRF ATM extend to −40 and 30 WÁm −2 , respectively. Rest of the land area experiences much lower values of DRF.
On the other hand, oceanic regions exhibit lower perturbation of the solar radiation due to lower aerosol loading, inferring an average value of 2.99 ± 0.04 WÁm −2 in the atmosphere associated with the cooling at the TOA by −5.24 ± 0.04 WÁm −2 and surface dimming of −8.23 ± 0.07 WÁm −2 (averaged from 2001 to 2017). Values of DRFs are higher near the west coast of Africa and east coast of Asia, similar to that reported by Zhang et al. (2005) where DRF (−5.3 WÁm −2 ) was estimated over cloud-free global oceans using the collocated MODIS and CERES data for the period November 2000-August 2001. Oceanic DRF TOA , DRF ATM and DRF SFC (−5.25, 2.94 and − 8.19 WÁm −2 , respectively) in 2001 was changed by 1.9, 1.4 and 1.7% (−5.15, 2.90 and − 8.05 WÁm −2 ), respectively in 2017. These changes are much lower compared to that of the land. Anthropogenic aerosols are found to contribute significantly on the total aerosol loading over the oceanic region predominantly over the regions where the values of DRFs are found higher (Boucher and Haywood, 2001;Bellouin et al., 2005;Li et al., 2014). Strong dust plumes originating from Asia reaches northern Pacific Ocean and west coast of North America. These regions exhibit annual DRF SFC of~− 14 WÁm −2 and DRF ATM of  Figure S1 ~8 WÁm −2 during 2001 (similar to Zhang et al., 2005). Fire activities are persistent near the west and east coast of South America resulting to atmospheric forcing of~12 WÁm −2 . The aerosol forcing has a strong seasonal and regional dependence which are addressed in various studies in the past (Hatzianastassiou et al., 2004;Quaas et al., 2008;Randles and Ramaswamy, 2008;Pathak et al., 2016;Nair et al., 2016;Paulot et al., 2018;Subba et al., 2018;Yang et al., 2018) and is not the primary focus of this study.

| Trends in global aerosol direct radiative forcing
The above results clearly indicate that in recent years, there has been potential change of the radiative balance of the atmosphere due to the impact of aerosols. Considering 2001 as base year, the DRF TOA , DRF ATM and DRF SFC trend over the globe are found to be 0.02, −0.02, and 0.03 WÁm −2 Áyear −1 , respectively ( Figure 2). The global land area exhibits rapid change of atmospheric forcing by −0.03 WÁm −2 Áyear −1 associated with an increase in DRF TOA and DRF SFC by 0.02 and 0.05 WÁm −2 Áyear −1 . However, oceanic region shows relatively weaker change in the DRF trend which is consistent with the comparatively lower change in the underlying oceanic surface and aerosol loading over majority of oceanic region. General trend in atmospheric forcing and surface dimming over oceanic regions is −0.01 and 0.03 WÁm −2 Áyear −1 , respectively, almost half of the values over land regions.
Over land regions, trend of DRF TOA are decreasing mostly over eastern and western of Africa, and south Asia, as against slightly increasing trend seen over North America and Europe (Figure 3). There cannot be one to one correlation between the radiative effect and the aerosols over every region. Especially over India, east and central China, South America, Africa due to the heterogeneity in the aerosol type (Paulot et al., 2018). Anthropogenic emissions over India and China are expected to be more uncertain than those over America and Europe (Pan et al., 2015;Saikawa et al., 2017). For example, India has been experiencing changes in surface albedo due to underlying surface properties (Paulot et al., 2018). Increase in rainfall and decrease in dust emissions (Pandey et al., 2017) have caused more greening over the north-western regions (Jin and Wang, 2018). The decrease in surface albedo over India masks the effect of the increasing anthropogenic emissions on the outgoing radiation (Paulot et al., 2018). This can change the signs of TOA forcing (Satheesh, 2002).
Since, aerosol forcing is dependent on the regional heterogeneity of land surface and aerosol loading, we have selected six distinct regions (R1-R6; R1-North-East America, R2-South America, R3-Europe, R4-Africa and Middle East countries, R5-south Asia and R6-eastern China) having significant trends and greater influence in overall DRFs. The trend of DRFs over each of these regions reveals substantial changes ( Figure 3, Table 1). These changes are much greater than that of the global averaged changes, indicating the heterogeneous spread of aerosols. Trend over R2 and R4 are least (DRF TOA < 0.012 WÁm −2 Áyear −1 ), while showing weak increasing trends over R1, R3 and R6 (0.074, 0.114 and 0.038 WÁm −2 Áyear −1 ). However, DRF trend over R5 is unique indicating decreasing trend of TOA cooling (−0.076 WÁm −2 Áyear −1 ) and increasing trend of ATM warming (0.025 WÁm −2 Áyear −1 ). These changes can be attributed to change in both mass loading and type of aerosols over these regions. Decreasing trend of cooling at TOA is attributed to increase in the fraction of absorbing AOD and decrease in aerosol single scattering albedo, respectively. Growth in the small sized particles, mostly of anthropogenic origin enhancing the cooling at the TOA. These properties are discussed in the following section.

| Regional variability of DRF
The long-term change in DRFs indicate that there has been profound change in aerosol burden across the globe in recent decades. We find potential enhancement of atmospheric warming (>0.04 WÁm −2 Áyear −1 ) over central, eastern IGP and southern India, southern and central part of Africa and Middle East countries and their adjacent water bodies. Rest of the globe shows the opposite. To understand the behavior of aerosols, we have investigated the spatio-temporal distributions of aerosol optical depth (AOD) and Ångström exponent (α) obtained from MODIS onboard Terra and Aqua satellites. Since the change in the DRFs in the recent decade is mostly attributed to the anthropogenic emissions, we have also studied the change in columnar SO 2 retrieved from OMI onboard Aqua satellite (Fioletov et al., 2013;Krotkov et al., 2016). Although SO 2 is just one of the precursor gases of the sulfate aerosol and not always the representative of the anthropogenic aerosols, we have examined its global distribution and trends to understand the widespread emissions by fossil fuel combustion, which is a major source of anthropogenic aerosols. We see that there have been significant positive trends of AOD over the continental regions, such as central and eastern parts of the IGP (R5), Middle East countries and the central part of Africa (R4), along with their adjacent water bodies such as the Arabian Sea, the northern Bayof-Bengal and the east coast of Africa (Figure 4a). Whereas, continental regions such as eastern parts of the United States, South America, European countries, eastern China and their adjacent water bodies show negative trends in AOD (−0.01 to −0.04 AOD Áyear −1 ). In line with this, the values of α (a qualitative indicator of the size of aerosol particles following the inverse relationship [Angstrom, 1929]; higher values indicate fine mode particles and vice versa) also exhibits significant changes over the period 2003-2017; especially over eastern China, India, Europe, the Middle East countries and Africa (Figure 4b).  This reveals that not only the amount of global aerosol loading has changed; simultaneously the physicochemical properties of aerosols have also been changed. However, change in total columnar AOD is not just caused evenly by the change in the amounts of all types of aerosols, but is rather driven by the predominant contributors of the composite aerosols. This can be observed in terms of the change of α, which especially increase over the Indian subcontinent and central China. This is attributed to increasing anthropogenic activities in response to growing population and urbanization resulting in higher emission of small sized particles. Various studies in the past have shown that the aerosol loading over the Indian subcontinent is highly influenced by the anthropogenic aerosol loading/activities over the region Moorthy et al., 2013;Srivastava et al., 2017;Pathak et al., 2016). SO 2 is one of the major precursors of the anthropogenic aerosols (Smith et al., 2011) and falls within the sub-micron range. Investigation of long-term SO 2 indicates that the increase in α coincides with the increase in SO 2 loadings (Figure 4c can be attributed to various factors (e.g., varying source processes and transport pathways, etc.). Over the Indian region, very good association between the increasing trends of AOD and SO 2 over the eastern IGP indicates the growing role of anthropogenic activities over this region, while the rest of country experiences the influence of the aerosol loading which are driven by various other sources. These include the emissions from thermal power plants, industrial, residential, transportation and commercial divisions, constructional areas, livestock farming and fertilizer applications, biomass burning, wild forest fires, desert dust plumes, marine aerosols and biogenic emissions, etc. (IPCC, 2007;IPCC, 2013;Dey and Girolamo, 2011;Streets et al., 2013;Benedetti et al., 2014;Matthias et al., 2018). Among these, emissions of dust, sulfate and organic carbon enhancing their abundance in F I G U R E 4 The trends of (a) aerosol optical depth (AOD) at 550 nm, (b) angstrom exponent and (c) columnar sulfur dioxide (SO 2 ). Dotted areas are the regions having significant trend above the 95% confidence level. The red box represents different regions R1-R6 the troposphere is significant (Fadnavis et al., 2019) in altering the incoming solar radiation. The western part of India, where dust is the primary contributor of aerosol loading, exhibits decreasing trend of DRF SFC driven by decreasing trend of dust aerosols in the last decade (Pandey et al., 2017). Similarly, contributions of carbonaceous components to the long-term trend of DRF over the Indian region cannot be ignored. From the study of BC trend over the Indian region, Manoj et al. (2019) have indicated the presence of carbonaceous aerosols in the elevated lower free-tropospheric regions, even though a general decreasing trend of BC near the surface was revealed from the study. In addition to distinct source processes, influence of regional climate variables (e.g., seasonal and regional rainfall) on the variability of aerosol loading is important. Babu et al. (2013) have shown that increasing trend in AOD over the IGP can be attributed to the competing effects of aerosol influx and reduced wet removal impact of aerosols. Thus, the sharp increase in atmospheric warming over the Indian sub-continent can be attributed to the complex interaction between natural and anthropogenic forcing at different spatio-temporal scales, modulated by regional climate variables.
Based on the above observations, increasing trend in TOA cooling can be synchronous with the change in anthropogenic aerosol loading, where the emission of SO 2 hold significance over eastern America, parts of Africa and eastern China and the adjoining oceanic region. Zhao et al. (2017) attributed the decreasing trend of DRFs over eastern United States and Europe to decrease in SO 2 and NO x which are the precursors for the secondary aerosol formations, except for mineral dust and ammonia especially after 2011. Similarly, 25-50% of organic carbon (30% of PM) over US was found to decline between 1990 and 2012, with no significant change in biogenic aerosols (Ridley 2018;Blanchard et al., 2013Blanchard et al., , 2016Attwood et al., 2014). This is attributed to decline in the anthropogenic emissions predominantly from vehicular and residential fuel burning related to series of control measurements under the Clean Air Act (Xing et al., 2013).
Interestingly, R1 exhibits increasing trend until 2007 (DRF ATM trend: 0.13 WÁm −2 Áyear −1 ) and a decreasing trend after 2007 (DRF ATM trend: −0.15 WÁm −2 Áyear −1 ). Similar bimodal trend was observed by Paulot et al. (2018) over the same region. Decrease in DRFs is indicative of the positive result of the implementation of various air pollution acts (Wang et al., 2014). However, it is important to note that though there has been decrease in DRFs over East China (Figure 2), their magnitudes in the recent years are still higher over some of the locations compared to the rest of the world (Figure 1 and Figure 2). Hence, continuous mitigation of air quality is still important.

| CONCLUSIONS
Present study demonstrates an extensive utilization of the combination of space-borne radiation flux measurements and radiative transfer modeling for improved quantitative estimation of global aerosol DRF trend over the period from 2001 to 2017. Since there exist noticeable discrepancies in the estimation of aerosol trends from various satellite retrievals (e.g., MODIS, MISR, OMI), most likely due to the difference in spatio-temporal sampling, observing strategy and retrieval algorithms (Kahn et al., 2009;Zhao et al., 2017), consideration of direct radiation measurements by satellite sensor is unique and suitable for global radiation budget studies. Thus, utilization of highly accurate synoptic (SYN) datasets from Clouds and Earth's Radiant Energy System (CERES) provides updated information on the global distribution of DRF trend, which is otherwise difficult to obtain from the limited in-situ measurements and high-computation modeling. The major outcomes of this study include: 1. In the recent decades, there have been substantial changes in the regional and global distribution of DRFs due to aerosols. The global (−60 S to 60 N) values indicate a decreasing trend of TOA forcing (i.e., a reduced cooling effect) by 0.017 WÁm −2 Áyear −1 , along with a decreasing trend of atmospheric forcing (i.e., a reduced warming effect) by −0.017 WÁm −2 Áyear −1 having rapid change over the land (− 0.029 WÁm −2 Áyear −1 ) compared to that over the ocean (−0.012 WÁm −2 Áyear −1 ). While many distinct geographic regions (over land) of the globe depicted a general decreasing trend of the atmospheric forcing (warming) ranging from −0.053 WÁm −2 Áyear −1 in North-east America to −0.076 WÁm −2 Áyear −1 in Western Europe, Africa and South Asia showed increasing trend (i.e., an enhanced warming effect in the atmosphere) by 0.002 WÁm −2 Áyear −1 and 0.025 WÁm −2 Áyear −1 respectively.
2. Global map of AOD reveal significant positive trend over the continental regions, such as central and east Indo-Gangetic basins of India, Middle East countries and Africa, along with their adjacent water bodies such as Arabian Sea, North Bay of Bengal and east coast of Africa.
3. The global trend of DRF can be attributed largely to the anthropogenic aerosols. The increasing trend of SO 2 over eastern IGP, a small region of Africa, Middle East, and South America indicates increase in anthropogenic activities over these regions and may be a dominant factor in governing the DRF trend over these regions.

ACKNOWLEDGEMENT
This work was carried out as part of the Aerosol Radiative Forcing over India (ARFI) project of ISRO-GBP. The authors thank CERES, OMI and MODIS science team for providing scientific datasets, which have been obtained through the NASA Langley Distributed Active Archive Systems and GIOVANNI portal. Tamanna Subba was supported by ARFI project fellowship from ISRO-GBP and Kalam-Climate Doctoral Fulbright Fellowship program 2017-2018 under USIEF. Binita Pathak is a junior associate of ICTP, Italy.