Long-term spatiotemporal trends in aerosol optical depth and its relationship with enhanced vegetation index and meteorological parameters over South Asia

Satellite-based aerosol optical depth (AOD) is columnar light extinction by aerosol absorption and scattering and has become the most important variable for the assessment of the spatiotemporal distribution of aerosols at a regional and global level. In this paper, we have used AOD observations of multiangle imaging spectroradiometer (MISR) from September 2002 to May 2017, moderate resolution imaging spectroradiometer (MODIS) from September 2002 to December 2020, and sea-viewing wide field-of-view sensor (SeaWiFS) from September 2002 to December 2010 over South Asia. We have observed the association of AOD with enhanced vegetation index (EVI) and meteorological variables (temperature (TEMP), wind speed (WS), and relative humidity (RH)) acquired from Giovanni during the period September 2002–December 2020. The satellite observations of Terra-, MISR-, and SeaWiFS-AOD were also compared with Aqua-AOD. The findings show that AOD in eastern Pakistan is higher than in the western Pakistan due to increase in population density and biomass burning. Mean annual peak AOD (˃ 0.7) has been observed over the IGB region because of the significant increase in economical, industrial, and agricultural activities while AOD of ˃ 0.6 is observed over Bangladesh. The lowest mean annual AOD of ˂ 0.3 is observed over northeastern Afghanistan, western Nepal, and Bhutan whereas the AOD of 0.3 is seen over Sri Lanka. The highest seasonal mean AOD of 0.8 has been seen over Bihar, India, and AOD of ~ 0.7 is observed over Bangladesh while the lowest AOD is observed over Afghanistan, Sri Lanka, Nepal, and Bhutan during the winter season. However, the mean AOD over eastern Pakistan is maximum in both monsoon and post-monsoon season but relatively low in pre-monsoon and winter. The highest positive seasonal AOD anomalies were observed over South Asia in winter season followed by post-monsoon, pre-monsoon, and least being monsoon. The higher mean AOD anomaly value is found to be 0.2 over eastern Pakistan and western India. In northeastern Pakistan and central India, AOD and RH are positively correlated (r ˃ 0.54) while negatively correlated over Afghanistan, southwestern region of Pakistan, eastern India, Nepal, Bhutan, and Bangladesh. AOD is negatively correlated (r =  ~  − 0.3) with EVI over eastern Pakistan and western India. The highest correlation coefficient (r) obtained among Terra and Aqua is 0.97, MISR and Aqua is 0.93, and SeaWiFS and Aqua is 0.58 over South Asia.


Introduction
Atmospheric aerosols are solid or liquid particles ranging from 0.001 to 100 µm and suspended in the earth's atmosphere (Bilal et al. 2019;Ling and Han 2019). These particles have remarkable impacts on the air quality (ul-Haq et al. 2017;Jung et al. 2019), visibility (Han et al. 2013), climate system (Tariq and Ul-Haq 2018;Tariq 2020), and human health (Ren-Jian et al. 2015;Mehmood et al. 2021;Qayyum et al. 2021) through radiation forcing. Aerosols can also warm or cool down the surface depending on their size, source, and location. The emission of aerosols is from both natural and anthropogenic sources. Anthropogenic emissions of aerosols include biomass and crop residue burning, fossil fuels burning, transportation, and industries while natural emissions are desert dust, volcanic eruptions, sea spray, and forest fires. Therefore, an accurate understanding of aerosols variability, types, and sources are needed for air quality and climate modeling.
Non-uniformity in their sources and a short life span of few days in the atmosphere limit the understanding of the spatiotemporal distribution of aerosols (Logan et al. 2013). The spatiotemporal variability of aerosols is examined using continuous observations of aerosol optical depth (AOD). AOD is defined as columnar absorption and scattering of aerosol content in the atmosphere. For a better understanding of aerosol distribution in the atmosphere at high spectral and temporal resolutions, an in situ observational instrument aerosol robotic network (AER-ONET) was used globally for continuous observations of aerosol properties (Holben et al. 1998). However, these in situ observations provide accurate information on aerosols but are spatially limited (Tariq and Ali 2015). To better address this spatial limitation problem, satellite datasets have been used to monitor the spatiotemporal distribution of aerosols at territorial and global scale (Remer et al. 2005;Nichol and Bilal 2016;Tariq and Ul-Haq 2018). Currently, various satellites in space have instruments for the accurate estimation of AOD such as moderate resolution imaging spectroradiometer (MODIS) onboard Terra and Aqua platforms (Kaufman et al. 1997), advanced very high-resolution radiometer (AVHRR) (Riffler et al. 2010), ozone monitoring instrument (OMI) (Kumar et al. 2015), seaviewing wide of view sensor (SeaWiFS) (Sayer et al. 2012), and multiangle imaging spectroradiometer (MISR) . Therefore, utilizing different aerosol products from several satellite sensors can overcome the limitation of each sensor. Moreover, there has been significant improvement in estimation of AOD using single spectral band (Bilal et al. 2013).
In recent years, aerosol emissions from both human and natural activities have reduced the air quality of Asian countries (Bilal et al. 2021;Mhawish et al. 2020;Mhawish et al., 2021a, b;Singh et al. 2017). Pakistan, India, Nepal, Sri Lanka, Afghanistan, Bangladesh, Maldives, and Bhutan are the countries included in the South Asia and consist of ~ 40% population of the Asian region. Particulates of complex mixtures of light scattering and absorbing are found in Asian aerosols. Brown carbon, iron oxide, and soot or black carbon are three main light-absorbing species while vehicular emissions and industrial sources are the main light-scattering aerosols.
Several recent studies have been carried out on spatiotemporal and seasonal variations of aerosol characteristics using both satellite and ground-based observation over South Asia (Alam et al. 2016;ul-Haq et al. 2017;Kumar et al. 2018;Ali et al. 2020). Prasad and Singh (2007a, b) found strong seasonal variability of AOD over the Indo-Gangetic plains using MODIS and ground-based aerosol robotic network (AERONET) data during the period 2000-2004. Tariq and Ul-Haq (2018 analyzed aerosol properties over Karachi during the period September 2006-August 2014 and reported peak AOD in July due to higher concentration of coarse mode particle. Kumar et al. (2013) identified the homogeneous aerosols region in south Indian regions using MODIS-Terra-AOD from 2001 to 2008. Bilal et al. (2021) reported the higher concentrations of PM 2.5 over Pakistan using Copernicus Atmosphere Monitoring Service data with mean value of 54.7 μg/m 3 from 2003 to 2020. Mhawish et al. (2021a, b) found previously unidentified aerosols hotspots over eastern Indian coastal state of Odisha, West Bengal, and Bihar due to the dominance of small particles of aerosols using high spatial resolution data from MISR and MAIAC from 2000 to 2019 over South Asia. All the previous studies conducted so far over South Asia have limitations in terms of datasets range and analysis techniques.
The rapid increase in urbanization, industrialization, population density, and anthropogenic activities has been experienced in mega-cities of South Asia in the recent years. Spatiotemporal distribution of aerosols in South Asia is poorly understood because of limited monitoring networks and satellite observation at high resolution (Guo et al. 2011). The study on variability of aerosols and their association with enhanced vegetation index (EVI) and meteorological variables over Pakistan has been conducted by Tariq et al. (2021) using the long-term data from September 2002 to August 2015. The authors have used limited long-term datasets, different analysis techniques, and considered only the Pakistan region . The main aim of present study is to observe the spatiotemporal distributions of aerosols over South Asia using datasets of 18 years. For this purpose, we have analyzed the association of Aqua-AOD with EVI and meteorological variables during the period September 2002-December 2020 over South Asia. Furthermore, the seasonal anomalies of Aqua-AOD were also calculated to examine the change in aerosol loadings over South Asia from 2002 to 2020. Additionally, we have used Terra-MODIS, SeaWiFS, and MISR-AOD data for the validation of Aqua-MODIS over South Asia using the reduced major axis regression technique.

Study area
The rapid increase in population and economy of South Asia in recent years has caused a significant increase in energy utilization, which has drastically increased AOD over South Asian countries. Afghanistan, Nepal, Maldives, Bhutan, Sri Lanka, India, Bangladesh, and Pakistan are the main countries of South Asia as shown in Fig. 1. Geographically, it consists of valleys, rain forests, deserts, glaciers, and grasslands. The climate of South Asia varies from area to area, from temperate in the north to tropical monsoon in the south. The northern part of the Indo-Gangetic belt (IGB) has high temperatures in the summer season and relatively low in winter. It can be divided into five main regions: the Karakoram and Himalayan mountains in the north, the southern lowlands, the Baluchistan Plateau, peninsular India, and the island realm.

MODIS
MODIS sensor aboard the Aqua and Terra satellites has provided daily global measurement of aerosols with a swath width of ~ 2330 km and creating geophysical datasets. Aqua and Terra are polar-orbiting satellites at an altitude nearly equal to 700 km. MODIS sensor has 36 spectral bands (0.41 to 14 μm) at variable spatial resolutions (250 m, 500 m, and 1 km) and temporal resolution of 1-2 days. In this study, we have used deep blue (DB) Aqua-MODIS (MYD08_M3 v6.1) and Terra-MODIS (MOD08_M3 v6.1) monthly product at 1° × 1° spatial resolution and monthly mean EVI product (MYD13C2 v006) during the period September 2002-December 2020 acquired from Giovanni (http:// giova nni. gsfc. nasa. gov).

MISR
MISR onboard Terra satellite was launched by NASA (National Aeronautics Space and Administration) in 1999 at an elevation of 705 km and has a temporal resolution of 16 days as well as the spatial resolution of 250 m, 275 m, and 1 km. It is intended to monitor all patterns and kinds of vaporized particles created from anthropogenic and natural activities. We have used MISR-AOD (MIL3MAE v4) data product at 0.5° × 0.5° spatial resolution during the period September 2002-May 2017 acquired from Giovanni (http:// giova nni. gsfc. nasa. gov). Kahn et al. (2010) reported detailed data on MISR instrumentation, recovery algorithms, and methodology.
SeaWiFS SeaWiFS onboard Sea Star satellite was first introduced in 1997 and has a temporal resolution of 2 days with a swath width of 1,502 km. It can measure the small contribution of the ocean to the reflected solar radiance at top of the atmosphere with great accuracy and stable calibration (Jr. et al. 2011). Hsu et al. (2006) established DB aerosol retrieval algorithm for SeaWiFS over land. SeaWiFS makes its AOD more reliable because of its high calibration and long-term data range. Monthly AOD (SWDB_L3M10 v004) product at 1° × 1° spatial resolution available from September 2002 to December 2010 over South Asia was acquired from http:// giova nni. gsfc. nasa. gov.

Methodology
To investigate the spatiotemporal distribution of AOD over South Asia and to find its relationship with meteorological parameters, the following methodology was used in this study:

MODIS-Aqua retrieved Deep Blue Aerosol Optical
Depth 550 Land (MYD08_M3 v6.1) product was used to make the mean spatial and seasonal maps of AOD and to investigate the long-term variations of mean annual spatial distribution of AOD over South Asia during the period September 2002-December 2020. The mean annual spatial AOD map is created by averaging the Aqua-AOD data over South Asia from September 2002 to December 2020. Moreover, seasonal AOD anomalies over study region were calculated using Eq. 1. The selected spatial resolution of Aqua-, Terra-, and SeaWiFS-AOD is 1° × 1° while MISR-AOD is 0.5° × 0.5°. The errors and accuracy in the datasets used for the validation are calculated using relative mean bias (RMB), expected error (EE), correlation coefficient (r), and root mean square error (RMSE) (Bilal et al. 2017). The reduced major axis regression technique is used to calculate the slope ( ) and intercept ( ) between both the variables (Bilal et al. 2021).
The independent parameter used in this study is Aqua-AOD and the dependent parameters are SeaWiFS-AOD, MISR-AOD, and Terra-AOD. The lowest value of RMB, lowest value of RMSE, highest value of r, and highest percentage within the EE retrievals are used to evaluate the performance of Aqua-AOD. The relative percentage difference in number of collocated observations is calculated using Eq. 7.
The DependentparameterAOD and IndependentparameterAOD) represents the mean values of AOD. The RMB value ˂ 0 represents the underestimation while RMB greater than ˃ 0 represents the over estimation of AOD retrieval and RMB value of 0 represents no under and over estimations of AOD. The β represents the slope, α represents the intercept of reduced major axis regression, and the σ represents the standard deviation. The higher accuracy between two parameters can be achieved when the RMSE value is closed to the zero.
(2) EE = ±(0.05 + 0.20 × Parameter AOD) 3. The spatial correlation maps were made to better understand the correlation of Aqua-AOD with meteorological parameters and EVI over South Asia from 2002 to 2020. For the spatial correlation maps, both the parameters are required to pair in space and time. If the spatial resolutions of both the parameters are different then the fine resolution dataset is regridded to the coarse resolution using the application named lats4d. 4. The monthly time series of Aqua-AOD was made to asses the temporal variations and trend of AOD over South Asia from 2002 to 2020. Furthermore, monthly and annual time series of MODIS-Aqua-AOD were plotted over selected cities of South Asia to examine the long-term trends of aerosols. The annual percentage change in AOD over selected cities was calculated using Eq. 7. Figure 2 exhibits the seasonal mean of meteorological variables (specific humidity (SH), surface air temperature (SAT), and WS) over the study area during the time span of 2002 to 2020 using NCEP/NCAR reanalysis data. It can be seen from Fig. 2a that the Indian IGB region and Bangladesh has SAT of 20 C, SH of 6 g/kg, and WS of 2 to 4 m/s while northeastern and southwestern Pakistan has 15 C, SH of 4 g/ kg, and WS of 2 m/s during the winter season.

Variations of meteorological parameters
The lowest SAT is observed over Afghanistan with SH of 4 g/kg and WS of 2 m/s during December to January. The SAT over southwestern Afghanistan is 25 C, northeastern Pakistan, IGB, Sri Lanka, and Bangladesh is 30 C during March to May. Figure 2b exhibits that the eastern and southeastern region of Afghanistan, eastern part of Pakistan, and central India have SH of 6 g/kg whereas the Bangladesh and Sri Lanka have SH of 10 g/kg during the period March to May. The eastern Nepal and Bhutan have SH of 15 g/kg with WS of 4 m/s during March to May. The WS over Afghanistan, eastern Pakistan, and Sri Lanka is 2 m/s while over India and Bangladesh is 4 m/s as shown in Fig. 2c. During June to August, the SAT over southern Afghanistan, eastern Pakistan, and IGB region is found to be 35 C while over Bangladesh and Sri Lanka is 30 C. The SH of 14 g/kg is found over Bangladesh and India whereas the SH of 12 g/kg is observed over Pakistan and Sri Lanka. Highest SH of ˃ 18 g/kg is observed over Nepal and Bhutan during the period June-August. In Fig. 2c, WS of In September 2002-November 2020, the SAT of 30 C is observed over eastern Pakistan, India Bangladesh, and Sri Lanka while 25 C is seen over southwestern Pakistan and southeastern Afghanistan. The SH over central Afghanistan, eastern Pakistan, and western India is 8 g/kg while over Bangladesh is 10 g/kg and Sri Lanka is 12 g/kg during the period September 2002-November 2020. The WS of 4 m/s is observed over southeastern Pakistan, IGB, Nepal, Bhutan, and Bangladesh whereas the WS of 2 m/s is seen over Afghanistan, northern Pakistan, and Sri Lanka as shown in Fig. 2c.  obtained between Aqua-and Terra-AOD due to the same aerosol retrieval algorithm for both sensors while RMB of 0.43% represents the underestimation of Terra-AOD as compared to Aqua-AOD as shown in Fig. 3a.

Validation of Aqua-AOD with Terra, MISR, and SeaWiFS retrievals
The value of the y-intercept is 0.03, slope is 0.88, and RMSE is 0.02. The EE value indicates that the 13.05% retrievals of Aqua-and Terra-AOD are within the expected error. The r of 0.93 is obtained among MISR and Aqua-AOD with slope of 0.95 and y-intercept of 0.003. The RMB of − 3.78% shows the underestimation of MISR with 11.23% of the data in EE. The RMSE value obtained between MISR and Aqua-AOD is 0.03. Martonchik et al. (1998) reported that the different aerosol retrieval algorithm and wavelength of MISR and Aqua-MODIS cause this underestimation. The lowest r of 0.58 is obtained among SeaWiFS and Aqua-AOD over South Asia while 11.57% of the AOD retrievals are within the EE. The RMB is − 10.68% which represents the underestimation of SeaWiFS-AOD over South Asia. This underestimation of AOD is due to the calibration of instruments (Remer et al. 2005). The RMSE value of SeaWiFS and Aqua-AOD over South Asia is 0.07 with slope of 0.87 and y-intercept of 0.06. However, Terra performs better with Aqua as compared to MISR and SeaWiFS in terms of RMSE, the percentage of AOD retrievals within the EE, and correlation. These findings depict the better performance of the Terra-AOD product as compared to MISR and SeaWiFS over South Asia. The % relative difference of EE obtained between Terra, MISR, and SeaWiFS with Aqua is − 13.25%, − 44.60%, and − 43.27% respectively over South Asia. Figure 4 represents the distributions of annual mean Aqua-AOD over South Asia during the period September 2002-December 2020. The mean annual AOD values (0-1) retrieved from Aqua-MODIS show significant spatial variability over South Asia. The findings from Fig. 4 show that AOD had a marked impact on countries in South Asia, namely Afghanistan, Nepal, Maldives, Bhutan, Sri Lanka, India, Bangladesh, and Pakistan. High AOD (~ 0.7) is observed over eastern Pakistan while AOD ˃ 0.6 exists over southeastern Pakistan. This includes an industrial area of Lahore, Gujranwala, Sialkot, Faisalabad, Jacobabad, Larkana, Hyderabad, and Nawabshah. Higher levels of AOD over eastern region of Pakistan are not only due to the local emissions but the transportation of aerosols due to crop residue burning from Indian region are also responsible for low visibility and destruction of air quality over Pakistan (Tariq and Ali 2016). The aerosol emissions from burning of biomass, fuel combustion, and crop residue burning in Indian region are transported to northeastern region of Pakistan and are also supported by wind pattern (Tariq and Ul-Haq 2018). Sharif et al. (2015) also reported high AOD over the southeastern region of Pakistan using MODIS-AOD from 2001 to 2013. The highest mean annual AOD (˃ 0.7) has been observed over the IGB region. Recent studies have reported that high AOD over the IGB region is due to the increased industrialization, burning of crop residue, high population density, burning of biomass, heavy urbanization, burning of coal, and meteorological conditions (Jethva et al. 2005;ul-Haq et al. 2017;Ojha et al. 2020;Bilal et al. 2021). The annual mean AOD over Bangladesh is ˃ 0.6 due to the emissions from brick kilns, transportation of sea salts from Arabian Sea, and desert dust particles from Indian region (Mainul et al., 2015;Mamun 2014). In addition, low population areas in this study region also contribute to aerosol emissions from garbage and wood burning. The low AOD values over northwestern region of Pakistan are due to less population and reduced agricultural activity. Low AOD values (˂ 0.3) are observed over northeastern Afghanistan, western Nepal, and Bhutan while AOD of 0.3 is observed over Sri Lanka. Mean annual AOD values of ˃ 0.3 are observed consistently over central parts of India is due to the local emissions and long-range transportation of aerosols from northwestern region of IGB. Similar to the outcomes acquired by Nizar and Dodamani (2019), a significant decline in aerosol loading is seen over the Himalayas and Tibetan Plateau as compared to the IGB. Figure 5 represents the seasonal distribution of averaged Aqua-AOD at a wavelength of 550 nm over South Asia from 2002 to 2020. The seasonality of AOD, despite some variations observed from area to area, in general, demonstrates lower AOD observed during pre-monsoon and winter season while higher AOD was seen during the period of post-monsoon and monsoon, depending upon the area. During the winter season, the highest AOD (0.8) is observed over the Bihar state of India as shown in Fig. 5a. The AOD of ˃ 0.6 is seen over the IGB region. The emissions from anthropogenic aerosols in IGB regions are responsible for high AOD over central IGB and are transported by northwestern winds over the whole IGB (Mhawish et al. 2020). Moreover, IGB region is often surrounded by haze activity during the winter season which resulted in heavy aerosol loading over the region (Badrinath et al., 2009). Similarly, the 0.5 AOD is seen in eastern and southern region of Pakistan. Bilal et al. (2021) also find low values of AOD during winter season as compared to the summer season over Pakistan during the period 2003-2017. It is evident from Fig. 5a that the areas of northwestern and southwestern Pakistan show the minimum AOD of ˃ 0.2. The mean highest AOD (~ 0.7) is observed in winter over Bangladesh. The lowest AOD values are observed over Afghanistan, Sri Lanka, Nepal, and Bhutan as shown in Fig. 5a. In the pre-monsoon, the highest AOD value of ˃ 0.6 is seen over the IGB region and Bangladesh as shown in Fig. 5b. Mhawish et al. (2021a, b) reported that in premonsoon season, IGB is surrounded by dust aerosols from desert regions mainly over central and northwestern Indian region. Recent studies have reported the transportations of mineral dust aerosols from Thar and Arabian desert lead to high AOD over IGB region (Dey et al. 2004;Sijikumar et al. 2016). The 0.5 AOD is observed in the eastern region of Pakistan and ˃ 0.3 is seen over central India. Higher values of AOD over eastern region of Pakistan and western region of India as shown in Fig. 5b are because of the presence of dust particles and crop residue in the atmosphere and these aerosol particles are transported to far distances by higher wind speed (Ali et al. 2014). Tariq et al. (2021) reported that the dust aerosols over Pakistan in pre-monsoon season are originating from south and southwestern region. The lowest AOD of 0.2 is observed in pre-monsoon over northwestern Pakistan, Bhutan, and northern parts of Afghanistan. Low AOD values over northwestern Pakistan are associated with wide agricultural area and less population density (Tariq and Ul-Haq 2018;Tariq et al. 2021).

Spatial distribution of AOD
In the monsoon season, high AOD ranging from ˃ 0.6 to 0.9 has been observed over the IGB region, Bangladesh, and eastern Pakistan. The high AOD during monsoon over IGB and eastern region of Pakistan is due to dust storms activity (Bilal et al., 2021;Mhawish et al., 2021a, b). Dust aerosols from Arabian desert due to the strong westerly winds are responsible for enhanced values of AOD over northwestern region of India (David et al. 2018). The 0.8 AOD value is also observed in the southeastern region of Pakistan. Tariq and Ali (2015) reported that higher concentrations of AOD over southern region of Pakistan are due to the high wind speed and temperatures. The high levels of AOD over Pakistan are also due to the long-range transport of desert dust aerosols from African and Arabian regions (Prasad and Singh 2007a). Similarly, ˃ 0.6 AOD value is noticed in the eastern region of India, southern part of Afghanistan. It is evident from Fig. 5c that the minimum AOD of ˃ 0.3 is seen in northwestern Pakistan, Sri Lanka, and the southern region of India. The high AOD levels observed during monsoon are due to the anthropogenic emissions and long-range transport of aerosols by the southwestern winds, even though rainfall is expected to play remarkable role in decrease of AOD values over IGB. In the post-monsoon season, ˃ 0.7 AOD value is detected over few areas of the IGB region as shown in Fig. 5d. The high AOD over IGB region is due to crop residue burning in Pakistan and Indian regions (Mhawish et al. 2020). Approximately 0.5 AOD value is noticed in central India and Bangladesh. The dust transport to Indian region is reduced during post-monsoon season due to the weak strength of the southwesterly winds. The northwesterly winds transport the aerosols over the central and eastern regions of IGB during the post-monsoon season. Similarly, ˃ 0.3 AOD value is observed in the southwestern region of Pakistan. Finally, it is noticed from Fig. 5d that the northern region of Pakistan, eastern area of Afghanistan, western part of Nepal, and Bhutan show a minimum AOD of 0.2. Chung et al. (2005) analyzed that anthropogenic aerosols contribute about 70% to the overall AOD over South Asia.
Seasonal distributions of AOD anomalies over South Asia from 2002 to 2020 are shown in Fig. 6. During the winter season, a high anomaly value of 0.2 is observed over few areas of eastern Pakistan and western India as shown in Fig. 6a. Major sources of high AOD anomaly values over these regions are due to excessive crop residue burning, urbanization, and industrialization. The lowest anomaly values were observed over eastern IGB and Bangladesh. The anomaly value of 0.1 is seen over Sri Lanka, Nepal, southern Afghanistan, and Pakistan in winter season. In all the seasons, regions with low AOD anomaly values were situated in less-developed and rural regions of southwestern and northeastern Pakistan. In the pre-monsoon season, the anomaly of 0.1 is seen over the IGB region and eastern Pakistan while negative anomaly values are observed over Bhutan and Nepal. In Afghanistan and Bangladesh, the anomaly value of less than 0.1 is observed during the premonsoon. The most noticeable feature is that (only during  Fig. 6c. In terms of the magnitude of AOD anomalies and area coverage, it is worth presenting that positive AOD anomaly centers were observed during winter, followed by post-monsoon, pre-monsoon, and least being monsoon. In post-monsoon, the positive AOD anomaly is observed over Afghanistan, southeastern Pakistan, western areas of India, Nepal, Bangladesh, and Sri Lanka while lowest values are observed in central India as shown in Fig. 6d. Figure 7 represents spatial correlation maps of AOD with meteorological variables (TEMP, RH, and WS) and MODIS-Aqua retrieved (EVI) over South Asia during the period September 2002-December 2020. A correlation between AOD and TEMP of ˃ 0.26 is observed over the eastern region of Pakistan and the western region of India whereas the highest correlation ˃ 0.40 is seen over the eastern Afghanistan, southwestern region of Pakistan, and western Nepal as shown in Fig. 7a. Deng et al. (2012) reported that meteorological variables such as pressure, WS, RH, TEMP, and wind direction change the chemical composition and path of aerosols. The lowest correlation between AOD and temperature is observed over Bangladesh and few parts of IGB. The highest positive correlation of ˃ 0.54 is obtained among AOD and RH over the northeastern region of Pakistan and the central part of India while a negative correlation is observed over Afghanistan, southwestern region of Pakistan, eastern India, Nepal, Bhutan, and Bangladesh as shown in Fig. 7b. The positive correlation of RH will lead to the increase in hygroscopic growth of aerosols and thus reduce the visibility . The decline in RH correlation on the western side of the Indian region accompanied by increase in RH over eastern regions of India. The highest correlation among AOD and WS is observed over the southeastern region of Pakistan, central parts of Afghanistan, and the western region of India while the negative correlation is observed over southern India, IGB, and Bangladesh. Tiwari et al. (2018) reported that the high WS carries the aerosols over Pakistan from Middle East, Afghanistan, Sahara, and Thar deserts. The maximum negative correlation of − 0.58 is observed among the AOD and EVI over the western region of India and southern region of Pakistan as shown in Fig. 7d. EVI can provide evidence about the changes in vegetation content. Chi et al. (2019) reported that the long-range transport of aerosols is restricted to the areas with high vegetation cover because vegetation provides resistance to the wind. The negative correlation shows that an increase (decrease) in AOD results in a decrease (increase) in EVI over the area. In eastern Pakistan and western India, a negative r (~ − 0.3) among the AOD and EVI is seen due to the presence of high aerosols loading and low EVI. Tariq et al. (2021) also indicated a negative correlation among AOD and EVI over Lahore, Multan, Faisalabad, Jacobabad, and Karachi. The highest positive r of 0.54 is observed over Madhya Pradesh, Maharashtra, India, and Afghanistan while negative correlation is observed over Nepal and Bangladesh. Anthropogenic emissions due to rapid urbanization, more fuel consumption, and biomass burning are responsible for high aerosol loadings over South Asia. Prasad and Singh (2007b) reported high AOD over the IGB region in the summer as compared to other seasons. The slope value of 0.0002 represents an increasing trend of monthly mean AOD and the y-intercept is found to be 0.29. Ali et al. (2014) investigated aerosol loading over Lahore and obtained high AOD in the summer season due to high temperatures and dust activity. Figure 9 demonstrates the monthly average AOD and their standard deviations over the study area from 2002 to 2020. The highest AOD (~ 0.44) is detected in June while the second-highest AOD (~ 0.43) is seen in July. The minimum AOD (~ 0.22) is obtained in December. The relative change between January and June AOD is ~ 76% as shown in Fig. 9.

Temporal variations in AOD
This depicts the impact of the high temperature during the summer months over the study region. The interannual variability of AOD is higher in the summer season than in winter as noticeable from the standard deviation of AOD as shown in Fig. 9. Tariq et al. (2021) also reported high AOD over Pakistan in summer as compared to the winter season.
The Aqua-AOD measured at 550 nm ranging from 0.02 to 0.66 with a mean value of 0.34 ± 0.16 over Pakistan, 0.08 to 0.72 with a mean of 0.36 ± 0.13 over India, 0.26 to 0.57 with a mean of 0.45 ± 0.11 over Bangladesh, 0.02 to 0.55 with mean of 0.25 ± 0.15 over Nepal, 0.05 to 0.07 with mean of 0.06 ± 0.006 over Bhutan, and 0.24 to 0.33 with mean of 0.27 ± 0.03 over Sri Lanka as shown in Table 1. The descriptive statistics show that the aerosol loading was relatively higher over Pakistan, India, and Bangladesh than Nepal, Afghanistan, Bhutan, and Sri Lanka. Figure 10 depicts the time series of MODIS-Aqua retrieved AOD over Kabul, Lahore, Karachi, New Delhi, Bangalore, Varanasi, Kathmandu, and Dhaka from September 2002 to December 2020. The slope value of − 2E-06 depicts the decreasing trend of AOD over Kabul with an intercept value of 0.5915 as shown in Fig. 10a.
The highest AOD of ~ 0.28 over Kabul was observed in March 2012. Figure 10b shows the increasing trend of Aqua-AOD over Lahore with a slope of 2E-05 and intercept of 0.5915 while the highest AOD of ~ 1.36 was observed in November 2014. The natural activity contributing to high AOD over Lahore is dust storms with significant increase in temperature and high wind speed. The slope value of 5E-06 represents an increasing trend of AOD over Karachi as shown in Fig. 10c. The highest AOD of ~ 1.0 over Karachi was observed in August 2010. Dust aerosols from Cholistan and Arabian Sea contribute to high AOD over Karachi along with local vehicular emissions (Alam et al. 2014). In New Delhi, the positive slope value of 0.0006 shows the increasing trend of AOD while the highest AOD (~ 1.4)  Figure 10d shows slightly high AOD over New Delhi throughout the study period due to the high aerosol loading over the city from the IGB region. The long-range transport of dust aerosols from Middle East and the Thar Desert contributes to aerosol loadings over New Delhi (Tiwari et al. 2016). This increase in AOD affects the climate of New Delhi badly for the last few years. Panday and Prinn (2009) found increasing trends of AOD over New Delhi during the pre-monsoon and winter because of dust activity and burning of biomass.
The highest AOD of ~ 0.51 over Bangalore was observed in May 2018 as shown in Fig. 10e. The slope value of 3E-05 and intercept of 0.119 represent rising trend of AOD over Bangalore. Ramachandran et al. (2012) found an increasing trend of AOD over Hyderabad and Bangalore due to the rapid increase in urbanization. Sreekanth (2013) reported presence of marine aerosols over Bangalore due to the presence of high AOD and low values of fine mode fraction. Figure 10f shows the increasing trend of AOD over Varanasi, India. The slope and y-intercept were 4E-05 and  Fig. 10g. Previous studies showed that the katabatic winds and topography of Kathmandu cause pollutants to be trapped under inversion layers near the valley (Panday and Prinn 2009). Figure 10h shows the increasing trend of AOD over Dhaka with the highest AOD of ~ 1.2 in June 2018. The y-intercept and slope from the linear model are 0.4204 and 4E-05, respectively. The main sources of natural aerosol emissions over Bangladesh are desert dust from Indian region and sea salt from Arabian Sea (Mamun 2014). Figure 11 shows the time series of annual average Aqua-AOD over Kabul, Lahore, Karachi, New Delhi, Bangalore,      Fig. 11. Tiwari et al. (2016) reported that the aerosol mixtures present in the atmosphere of Delhi are major causes for the reduction of air quality and destruction of climate. In Bangalore, the AOD has several change points starting from 2002 to 2015 later on it rises till 2017 and then decreases onwards. Ramachandran et al. (2012) have increasing annual AOD trend over most of the Indian cities including Bangalore using 10 years datasets of MODIS. They reported increase in AOD of ~ 0.013 per year due to rapid increase in vehicles number and urbanization. The Varanasi, Kathmandu,

Conclusion
This study shows the mean spatial distribution of AOD over South Asia and its association with EVI and meteorological variables from September 2002 to December 2020. The seasonal spatial distribution of AOD, seasonal anomalies of AOD, mean monthly and annual trends of AOD over selected cities, and mean monthly trends of AOD over South Asia were also examined. The mean annual highest AOD has been observed over the IGB region due to an increase in economical and agricultural activities while the lowest AOD was observed over northwestern Pakistan due to the less population density and limited agricultural activities. The highest AOD was observed during the monsoon season over IGB, Bangladesh, and southern Afghanistan region. The main reason for high AOD over IGB and eastern Pakistan is dust storms. High AOD was also observed during the winter season over IGB region because of the anthropogenic emission while low AOD levels were seen over Pakistan and Sri Lanka as compared to the monsoon season. The highest seasonal anomaly value is observed over eastern Pakistan and western India during winter. In eastern Afghanistan, southwestern Pakistan, and western Nepal, the AOD is highly correlated with temperature. The highest positive correlation among AOD and RH is noticed over northeastern Pakistan and central India. The negative correlation of AOD and RH is seen over Afghanistan, southwestern Pakistan, Bhutan, Nepal, and Bangladesh. The AOD is negatively correlated with EVI over Nepal and Bangladesh while positively correlated over Madhya Pradesh and Maharashtra India and Afghanistan. The positive AOD trends have been observed over Lahore, Karachi, New Delhi, Bangalore, Varanasi, Kathmandu, and Dhaka during the period September 2002 to December 2020. The highest annual change in AOD of 105.25% over Bangalore while the lowest change of − 12.12% is observed over Karachi from 2002 to 2020.