Next Article in Journal
High-Accuracy Real-Time Kinematic Positioning with Multiple Rover Receivers Sharing Common Clock
Next Article in Special Issue
Spatiotemporal Investigations of Multi-Sensor Air Pollution Data over Bangladesh during COVID-19 Lockdown
Previous Article in Journal
Diurnal Cycle of Passive Microwave Brightness Temperatures over Land at a Global Scale
Previous Article in Special Issue
Inferring Near-Surface PM2.5 Concentrations from the VIIRS Deep Blue Aerosol Product in China: A Spatiotemporally Weighted Random Forest Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Interdecadal Changes in Aerosol Optical Depth over Pakistan Based on the MERRA-2 Reanalysis Data during 1980–2018

1
Collaborative Innovation Centre on Forecast and Evaluation of Meteorological Disasters, Key Laboratory of Meteorological Disaster, Ministry of Education (KLME), International Joint Laboratory on Climate and Environment Change (ILCEC), Key Laboratory for Aerosol-Cloud-Precipitation of China Meteorological Administration, School of Atmospheric Physics, Nanjing University of Information Science and Technology, Nanjing 210044, China
2
Department of Physics, Higher Education, Government of Khyber Pakhtunkhwa, Peshawar 25000, Pakistan
3
Department of Physics, Koneru Lakshmaiah Education Foundation (KLEF), Vaddeswaram, Guntur, Andhra Pradesh 522502, India
4
School of Geographical Sciences, Nanjing University of Information Science and Technology, Nanjing 210044, China
5
Royal Netherlands Meteorological Institute (KNMI), R & D Satellite Observations, 3731 GA De Bilt, The Netherlands
*
Author to whom correspondence should be addressed.
Submission received: 22 December 2020 / Revised: 14 February 2021 / Accepted: 16 February 2021 / Published: 23 February 2021

Abstract

:
The spatiotemporal evolution and trends in aerosol optical depth (AOD) over environmentally distinct regions in Pakistan are investigated for the period 1980–2018. The AOD data for this period was obtained from the Modern-era retrospective analysis for research and applications, version 2 (MERRA-2) reanalysis atmospheric products, together with the Moderate-resolution imaging spectroradiometer (MODIS) retrievals. The climatology of AODMERRA-2 is analyzed in three different contexts: the entire study domain (Pakistan), six regions within the domain, and 12 cities chosen from the entire study domain. The time-series analysis of the MODIS and MERRA-2 AOD data shows similar patterns in individual cities. The AOD and its seasonality vary strongly across Pakistan, with the lowest (0.05 ± 0.04) and highest (0.40 ± 0.06) in the autumn and summer seasons over the desert and the coastal regions, respectively. During the study period, the annual AOD trend increased between 0.002 and 0.012 year−1. The increase of AOD is attributed to an increase in population and emissions from natural and/or anthropogenic sources. A general increase in the annual AOD over the central to lower Indus Basin is ascribed to the large contribution of dust particles from the desert. During winter and spring, a significant decrease in the AOD was observed in the northern regions of Pakistan. The MERRA-2 and MODIS trends (2002–2018) were compared, and the results show visible differences between the AOD datasets due to theuseof different versions and collection methods. Overall, the present study provides insight into the regional differences of AOD and its trends with the pronounced seasonal behavior across Pakistan.

Graphical Abstract

1. Introduction

With the rapid urbanization and economic growth, atmospheric aerosols and their effects on the earth’s climate system have become a growing concern over the past few decades [1]. Aerosols influence atmospheric stability and climate change directly through interaction with solar radiation and indirectly through their effects on clouds and the regional hydrological cycle [2]. They also influence the atmospheric environment [3,4], surface temperature, and human health [5]. The atmospheric residence time of aerosols varies from hours to days, depending on particle size and removal processes, and their concentrations vary strongly due to sources and sinks, including atmospheric processes [6]. However, the effects of aerosols are still one of the largest uncertainties in climate research [1].
The most common property of atmospheric aerosols related to their column-integrated amount is the aerosol optical depth (AOD). AOD is the integral of the aerosol extinction coefficient over the whole atmospheric column, i.e., “an optical property, which can be retrieved from measurements of radiation using satellite- and ground-based sensors”. AOD is considered a crucial parameter for climate change assessment and long-term analysis of aerosols on local, regional, and global scales [6,7,8,9,10,11,12,13,14]. Since the end of the 1970s, “with the development of remote sensing technology (TOMS (total ozone mapping spectrometer) and AVHRR (advanced very high-resolution radiometer))”, spaceborne aerosol products have been developed to investigate the direct radiative effect of aerosols, spatial and temporal distributions, and their trends on various scales (spatial and temporal) across the world [15,16,17,18,19,20,21,22,23]. As opposed to ground-based observations (which provide detailed information with limited spatial representativeness), space-borne sensors are capable of providing knowledge about theaerosol spatiotemporal distributions and trend variations over large areas in a continuous manner, with considerably high temporal resolution. However, satellite-based retrievals of AOD are now much more accurate compared to the early versions of the products. In addition, the satellite retrieval of aerosol properties is now possible from observations in all-sky conditions (e.g., Aura/OMI sensor) developed in the last decade [24]. The Modern-era retrospective analysis for research and applications, version 2 (MERRA-2), is a reanalysis tool, incorporating satellite and model data, which provides information on aerosol concentrations, composition, and optical properties on regional to global scales with complete temporal and spatial coverage [23,25,26].
MERRA-2 data havebeen widely employed so far to investigate the long-term annual and seasonal trends of AOD at global and regional scales (See Table 1) due to its wide spatial and/or temporal coverage and long-term variations for the ground-based and satellite validation. Klingmüller et al. [27] studied regional AOD trends from MODIS data for the period 2001–2012 using the multiple linear regression statistical technique over different parts of the Middle East and revealed regional differences in the trend. Che et al. [28] analyzed AOD trends using MERRA-2 reanalysis data at the global level and found a significant decrease in AOD trend over Europe and Eastern US, and increasing trends were noticed over China and South Asia during the most recent decades.
The main focus of this study is to describe the spatial and temporal variations of aerosols over Pakistan, which, located in Southeast Asia (Figure 1), is one of the most important regions globally, bit, which have been explored to a minimum extent. Pakistan experiences a growing population, urbanization, and/or industrialization and has a unique topography covering different climate regions [15,29,34]. The aerosol properties in Pakistan have been documented well from several studies [14,16,29,34,35,36] by using ground and space-borne observations in different periods. For example, an earlier study by Alam et al. [15] reported an increasing trend in AOD from 1979 to 2001 for most of the regions in Pakistan using data from the Total ozone mapping spectrometer (TOMS). Gupta et al. [30] noticed an overall decreasing trend (AOD) in Lahore with a negative slope (−0.0006) during the study period 2001–2010. Kumar et al. [9] analyzed the AOD climatology over the Indo-Gangetic Plain (IGP) during 2006–2015, utilizing MODIS collection 6 (C6) level 2 Deep Blue (DB) products [37,38,39,40,41]. Gupta et al. [30] and Kumar et al. [9] reported an increasing trend in AOD over urban stations in Pakistan, such as Lahore, Multan, and Karachi. Khan et al. [29,34] extensively investigated the AOD using the ground-based aerosol robotic network (AERONET) sunphotometer data over the megacities, Lahore and Karachi, between 2001 and 2018 and analyzed implications of different aerosol types for radiative forcing. Their results showed a positive AOD trend in Lahore of 0.013 year−1 and a negative AOD trend in Karachi of −0.02 year−1. They also found a close relationship between AOD and regional meteorology. For instance, in a recent study using the long term (2000 to 2016) MODIS and MISR data, Khalid et al. [42] found that Western Pakistan and Eastern urban regions experienced a decreasing AOD trend in the summer and an interannual uniform increasing trend for the Eastern route of the CPEC regions (Pakistan). Recently, Ali et al. [14] explored the spatiotemporal variations of AOD trends over urban regions in Pakistan using the MODIS-Terra C061 DB products for 2003–2017. These authors observed trends varying between 0.0002 year−1 and 0.0047 year−1, which were attributed to the dynamic role played by meteorology in the region. In short, although the trend analysis was applied for spatiotemporal time-series in earlier studies, so far, they focused on the data recorded in the major cities (Lahore and Karachi) of Pakistan. A long-term AOD climatology over Pakistan is still lacking, especially in the major urban regions. Therefore, there is a growing concern for large-scale spatiotemporal analysis of aerosol optical property (AOD), which became the main motivation of undertaking the analysis in the present work.
In the present paper, the multi-decadal changes in AOD are investigated using 39 years of MERRA-2 reanalysis data (1980–2018) together with MODIS-Aqua (2002–2018) data. TheseAOD data sets are compared for the overlapping period. The seasonal and annual variations of AOD on regional and local scales are analyzed with different characteristics for the selected cities and regions located in Pakistan. The spatial variations and interannual trends in AOD are also investigated in connection with the local and regional emission sources.

2. Data and Methods

2.1. Study Region

Pakistan is a country located in Southeast Asia, with a large spatial variability of meteorological conditions, topography, and population. The different regions are characterized by differences in population, industrialization, geography, and meteorological conditions, with multifarious and complex topography. However, the division in different sub-zones makes it possible to analyze the spatiotemporal climatology of the AOD over Pakistan. The spatiotemporal variation of satellite and model-based AOD wasderived over twelve selected cities located in six different geographical regions across Pakistan (Figure 1) to account for a comprehensive spatial and temporal coverage (including estimation of trends), seasonal and annual heterogeneities in aerosol properties, emission sources, and meteorological conditions.
The cities were chosen in such a way that they represent the regional climate and emission sources across the entire study domain. The cities considered are Gilgit (Gil) and Muzafarabad (Mzf), located in the North Dry Region (NDR); Peshawar (Psh) and D.I.Khan (DIK) in the Katawaz Basin; Lahore (Lhr) and Multan (Mlt) in the Indus Basin region; Karachi (Kr) and Chhor (Cr) in the coastal region; Quetta (Qt) and Sibbi (Sb) located in the BalochistanPlateau; and Khuzdar (Kz) and Panjgur (Pj) in the Kharan Desert Region (KDR). Figure 1a represents the various land cover classes across the study domain. Information on all cities related to the geographical environment, i.e., location, elevation, population density, precipitation, and climate, is presented in Table 2. Figure 1b shows the location of the main study domain (Pakistan) in south Asia, with its bordering countries and potential source desert regions, including the directions of air masses towards the study region.
Gilgit (Gil; 35.92° N, 74.33° E) and Muzafarabad (Mzf; 34.37° N, 73.48° E) are located in the North Dry Region (NDR), which is considered the least polluted rural site in Pakistan. The main sources of aerosols in this cold and semi-arid climate zone include soil/road dust, transportation, and anthropogenic emissions from vehicular and industrial sectors. Peshawar (Psh; 34.01° N, 71.52° E) is the capital city of Khyber Pakhtunkhwa province with a tropical climate year-round. There is a variety of aerosol sources, with the dominance of dust during the spring (MAM) and fine anthropogenic aerosols in the autumn (SON) and winter (DJF) seasons [16]. D.I.Khan (DIK; 30.03° N, 70.64° E) is located close to the desert region; therefore, dust aerosols are persistent [15], along with aerosols emitted from local sources. Multan (Mlt; 30.19° N, 71.47° E) is situated in the central part of the Punjab province, with vehicular emissions and crustal sources as the major contributions to aerosols [15].
Strong seasonality in the aerosol-loading and type was reported in Lahore (Lhr; 31.54° N, 74.32° E), including vehicular and industrial emissions, biomass burning, transported specks of dust, and fossil fuel combustion-related particles, recognized as the principal contributors to the aerosol column load [17,27,28] in the recent past. Quetta (Qt; 30.18° N, 67.01° E) is the capital city of Baluchistan province, which experiences high aerosol-loading during the spring and winter seasons [43], originating from both continental and anthropogenic sources as thermal power plants, brick kilns, and coal-rich regions of the country [43]. Sibbi (Sb; 29.33° N, 67.55° E) is a mountainous region in the Baluchistan foot belt, with comparatively less pollution; major sources include soil dust and vehicular/industrial emissions. Karachi (Kr; 24.87° N, 67.03° E) is a coastal city in a typical urban setting with a variety of pollution sources, including industrial and vehicular emissions, dust storms, and sea spray aerosol [34]. Chhor (Cr; 25.31° N, 69.47° E) is an arid, suburban region located in southern Pakistan. It is considered a hotspot of dust aerosols because of its proximity to the Thar and Cholistan deserts, as well as the intrusion of sea spray aerosols from the nearby Arabian Sea. Both Khuzdar (Kz; 27.5° N, 66.38° E) and Panjgur (Pj; 26.58° N, 64.10° E) in the Balochistan province are located at high elevation (4058 ft), with a subtropical hot desert climate along with biomass burning and marine aerosols. Khuzdar and Panjgurdiffer from each other mainly by elevation, population, and sources of pollutants.

2.2. The AOD Data

2.2.1. The MERRA-2 Reanalysis Data

The monthly AOD data over Pakistan were collected from the latest MERRA-2 atmospheric reanalysis data for the period 1980–2018, created with an emphasis on the assimilation of various satellite measurements. MERRA-2 uses the new version of the Goddard Earth Observing System Model, version 5(GEOS-5) data assimilation system to synthesize regular time-series of gridded data, with a spatial resolution of 0.5° × 0.625° (latitude × longitude) at 72 pressure levels (from the surface up to 0.01 hPa) both instantaneous and time-averaged (hour, 3h, and month). Among the advantages of MERRA-2, in comparison with its previous version, are the joint assimilation of aerosol and meteorological fields, accounting for their interaction through the radiative effects of atmospheric aerosol [44]. The aerosol assimilation system in MERRA-2 uses AOD measurements from satellite instruments (MODIS, MISR; AVHHR only over the ocean) and also from the ground-based AERONET network stations. Over the land, the MERRA-2 is constrained by assimilation of MODIS dark target (DT) C5 AOD, and when the surface reflection exceeds 0.15, MISR AOD products are used [26]. MODIS and MISR are only available from early 2000, and hence MERRA-2 AOD over land for earlier years is not constrained by satellite data. MERRA-2 AOD was evaluated by comparison with ground-based sunphotometer data from AERONET and further with satellite data [26,45]. For the pre-EOS era, when satellite-retrieved AOD and AERONET data were not available, there are only a few alternatives, such a sun duration measurements. The comparison with sun duration measurements shows that MERRA-2 overestimates the AOD [46]. The better comparison with the sunphotometer, as compared with sun duration, data may be the effective cloud screening in the sun photometer processing chain, which biases the validation with sun photometer data to clear sky only. Likewise, the comparison with satellite data shows the underestimation of high AOD and the overestimation of low AOD [11]. The MERRA-2 reanalysis data are open access and free for the public (https://disc.sci.gsfc.nasa.gov/ (accessed on 1 December 2020)). An overview of the MERRA-2 modeling system and more detailed descriptions of aerosol treatment in the MERRA-2 system can be found in [26,44,45].

2.2.2. The MODIS Sensor

The MODIS instrument has been flying aboard the Aqua satellite, which is an important part of the NASA Earth Observing System, since July 2002. It provides several aerosol products over the ocean and land [18,40]. The sensor measures radiances at spatial resolutions of 0.25, 0.5, and 1.0 km depending upon the wavelength band for 36 spectral channels from 0.415 to 14.235 μm, with a viewing swath of 2330 km. The MODIS team provides three modes of Aqua-MODIS Collection-6.1 (C061) AOD products retrieved using the dark target (DT) [18] or deep blue (DB) [38] algorithms. The choice ofthe MODIS algorithm (DT or DB) used for aerosol retrieval depends on the properties of the underlying surface, following the criteria outlined in Levy et al. [18] and Sayeret al. [41]. For better coverage, the DT and DB AOD products in C060 are combined in a merged product [18,40], which in this paper is referred to as DTB. C061 is the latest version of the Collection of MODIS aerosol retrieval products. The MODIS-Terra C061 AOD product over China was evaluated by Sogacheva et al. [23], and the MODIS-Aqua C061 AOD product over Pakistan was evaluated by Ali et al. [14]. The C061 products for both MODIS sensors are significant improvements over earlier versions like C005, C051, and C006 (e.g., the increased AOD over land and reduced AOD over the ocean from Aqua). Overland, the uncertainty in the MODIS AOD is 0.05 ± 0.15 τ; over the ocean, it is 0.03 ± 0.05τ (Kaufman et al. 1997; [18,47], where τ is the AOD. In this study, the AOD at a wavelength of 550 nm (AOD550) Level 3 data obtained from the MODIS-Aqua Collection 6.1 (C061), DB product [39,40,41], is used with a spatial resolution of 1° × 1°. The MODIS AOD product was downloaded from http://modis.gsfc.nasa.gov/ (last accessed on 19 January 2020). In this manuscript, only AOD550 is used and, therefore, from here on, referred to as AOD.

2.3. Statistical Analysis and Modeling

For the detection of a statistically significant trend in the AOD time-series, the Mann–Kendall (MK) test can be used [48,49,50]. The MK test is simple and robust against outliers, missing values, normal distribution, and has low sensitivity to abrupt breaks in time-series [51]. However, the presence of autocorrelation would add a significant trend to the time-series, even when there is no trend, and affect the outcomes of the MK test [52,53]. Therefore, the occurrence of autocorrelation should be checked and removed before the application of the Sequential-MK (Seq-MK) test to the time-series data set [49].
To avoid the effects of serial correlation in the time-series on the results of the MK test [11], the Seq-MK technique can be applied [52,53]. Furthermore, the MK test was used to detect monotonic trends, whereas the Seq-MK test was adopted to detect a sudden change or shift in the trend over time (see Supplementary Material (SM)). The relationship between the AOD trend and elevation was assessed by linear regression. These statistical methods and models are briefly discussed in the following subsections.

2.3.1. Mann–Kendall Test

The MK test is a nonparametric method that can be used to identify a statistically relevant pattern in the AOD time-series [48,50]. The MK test has low sensitivity to sudden breaks in time-series and is capable of outliers, incomplete values by checking the regular distribution of data [51]. The test’s null hypothesis (H0) claims that there is no pattern over time in the AOD time sequence, while the alternate hypothesis (H1) suggests that there is a monotonous pattern over time (increasing or decreasing). The mathematical equations for computing Mann Kendall statistics is defined by (Hirsch and Slack [50] in Equation (1)
S =   i = 1 n 1 j = i + 1 n sgn   x j   x i
where sgn (xj-xi) can be computed by Equation (2)
sgn   x j   x i = + 1 ,   if   x j x i > 0 0 ,   if   x j x i = 0 1 ,   if   x j x i < 0
where Xj and Xi are the observed time-series of length n. For identically and independently distributed datasets with zero mean, the variance St of the statistics S is given by Equation (3) [54].
S t =   1 18 n n 1 2 n + 5 p = 1 q t p t p 1 2 t p + 5
where tp is the number of ties for the pth value, and q shows the number of tied values. The MK test statistic, z, can be computed using Equation (4):
z =   s 1 s t   if   s t > 0 0   if   s t = 0 s + 1 s t   if   s t < 0
where statistical z values beyond the confidence interval of ±1.96 suggest that the phenomenon has a statistical significance of 95% [52,53]. Positive z values indicate a positive trend and vice versa.

2.3.2. Autocorrelation Test

Seasonality typically impacts the time-series of atmospheric parameters and should be removed in the analysis of long-term patterns. Autocorrelation or serial correlation is a statistical tool used in testing and interpreting time-series data [55]. Further, autocorrelation analysis is suitable to understand and check on serial correlation caused by the seasonality within a long-term data set, such as AOD. The following steps were adopted to check on autocorrelation in different time-series.
I. Calculate the coefficient of lag-1 serial correlation of the different time scales concerning the AOD series as provided in Equation (5)
r 1 =   1 n 1 i = 0 n 1 ( x i     x ¯ ) ×   x i + 1     x ¯ / 1 n i = 0 n 1 x i       x ¯ 2
where r1, xi, and x ¯ are the correlation coefficients of lag-1, AOD time-series, and its mean, respectively.
II. The confidence interval (CI) of the correlation coefficient lag-1(r1) at the 5% significant level can be computed using Equation (6) [4].
r 1   ( 5 % )   =   ( 1   ±   1.96   n   2 ) /   ( n 1 )
Moreover, the MK check can be implemented after the serial connection has been eliminated [56], which can be computed as:
(x2 − r1 x1, x3−r2 x2 ….xn − rn−1xn−1)

2.3.3. Linear Regression

Trends in AOD time-series and population were analyzed by linear regression expressed in Equation(8):
Y(x) = ax + b
where Y(x) is the AOD value at each value of the variable x, a and b are the slope and intercept, respectively, of the linear regression.

3. Results

3.1. Comparison of Datasets

To evaluate the quality of the MERRA-2 AOD over Pakistan, the data were compared with the three Aqua-MODIS C6.1AOD products (DB, DT, and DTB) for the years 2002–2018. For this comparison, the MERRA-2 data (0.5° × 0.625°) were re-gridded to the equivalent resolution (1° × 1°) of the MODIS-Aqua monthly AOD data. The scatter plots of the MERRA-2 AOD versus MODIS-Aqua AOD are presented in Figure 2a–c.
The results show the overall agreement between the MERRA-2 and the three MODIS Aqua AOD products. For the comparison with DB the bias is very small (0.015) and other statistical metrics are R = 0.88, slope = 0.718 and RMSE = 0.056. For the comparison with DTB the bias is 0.04 and the other metrics are R = 0.875, slope = 0.632 and RMSE = 0.052. The comparison with DT is less favorable with a higher bias of 0.22 and statistical metrics are R = 0.857, slope = 0.548 and RMSE = 0.051. A fair comparison with the DB AOD reflects the accuracy of the MERRA-2 data set (Since MERRA-2 assimilates MODIS AOD), which is consistent with the recent findings of Tsidu et al. [57] over the African continent.
The degree of auto-correlation of the MERRA-2 reanalysis data was checked by the application of the technique described in Section 2.3. It is important to mention here that the incomplete sampling from the retrieval source (satellite/ground) or data set induces biases in the long-term trend analysis [52]. The autocorrelation technique is capable of analyzing values of repeated data in a series data set and is based on discrete-time-series values [58]. The autocorrelation coefficients of the annual mean AOD time-series in the six regions are presented in Figure 3 as a function of the lag (for lag up to 20). The autocorrelation is most prominent over the Indus Basin, the Katawaz Basin, the KDR, and the Balochistan Plateau. In contrast, over the NDR and the Coastal region, the fluctuations indicate that the autocorrelation is not significant. Therefore, the presence of lag 1 and autocorrelation onward (lags) at different sites was removed, and seasonality was nullified overall study sites with the help of a pre-whitening approach to make the data suitable for further analysis (Seq-MK test). More detail about this approach (pre-whitening) is given in the studies by Von Storch. [59] and Douglas et al. [60] and are not given herein to avoid repetition.

3.2. Spatiotemporal Distributions of AOD

3.2.1. Spatial Variations

The spatial and seasonal variations of the AOD over the study domain were studied using the long-term MERRA-2 (1980–2018) and MODIS-DB (2002−2018) data sets. The maps of the annual and seasonal mean AOD overall years during the overlapping period are presented in Figure 4. The spatial distributions of the MERRA-2 and MODIS AOD are similar; however, the MERRA-2 AOD values are lower than those from MODIS, which is attributed to the spatial resolution of the MERRA-2 AOD data. The aerosol patterns for each season are similar, with high AOD observed in the east over the Indus basin and coastal regions and low AOD in the west and north of the study domain. The spatial variation of the aerosols across the study area is related to the occurrence of distinct aerosol sources (natural and anthropogenic) (see Figure S4 of SM).
The seasonal variation is due to meteorological effects and anthropogenic emission sources associated with the domestic use of energy, agricultural practices, biomass burning, natural seasonal emissions, such as dust, and long-range transport [15,17,34]. The strongest seasonal variation is observed in the East of Pakistan, with the highest AOD(>0.8) over the Indus basin and coastal area during the summer and the lowest AOD in the spring (~0.45). The seasonal variation is different between the Indus Basin and the coastal regions, indicating different reasons. In the coastal regions, where transport of aerosols and gases are not obstructed by surrounding mountains, the occurrence of high AOD is attributed to the enhanced atmospheric circulation (driven by solar radiation), which affects the regional aerosol-loading from different major dust sources (Central Asia, northern hemisphere, and the Thar Desert).
In the Katawaz Basin and the Balochistan Plateau, the AOD varies between 0.2 and 0.5 for the study period in different seasons, with the lowest ( <0.1) in winter. The reason is the stable atmosphere with weak vertical mixing, from which aerosol particles are removed relatively fast (Figure 4). However, the additive increment of aerosol-load during the summer over the Katawaz Basin and the Balochistan Plateau indicates the transport of pollution due to strong convective winds; including dust and long-range transport of other aerosol types (see Figure S4 of SM), which leads to an increase in column aerosol load. Nearly identical results were previously observed by Alam et al. [15,16] based on satellite data over different regions of Pakistan. However, the AOD values are substantially larger than the values reported by Kumar et al. [9] for the cities Karachi, Lahore, and Multan in the IGP region. The AOD was low ( <0.1) over the rural and less populated areas of the NDR throughout the study period, which are attributed to its high elevation, regional meteorology, the absence of anthropogenic emission sources, and low population density.

3.2.2. Temporal Variations

Monthly Variation

Time-series of monthly mean AOD over all the selected study cities averaged annually for the period 2002–2018 are presented in Figure 5 for both data sets (MERRA-2, MODIS). The highest AODMERRA was found during July over the DIK (0.82 ± 0.13) followed by Lahore (0.82 ± 0.16), Multan (0.83 ± 0.13), Karachi (0.75 ± 0.08), Chhor (0.74 ± 0.09), and Khuzdar (0.56 ± 0.07). This is attributed to enhanced anthropogenic activities leading to the increased emissions (smoke from biomass burning and seasonal local/regional transport from desert dust) during the summer. Further, in some cities in the study domain, the MERRA AOD is a factor 2–3 higher than the MODIS AOD, which is attributed to spatial resolution and uncertainties associated with each retrieval source (discussed in detail in Section 2).
Notably, the lower AODMERRA is evident over the NDR and Balochistan Plateau, especially during the winter (December–February) with mean values of 0.07 ± 0.01, 0.1 ± 0.01 in Gilgit and Muzaffarabad, followed by Quetta (0.13 ± 0.03), Khuzdar (0.13 ± 0.02), and Panjgur (0.12 ± 0.01), respectively. These lower values are attributed to the high elevations of these sites as compared to the other sites included in this study (Table 2, Figure 5). However, the lower AOD found at Sibbi is mainly interlinked with its smaller population and comparatively clean environment (Table 2). The time-series analysis of both AOD data sets (MODIS and MERRA-2) indicates almost similar patterns in individual cities, but with much smaller AODMODIS than AODMERRA, which may be due to the coarser spatial resolution of MERRA-2 than that of MODIS.

Seasonal and Annual Variations

The seasonal and annual MERRA-2 AOD, averaged over all years during the study period (1980–2018), for the six study domains in Pakistan are presented in Table 3, and the trends are presented in Table 4. The annual mean AOD is highest in the economically (highly urbanized) and industrially developed coastal area (0.41 ± 0.14), followed by the Katawaz Basin (0.40 ± 0.15), KDR (0.38 ± 0.12), and Indus Basin (0.33 ± 0.11), with annual increases of 0.49%, 2.29%, 1.56%, and 1.87%, respectively (Table 4). The lowest annual mean AOD is observed in the rural and less developed areas of the NDR (0.12 ± 0.04) and regions with a low population density (Balochistan Plateau, 0.24 ± 0.07), with an annual increase of 0.29% and 2.12%, respectively (Table 3, Table 4 and Table 5).

3.3. Trend Analysis of AOD

3.3.1. Temporal Trend

To investigate the interannual, decadal variation in aerosol-loading for different cities across the study domain, AOD trends were computed using the Seq-MK test with linear statistical regression (Figure 6 and Figure 7, Table 6). Strong seasonality in AOD was observed in the coastal region and the Kharan Desert Region (KDR), with high values during the summer (0.54 ± 0.13) followed by KatawazBasin (0.52 ± 0.16) and Indus Basin (0.45 ± 0.12). The seasonally averaged AOD was low in the winter over the North Dry Region (NDR; 0.08 ± 0.05), followed by Balochistan Plateau (0.14 ± 0.06) and Indus Basin (0.23 ± 0.12). Seasonally, the occurrence of high AOD during the summer over the Indus basin is mainly associated with the local and long-range transported dust aerosol from the Thar Desert in the East, the Kharan Desert in the West, and the Taklimakan Desert in the Northwest of Pakistan [14]. In addition, the atmospheric solar radiation and vertical uplifting of soil particles from the surfaceby high winds play a significant role in column aerosol-loading in the regional atmosphere during the summer and autumn seasons [16].
Generally, the annual mean AOD time-series showed a decreasing trend in all the urban and industrial sites during the first study period (1980–2000). For most cities, a sudden increase in the AOD occurs around 2001 with maximum slopes of 0.0050, 0.0042, 0.0041, 0.0040, and 0.0032 year−1 in Khuzdar, Sibbi, Quetta, Multan, and Panjgur, respectively, except in DIKhan, where the AOD decreased with−0.001 year−1. During 1985, 1991–1996, and 2012, an abnormal increase of AOD is evident in various regions of the study domain, probably due to large amounts of smoke (e.g., from peat and forest fires) transported to other Asian regions from the surroundings [61]. Very small positive trends were observed in almost all cities during the first half of the study period (1980–2000), with a slope of 0.001 year−1 associated with increases inpopulation density and pollution sources. Furthermore, in other cities (such as Karachi, Khuzdar, Chhor, and Panjgur), gradually decreasing trends were observed of −0.001 and −0.002 year−1.
Further, increasing trends are observed acrossall stations located on the edge of the Basin regions (Indus and Katawaz) and the Balochistan Plateau, with thelargest values in Multan (0.0094 year−1) and Peshawar (0.0083 year−1), which are attributed to emissions from natural and anthropogenic sources and regional dynamics (associated with the local meteorology) (Figure 6, Table 4). The observed results and trends are consistent with the conclusions in the recent reports by Ali et al. [14] over Karachi, Lahore, and Multan, but which are slightly larger than the findings of Kumar et al. [9]. In contrast, Gupta et al. [30] found quite different (reverse) trends from their analysis of the long-term climatology based on MODIS data retrieved from 2001 to 2010. However, Gupta et al. [30] used MODIS C5 L2 DT AOD data, whereas we have used C6.1 DB AOD data in the present study. The differences between recent MODIS data sets, including C5 DT, C6, and C6.1, were outlined by, e.g., de Leeuw et al. [21] and Sogacheva et al. [23]. The differences between DT, DB, and DTB for C6.1 over Pakistan are illustrated in Figure 2, showing the better quality of DB. Hence, it is not surprising that differences occur between the current results and those from Gupta et al. [30].
Figure 7 shows the time-series of the seasonal mean AOD for the six regions during 1980–2018. Similar to several cities in Figure 6, some of the regional time-series also show initially fluctuating AOD with a sudden rise in 2001, after which the AOD remains high (that is the reason for separating the two-time-series in 2001). Generally, a seasonally positive trend is noticed in all the regions over the multi-decadal study period. However, a decrease in trend is observed at almost the end of the study period, such as in the summer for all regions, except the NDR. The overall increasing trend occurred during the autumn for all study regions, with apositive slope (percentage deviation) of 0.012 year−1 (3.51%) in the Katawaz Basin and 0.011 year−1 (3.23%) in the Balochistan Plateau (Table 4). In the winter, the trends over the Indus Basin and KDR region were 0.0065 year−1 (2.87%), followed by KDR with 0.004 year−1 (1.09%) in spring. Further, the lowest trend of 0.00064 year−1 (0.19%) is observed during the summer over the coastal region. Mehta et al. [62] reported a significantly increasing AOD trend during the winter over the Indian subcontinent, including Pakistan, which is consistent with our results (in the Indus Basin and the Coastal region).

3.3.2. Spatial Variation of AOD Trends Across Pakistan

The AOD trends and their corresponding significance levels are estimated using the linear regression nonparametric Seq-MK statistical technique. The map of the trends (in% per year) was determined during the years of the study 1980–2018, 1980–2000, and 2002–2018 (MERRA-2) and for the period 2002–2018 (MODIS-Aqua), respectively. The results are shown in Figure 8a–d, where the pixels at 95% significant level (p < 0.05) are marked with a black dot on both the annual and seasonal maps. In the first period (Figure 8a), the MERRA-2 AOD showed a significant upward trend (p < 0.05) observed in the southeast and southern coastal domain of Pakistan that was related to the dust transferring from neighboring Thar’s desert almost round the year. The seasonal AOD trends had no obvious differences in terms of their spatial pattern in autumn and winter compared with the annual trends. In the second period (1980–2000), the annual and seasonal AOD trends showed no significant downward trend over the country, including northwest and south coastal areas, while, during the third period 2002–2018, the spatial distribution of annual and seasonal (winter, spring) trend had a greater increase (significant) compared with the second period. In the summer and autumn (Figure 8c), no significant downward trend could be observed in almost the whole country. Furthermore, a significant positive trend was noticed seasonally (winter, spring) in the southeast of Pakistan. “Thal” and “Thar” deserts and the Gangetic plain are the sources of dust aerosol in this area. In the fourth period 2002–2018, MODIS-DB AOD showed a slightly different spatial pattern comparatively (Figure 8d). The seasonal AOD trend had obvious differences in terms of spatial distribution in the spring and autumn, compared with the annual trend. Unlike in other periods, there was no significant upward and downward annual/seasonal AOD trend in the study region. Using a different version of MODIS (C006) AOD in this study can be one of the reasons for this difference, as the MERRA-2 assimilated MODIS-5 version (as discussed earlier) is different from the one (MODIS-C006) used in this study (Figure 8d). However, adopting different measuring techniques [26], the number of sample data and others may also be the reason for this deviation.

3.4. Variation in AOD with Elevation and Population Density

To study the variations in AOD over time as a function of elevation and population density, the yearly mean AOD for each of the years 1980 to 2018 was averaged in longitudinal bins of 0.5° for each region. Figure 9 shows the variation of the annual mean AOD with longitude for each of the years 1980 to 2018 for the six study regions.
The data in Figure 9 clearly shows the strong increase in AOD around 2000 over the western part of the NDR, the IB, and the coastal region for all longitudes. Over the other three regions, the increase after 2001 is not as strong and occurred mainly in the eastern parts. The increase over the years is expected from the increased industrialization. Moreover, the findings (Figure 9) suggest that the increasing AOD in distinct study zones with changing topography is mostly associated with the varying tendency both in the distribution and change in longitude. These fluctuating shifts in aerosol distribution with elevation may also be related to the regional/topography, meteorology, and climate dynamics in the study region.
Figure 10 shows the time-series of the annual mean AOD in 2000–2018, in steps of five years (red data points), together with the population density (plotted along the secondary y-axis). The reason for the 5yearsteps is that population density is only available every 5 years. However, AOD has a steep increase from 2000 to 2005, which has the same reason as in the time-series discussed in Section 3.3. Further, the lowest annual change in AOD was observed in Balochistan Plateau throughout the study period (Figure 10).
For each region, the population density (ρ) increased smoothly with time. Whereas in the Indus and Katawaz Basins and the Balochistan Plateau, the AOD strongly increased between 2000 and 2003 and then continued to increase onward, but with a much slower range associated with the population increase and associated emissions. The population density (AOD) was found the highest in the Indus Basin with 550 km−2 (0.54 ± 0.13) followed by the Katawaz Basin with 340 km−2 (0.48 ± 0.40), the Coastal regionwith 225 km−2 (0.42 ± 0.30), NDR with 55 km−2 (0.28 ± 0.30), and the least over the KDR with 15 km−2 (0.36 ± 0.21). Figure 10 shows that an increase in population density results in emissions from anthropogenic activities, such as heating, cooking, transport, etc., which is accompanied by an increase in the AOD. However, the close relationship between both the parameters (AOD and ρ) is obvious, with a similar variation in most of the study areas (See Figure 10). In general, the increase in population density by x% leads to an increase in AOD by F(x)%. Hence other factors, such as industrialization and the type of industry, which often also leads to an increase in population and transport, results in increased emissions, which may be different from one city or region to another. Furthermore, changes in land use may lead to a shift in emissions and their effects on AOD.

3.5. Probability Distribution Function

To investigate where and in which region the AOD peaks most during the study period, probability distribution functions (PDF) of the monthly mean AOD over the six regions were plotted for 4 different periods, i.e., 1980–1990 (period1), 1991–2000 (period2), 2001–2010 (period3) and 2011–2018 (period4) (Figure 11). The results show that between the first 2 periods, there is not much change in the AOD distributions for most of the regions (b, c, d, f), but after 2000 both the peak and width of the PDFs changed substantially. The pdfs show a shift to higher AODs after 2000, i.e., low AOD values occurred less frequently, resulting in overall higher AOD during the whole year and much larger peak values after 2000. The trend continues in period4, resulting in higher peak values at somewhat higher AOD, but the highest values are not occurring more frequently, although the lower AOD values occur less frequently.
The behavior over NDR and Coastal region is somewhat different, with a quite large difference between periods1 and 2, wherein period 2 the pdf is much wider and a second mode occurs at AOD of 0.4, but the smaller AODs are still occurring as well (see Figure 11a,e). After 2000, there is a definite shift to higher AODs at the cost of the smallest values, and also the peak values are higher and larger AODs occur more frequently than in period1. For the coastal region, the pdf is bimodal, and the distributions are similar in periods 1 and 2, whereas, in periods 3 and 4, the distributions are substantially wider, with lower peak values and a shift to larger AOD. The two different modes indicate the presence of different AOD types. The high AOD values, especially in the recent decade, are interlinked with the natural haze and dust episodes that lead to distinct aerosol transport from deserts around the study area. Desert dust emissions peak during spring and late autumn [16]. Anthropogenic emissions associated with rapid urbanization and industrialization are the second reason for the changing PDF patterns, as most of the heavy industries are located in the southeast of Pakistan (Figure S4 of SM, Table 2).

4. Discussion

The spatiotemporal evolution and trend variations in aerosol optical depth (AOD) over environmentally distinct regions in Pakistan were presented in the above sections. The results show that the aerosol level increased since 1980, and the increase accelerated from the start of the 20thcentury. (It is not clear whether this is in part a result of the MODIS AOD assimilation, which is likely having an impact on MERRA-2 data product (i.e., post-2001)). Before that period, a sudden increase in AOD throughout the study regions can also be seen around 1993, which may be related to aerosols emitted from burning the oil wells during the Gulf war during 1991–1992. The time-series analysis of the MODIS and MERRA-2 AOD data shows similar patterns in individual cities. The strongest seasonal variation is observed in the East of Pakistan, with the highest AOD over the Indus basin and coastal area during the summer and lowest AOD in the spring, indicating different reasons. The high AOD in the Indus Basin is mainly attributed to the geographical situation of the basin region surrounded by mountains, which prevent horizontal transport of aerosols/precursor gases out of the study area, which thus accumulates in the region. In the coastal regions, where transport of aerosols and gases are not obstructed by surrounding mountains, the occurrence of high AOD is attributed to the enhanced atmospheric circulation (driven by solar radiation), which affects the regional aerosol-loading from different major dust sources (Central Asia, and the Thar Desert). During the summer, particle formation is enhanced by solar radiation and warm weather, and the growth of these new particles and their accumulation in stagnant weather (low wind speed) leads to an increase in the regional aerosol load [7,8]. Furthermore, the increase of the planetary boundary layer height (PBLH) during the summer results in high aerosol concentrations, which in turn results in elevated AOD. However, the additive increment of aerosol load during the summer over the Katawaz Basin and the Balochistan Plateau indicates the transport of pollution due to strong convective winds, including dust [36,42] and long-range transport of other aerosol types, which leads to an increase in column aerosol load. The low AOD found over the rural and less populated areas of the NDR isattributed to its high elevation, regional meteorology, the absence of anthropogenic emission sources, and low population density. During the study period, the annual increase in the AOD trend throughout the study region is attributed to the increases in population and emissions from natural and/or anthropogenic sources. However, the general increase in the mean annual AOD values over the central to lower Indus Basin is ascribed to the large emission of dust particles from the nearby deserts. Further, in some cities n the study domain, the MERRA AOD was a factor 2–3 higher than the MODIS AOD, which is attributed to uncertainties in both the data sets.
Further, the increasing trends noticed over all the stations located on the edge of the Basin regions (Indus and Katawaz) and the Balochistan Plateau are mainly caused by natural and anthropogenic sources and regional dynamics associated with the local meteorology. The observed results and trends are consistent with the conclusions in the recent reports by Kumar et al. [9] and Ali et al. [14] over Karachi, Lahore, and Multan. In contrast, Gupta et al. [30] found quite different (reverse) trends from their analysis of the long-term climatology based on MODIS data retrieved from 2001 to 2010. However, Gupta et al. [30] used MODIS C5 L2 DT AOD data; whereas, in the present study, we have used C6.1 DB AOD data. The differences between recent MODIS data sets, including C5 DT, C6, and C6.1, were outlined by, e.g., de Leeuw et al. [21] and Sogacheva et al. [23].
In addition, the AOD increases with increasing population density and changes in elevation. In the Indus and Katawaz Basins and the Balochistan Plateau, the AOD strongly increased between 2000 and 2002 was associated with the population increase and associated emissions. These parameters may not be independent, and the observed tendencies are likely due to the relations of these parameters with economic activity and associated population density resulting in domestic activities and transport. In addition, the aerosol-loading in the urban areas grows faster than that of the rural areas, which is similar to the conclusion obtained by Wang et al. [63] observed over East China.

5. Conclusions

The spatiotemporal variation of the AOD in different regions of Pakistan over the last four decades is studied using the MERRA-2 (1980–2018) gridded reanalysis and MODIS-Aqua DB (2002–2018) datasets. The study domain considered in this work covers different scales in Pakistan, from the whole country to six regions and 12 cities selected in these six regions. The obtained results agree well with the occurrence of natural phenomena and the evolution of population and anthropogenic activities over the study domain. The main conclusions drawn from the present study are summarized below.
  • The MERRA-2 and MODIS (DT, DB, and DTB) AOD datasets were compared, and the validation shows a strong positive correlation, especially with the MODIS DB AOD product. However, the MERRA-2 tends to underestimate MODIS AOD, especially at larger AOD levels;
  • The spatial distribution of AOD shows high concentrations over economically, and industrialized urban areas of the Indus Basin; low AOD is observed over cities in the NDR and Kharan Desert regions, which are high-altitude arid regions with the sparse population;
  • The AOD shows a clear seasonal variation in the coastal region and the Indus Basin with the highest AOD during the summer ( >0.6) and smaller values in the spring ( ~0.5), autumn (0.3), and winter ( <0.1);
  • The AOD changes with the elevation and increases with increasing population density. These parameters may not be independent, and the observed tendencies are likely due to the relation of these parameters with economic and industrial activity. In addition, some relation with population density may be due to domestic activities and transport phenomena;
  • The MERRA-2 and MODIS trends (2002–2018) were compared, and the results show differences between the AOD datasets as a result of using different versions and collection methods;
  • The annual and seasonal spatiotemporal trend analysis shows a statistically significant increase of the theAODMERRA-2 (at the 95% confidence level (p < 0.05) in all study regions, with the largest trends in the Katawaz Basin followed by Balochistan Plateau, Indus Basin, and KDR, with 2.29%, 2.12%, 1.87%, and 1.56% per year, respectively.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/2072-4292/13/4/822/s1. Figure S1: Correlation between MODIS-Aqua and MERRA-2-derived data set for the study period (2002–2018). The color scale is common for both AODs retrieved from MODIS and MERRA-2. Figure S2: Abrupt variation in annual mean AOD values during 1980–2018; over the study regions (a) NDR, (b)Katawaz Basin (c) Indus Basin (d) Balochistan Plateau (e) Coastal region and (f) KDR covering the study domain of Pakistan at a significance level of 95%. The statistics UF (solid red line) represent the trend variability in the forward direction, whereas UB (blue dot line) represents the trend in the backward direction with zero mean and unit standard deviation. Figure S3: Same as in Figure S2, but for the seasonal variations. Figure S4: Spatial variation of total natural/anthropogenic emission of black carbon (BC), organic carbon (OC), SO4, dust, and sea salt particles (unit: kg m−2) retrieved from MERRA-2 reanalysis during the study period 1980–2018 observed in Pakistan and the surrounding areas.

Author Contributions

Conceptualization, R.K., K.R.K., T.Z., G.d.L.; methodology, R.K. and W.U.; software, R.K.; validation, R.K. and K.R.K.; formal analysis, R.K. and W.U.; investigation, R.K.; data curation, W.U.; writing-original draft preparation, R.K.; writing-review and editing, K.R.K., T.Z. and G.d.L.; visualization, R.K.; supervision, K.R.K. and T.Z.; project administration, T.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Chinese Second Tibetan Plateau Scientific Expedition and Research (STEP) program (2019QZKK0105), the Natural Science Foundation of China (42030612), Fund Project of Central Asian Atmospheric Science Research (CAAS201913), and the Start-Up Research Grant (SRG) (File No. SRG/2020/001445) sanctioned by the Science and Engineering Research Board (SERB), a statutory body under the Department of Science and Technology (DST), Govt. of India. One of the authors, K.R.K., is grateful to the DST, Govt. of India, for the award of the DST-FIST Level-1 (SR/FST/PS-1/2018/35) scheme to the Department of Physics, KLEF.

Informed Consent Statement

Informed consent was obtained from all subjects involved in the study.

Data Availability Statement

The data that support the findings of this study are available open and free to the public and can be downloaded at https://gmao.gsfc.nasa.gov/reanalysis/MERRA-2/ (accessed on 1 December 2020) and https://modis.gsfc.nasa.gov/ (accessed on 1 December 2020).

Acknowledgments

We owe our sincere thanks to the Global Modeling and Assimilation Office (GMAO) for the MERRA reanalysis and the GES DISC for the dissemination of MERRA-2 data. Thanks are also due to the MODIS science team for providing and processing various satellite data sets. The authors would like to acknowledge the Editor of the Journal, and the four anonymous reviewers for their helpful comments and constructive suggestions towards the improvement of earlier versions of the manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Intergovernmental Panel on Climate Change (IPCC) 5th Assessment Report (AR5). Observed Climate Change Impacts Database, version 2.01; NASA Socioeconomic Data and Applications Center: Palisades, NY, USA, 2017. [Google Scholar]
  2. Charlson, R.J.; Schwartz, S.E.; Hales, J.M.; Cess, R.D.; Coakley, J.A.; Hansen, J.E.; Hofmann, D.J. Climate forcing by anthropogenic aerosols. Science 1992, 255, 423–430. [Google Scholar] [CrossRef]
  3. Dubovik, O.; Holben, B.N.; Eck, T.F.; Smirnov, A.; Kaufman, Y.J.; King, M.D.; Tanre, D.; Slutsker, I. Variability of absorption and optical properties of key aerosol types observed in worldwide locations. J. Atmos. Sci. 2002, 59, 590–608. [Google Scholar] [CrossRef]
  4. He, Q.Q.; Zhang, M.; Huang, B. Spatio-temporal variation and impact factors analysis of satellite-based aerosol optical depth over China from 2002 to 2015. Atmos. Environ. 2016, 129, 79–90. [Google Scholar] [CrossRef]
  5. Dockery, D.; Pope III, A.; Xu, X.; Spengler, J.; Ware, J.; Fay, M.; Ferris, B.; Speizer, F. An association between air pollution and mortality in six US cities. N. Engl. J. Med. 1993, 329, 1753–1759. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Liu, Q.; Wang, S.; Zhang, W.; Li, J.; Dong, G. The effect of natural and anthropogenic factors on PM2. 5: Empirical evidence from Chinese cities with different income levels. Sci. Total Environ. 2019, 653, 157–167. [Google Scholar] [CrossRef]
  7. Kaskaoutis, D.G.; Singh, R.P.; Gautam, R.; Sharma, M.; Kosmopoulus, P.G.; Tripathi, S.N. Variability and trends of aerosol properties over Kanpur, northern India using AERONET data (2001–2010). Environ. Res. Lett. 2012, 7, 024003. [Google Scholar] [CrossRef] [Green Version]
  8. Kumar, K.R.; Sivakumar, V.; Yin, Y.; Reddy, R.R.; Kang, N.; Diao, Y.; Adesina, A.J.; Yu, X. Long-term (2003–2013) climatological trends and variations in aerosol optical parameters retrieved from MODIS over three stations in South Africa. Atmos. Environ. 2014, 95, 400–408. [Google Scholar] [CrossRef]
  9. Kumar, M.; Parmar, K.S.; Kumar, D.B.; Mhawish, A.; Broday, D.M.; Mall, R.K.; Banerjee, T. Long-term aerosol climatology over Indo-Gangetic Plain: Trend, prediction and potential source fields. Atmos. Environ. 2018, 180, 37–50. [Google Scholar] [CrossRef]
  10. Boiyo, R.K.; Raghavendra, K.; Zhaoa, T. Spatial variations and trends in AOD climatology over East Africa during 2002–2016: A comparative study using three satellite data sets. Int. J. Climatol. 2018, 38, 1221–1240. [Google Scholar] [CrossRef]
  11. Sun, E.; Xu, X.; Che, H.; Tang, Z.; Gui, K.; An, L.; Lu, C.; Shi, G. Variation in MERRA-2 aerosol optical depth and absorption aerosol optical depth over China from 1980 to 2017. J. Atmos. Solar-Terrest. Phys. 2019, 186, 8–19. [Google Scholar] [CrossRef]
  12. Dey, S.; Di Girolamo, L. A decade of change in aerosol properties over the Indian subcontinent. Geophys. Res. Lett. 2011, 38, 14811. [Google Scholar] [CrossRef] [Green Version]
  13. Pan, H.; Wang, M.; Kumar, K.R.; Zhang, J.; Meng, L. A Decadal Global Climatology of Ice Cloud Fraction with Their Microphysical and Optical Properties Inferred from the CALIPSO and Reanalysis Data. Remote. Sens. 2020, 12, 3795. [Google Scholar] [CrossRef]
  14. Ali, G.; Bao, Y.; Ullah, W.; Ullah, S.; Guan, Q.; Liu, X.; Li, L.; Lei, Y.; Li, G.; Ma, J. Spatiotemporal Trends of Aerosols over Urban Regions in Pakistan and Their Possible Links to Meteorological Parameters. Atmosphere 2020, 11, 306. [Google Scholar] [CrossRef] [Green Version]
  15. Alam, K.; Qureshi, S.; Blaschke, T. Monitoring spatio-temporal aerosol patterns over Pakistan based on MODIS, TOMS and MISR satellite data and a HYSPLIT model. Atmos. Environ. 2011, 45, 4641–4651. [Google Scholar] [CrossRef]
  16. Alam, K.; Trautmann, T.; Blaschke, T.; Subhan, F. Changes in aerosol optical properties due to dust storms in the Middle East and Southwest Asia. Remote. Sens. Environ. 2014, 143, 216–227. [Google Scholar] [CrossRef]
  17. Alam, K.; Trautmann, T.; Blaschke, T.; Majid, H. Aerosol optical and radiative properties during summer and winter seasons over Lahore and Karachi. Atmos. Environ. 2012, 50, 234–245. [Google Scholar] [CrossRef]
  18. Levy, R.C.; Mattoo, S.; Munchak, L.A.; Remer, L.A.; Sayer, A.M.; Patadia, F.; Hsu, N.C. The Collection 6 MODIS aerosol products over land and ocean. Atmos. Meas. Tech. 2013, 6, 2989–3034. [Google Scholar] [CrossRef] [Green Version]
  19. Luo, Y.; Zheng, X.; Zhao, T.; Chen, J. A climatology of aerosol optical depth over China from recent 10 years of MODIS remote sensing data. Int. J. Clim. 2013, 34, 863–870. [Google Scholar] [CrossRef]
  20. Popp, T.; de Leeuw, C.G.; Bingen, C.; Brühl, V.; Capelle, A.; Chedin, L.; Clarisse, O.; Dubovik, R.; Grainger, J.; Griesfeller, A.; et al. Development, production and evaluation of aerosol Climate Data Records from European satellite observations (Aerosol_cci). Remote. Sens. 2016, 8, 421. [Google Scholar] [CrossRef] [Green Version]
  21. De Leeuw, G.; Sogacheva, L.; Rodriguez, E.; Kourtidis, K.; Georgoulias, A.K.; Alexandri, G.; Amiridis, V.; Proestakis, E.; Marinou, E.; Xue, Y.; et al. Two decades of satellite observations of AOD over mainland China using ATSR-2, AATSR and MODIS/Terra: Data set evaluation and large-scale patterns. Atmos. Chem. Phys. Discuss. 2018, 18, 1573–1592. [Google Scholar] [CrossRef] [Green Version]
  22. Sogacheva, L.; Popp, T.; Sayer, A.M.; Dubovik, O.; Garay, M.J.; Heckel, A.; Hsu, N.C.; Jethva, H.; Kahn, R.A.; Kolmonen, P.; et al. Merging regional and global aerosol optical depth records from major available satellite products. Atmos. Chem. Phys. Discuss. 2020, 20, 2031–2056. [Google Scholar] [CrossRef] [Green Version]
  23. Sogacheva, L.; de Leeuw, G.; Rodriguez, E.; Kolmonen, P.; Georgoulias, A.K.; Alexandri, G.; Kourtidis, K.; Proestakis, E.; Marinou, E.; Amiridis, V.; et al. Spatial and seasonal variations of aerosols over China from two decades of multi-satellite observations—Part 1: ATSR (1995–2011) and MODIS C6.1 (2000–2017). Atmos. Chem. Phys. 2018, 18, 11389–11407. [Google Scholar] [CrossRef] [Green Version]
  24. Torres, O.; Jethva, H.; Bhartia, P.K. Retrieval of Aerosol Optical Depth above Clouds from OMI Observations: Sensitivity Analysis and Case Studies. J. Atmos. Sci. 2012, 69, 1037–1053. [Google Scholar] [CrossRef]
  25. Mc Carty, W.; Coy, L.; Gelaro, R.; Huang, A.; Merkova, D.; Smith, E.B.; Sienkiewicz, M.; Wargan, K. MERRA-2input observations: Summary and assessment; NASA TM-2016- 104606; NASA Global Modeling and Assimilation Office: Phoenix, AZ, USA, 2016; Volume 46, p. 64. [Google Scholar]
  26. Randles, C.A.; Da Silva, A.M.; Buchard, V.; Colarco, P.R.; Darmenov, A.; Govindaraju, R.; Mirnov, A.; Holben, B.; Ferrare, R.; Hair, J.; et al. The MERRA-2 aerosol reanalysis, 1980 onward. Part I, System description and data assimilation evaluation. J. Clim. 2017, 30, 6823–6850. [Google Scholar] [CrossRef]
  27. Klingmüller, K.; Pozzer, A.; Metzger, S.; Stenchikov, G.L.; Lelieveld, J. Aerosol optical depth trend over the Middle East. Atmos. Chem. Phys. Discuss. 2016, 16, 5063–5073. [Google Scholar] [CrossRef] [Green Version]
  28. Che, H.; Gui, K.; Xia, X.; Wang, Y.; Holben, B.N.; Goloub, P.; Cuevas-Agulló, E.; Wang, H.; Zheng, Y.; Zhao, H.; et al. Large contribution of meteorological factors to inter-decadal changes in regional aerosol optical depth. Atmos. Chem. Phys. Discuss. 2019, 19, 10497–10523. [Google Scholar] [CrossRef] [Green Version]
  29. Khan, R.; Kumar, K.R.; Zhao, T. The climatology of aerosol optical thickness and radiative effects in Southeast Asia from 18-years of ground-based observations. Environ. Pollut. 2019, 254, 113025. [Google Scholar] [CrossRef] [PubMed]
  30. Gupta, P.; Khan, M.N.; Da Silva, A.; Patadia, F. MODIS aerosol optical depth observations over urban areas in Pakistan: Quantity and quality of the data for air quality monitoring. Atmos. Pollut. Res. 2013, 4, 43–52. [Google Scholar] [CrossRef] [Green Version]
  31. Wu, J.; Luo, J.; Zhang, L.; Xia, L.; Zhao, D.; Tang, J. Improvement of aerosol optical depth retrieval using visibility data in China during the past 50 years. J. Geophys. Res. Atmos. 2014, 119, 370. [Google Scholar] [CrossRef]
  32. Qin, W.; Liu, Y.; Wang, L.; Lin, A.; Xia, X.; Che, H.; Bilal, M.; Zhang, M. Characteristic and Driving Factors of Aerosol Optical Depth over Mainland China during 1980–2017. Remote. Sens. 2018, 10, 1064. [Google Scholar] [CrossRef] [Green Version]
  33. He, L.; Wang, L.; Lin, A.; Zhang, M.; Xia, X.; Tao, M.; Zhou, H. What drives changes in aerosol properties over the Yangtze River Basin in past four decade? Atmos. Environ. 2018, 190, 269–283. [Google Scholar] [CrossRef]
  34. Khan, R.; Kumar, K.R.; Zhao, T.; Ali, G. The contribution of different aerosol types to direct radiative forcing over distinct environments of Pakistan inferred from the AERONET data. Environ. Res. Lett. 2020, 15, 114062. [Google Scholar] [CrossRef]
  35. Bibi, H.; Alam, K.; Chishtie, F.; Bibi, S.; Shahid, I.; Blaschke, T. Intercomparison of MODIS, MISR, OMI, and CALIPSO Aerosol Optical Depth Retrievals for Four Locations on the Indo-Gangetic Plains and Validation against AERONET Data. Atmos. Environt. 2015, 111, 113–126. [Google Scholar] [CrossRef]
  36. Iftikhar, M.; Alam, K.; Sorooshian, A.; Syed, W.A.; Bibi, S.; Bibi, H. Contrasting aerosol optical and radiative properties between dust and urban haze episodes in megacities of Pakistan. Atmos. Environ. 2018, 173, 157–172. [Google Scholar] [CrossRef]
  37. Hsu, N.C.; Tsay, S.-C.; King, M.D.; Herman, J.R. Aerosol Properties over Bright-Reflecting Source Regions. IEEE Trans. Geosci. Remote. Sens. 2004, 42, 557–569. [Google Scholar] [CrossRef]
  38. Hsu, N.C.; Jeong, M.-J.; Bettenhausen, C.; Sayer, A.; Hansell, A.R.; Seftor, C.S.; Huang, J.; Tsay, S.-C. Enhanced Deep Blue aerosol retrieval algorithm: The second generation. J. Geophys. Res. Atmos. 2013, 118, 9296–9315. [Google Scholar] [CrossRef]
  39. Sayer, A.M.; Hsu, N.C.; Bettenhausen, C.; Jeong, M.J. Validation and uncertainty estimates for MODIS Collection 6 “deep Blue” aerosol data. J. Geophys. Res.-Atmos. 2013, 118, 7864–7872. [Google Scholar] [CrossRef] [Green Version]
  40. Sayer, A.; Munchak, L.A.; Hsu, N.C.; Levy, R.C.; Bettenhausen, C.; Jeong, M.J. MODIS Collection 6 aerosol products: Comparison between Aqua’s e-Deep Blue, Dark Target, and “merged” data sets, and usage recommendations. J. Geophys. Res. Atmos. 2014, 119. [Google Scholar] [CrossRef]
  41. Sayer, A.M.; Hsu, N.C.; Bettenhausen, C.; Jeong, M.; Meister, G. Effect of MODIS Terra radiometric calibration improvements on Collection 6 Deep Blue aerosol products: Validation and Terra/Aqua consistency. J. Geophys. Res. Atmos. 2015, 120. [Google Scholar] [CrossRef] [Green Version]
  42. Khalid, B.; Khalid, A.; Muslim, S.; Habib, A.; Khan, K.; Alvim, D.S.; Shakoor, S.; Mustafa, S.; Zaheer, S.; Zoon, M.; et al. Estimation of aerosol optical depth in relation to meteorological parameters over eastern and western routes of China Pakistan economic corridor. J. Environ. Sci. 2021, 99, 28–39. [Google Scholar] [CrossRef] [PubMed]
  43. Ilyas, S.Z.; Khattak, A.I.; Nasir, S.M.; Qurashi, T.; Durrani, R. Air pollution assessment in urban areas and its impact on human health in the city of Quetta, Pakistan. Clean Technol. Environ. Policy 2009, 12, 291–299. [Google Scholar] [CrossRef]
  44. Gelaro, R.; Mccarty, W.; Suárez, M.J.; Todling, R.; Molod, A.; Takacs, L.; Randles, C.A.; Darmenov, A.; Bosilovich, M.G.; Reichle, R.; et al. The Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2). J. Clim. 2017, 30, 5419–5454. [Google Scholar] [CrossRef]
  45. Buchard, V.; Randles, C.A.; Da Silva, A.M.; Darmenov, A.; Colarco, P.R.; Govindaraju, R.; Ferrare, R.; Hair, J.; Beyersdorf, A.J.; Ziemba, L.D.; et al. The MERRA-2 Aerosol Reanalysis, 1980 Onward. Part II: Evaluation and Case Studies. J. Clim. 2017, 30, 6851–6872. [Google Scholar] [CrossRef] [PubMed]
  46. Feng, F.; Wang, K. Does the modern-era retrospective analysis for research and applications-2 aerosol reanalysis introduce an improvement in the simulation of surface solar radiation over China? Int. J. Clim. 2018, 39, 1305–1318. [Google Scholar] [CrossRef]
  47. Kaufman, Y.J.; Tanré, D.; Remer, L.A.; Vermote, E.F.; Chu, A.; Holben, B.N. Operational remote sensing of tropospheric aerosol over land from EOS moderate resolution imaging spectroradiometer. J. Geophys. Res. Space Phys. 1997, 102, 17051–17067. [Google Scholar] [CrossRef]
  48. Mann, H.B. Nonparametric tests against trend. Econometrica 1945, 13, 245–259. [Google Scholar] [CrossRef]
  49. Sonali, P.; Kumar, N.D. Review of trend detection methods and their application to detect temperature changes in India. J. Hydrol. 2013, 476, 212–227. [Google Scholar] [CrossRef]
  50. Hirsch, R.M.; Slack, J.R. A Nonparametric Trend Test for Seasonal Data with Serial Dependence. Water Resour. Res. 1984, 20, 727–732. [Google Scholar] [CrossRef] [Green Version]
  51. Salman, S.A.; Shahid, S.; Ismail, T.; Chung, E.-S.; Al-Abadi, A.M. Long-term trends in daily temperature extremes in Iraq. Atmos. Res. 2017, 198, 97–107. [Google Scholar] [CrossRef]
  52. Hamed, K.H.; Ramachandra, R. A modified Mann-Kendall trend test for auto correlated data. J. Hydrol. 1998, 204, 182–196. [Google Scholar] [CrossRef]
  53. Partal, T. Wavelet transform-based analysis of periodicities and trends of Sakarya basin (Turkey) streamflow data. River Res. Appl. 2009, 26, 695–711. [Google Scholar] [CrossRef]
  54. Nasri, M.; Modarres, R. Dry spell trend analysis of Isfahan Province, Iran. Int. J. Clim. 2009, 29, 1430–1438. [Google Scholar] [CrossRef]
  55. Ahmad, I.; Zhang, F.; Tayyab, M.; Anjum, M.N.; Zaman, M.; Liu, J.; Farid, H.U.; Saddique, Q. Spatiotemporal analysis of precipitation variability in annual, seasonal and extreme values over upper Indus River basin. Atmos. Res. 2018, 213, 346–360. [Google Scholar] [CrossRef]
  56. Rashid, M.; Beecham, S.; Chowdhury, R.K. Assessment of trends in point rainfall using Continuous Wavelet Transforms. Adv. Water Resour. 2015, 82, 1–15. [Google Scholar] [CrossRef]
  57. Tsidu, M.G.; Bililign, S. Comparison of MERRA-2 and Collection 6 Modis AODs at 0.55 Micron over Diverse Surfaces of Southern Africa during 2000-2017. In Proceedings of the AGU Fall Meeting, Washington, DC, USA, 10–14 December 2018. 51J-2.2301. [Google Scholar]
  58. Ahmed, K.; Shahid, S.; Chung, E.; Ismail, T.; Wang, X. Spatial distribution of secular trends in annual and seasonal precipitation over Pakistan. Clim. Res. 2017, 74, 95–107. [Google Scholar] [CrossRef]
  59. Von Storch, H. Misuses of statistical analysis in climate research. In Analysis of Climate Variability; Springer: Berlin/Heidelberg, Germany, 1999; pp. 11–26. [Google Scholar]
  60. Douglas, E.; Vogel, R.; Kroll, C. Trends in floods and low flows in the United States: Impact of spatial correlation. J. Hydrol. 2000, 240, 90–105. [Google Scholar] [CrossRef]
  61. Rupakheti, D.; Kang, S.; Rupakheti, M.; Cong, Z.; Panday, A.K.; Holben, B.N. Identification of absorbing aerosol types at a site in the northern edge of Indo-Gangetic plain and a polluted valley in the foothills of the central Himalayas. Atmos. Res. 2019, 223, 15–23. [Google Scholar] [CrossRef]
  62. Mehta, M.; Singh, R.; Singh, A.; Singh, N. Anshumali. Recent global aerosol optical depth variations and trends—A comparative study using MODIS and MISR level 3 datasets. Remote. Sens. Environ. 2016, 181, 137–150. [Google Scholar] [CrossRef]
  63. Wang, J.; Wu, S.; Zhou, S.L.; Tang, J. Relationship between aerosol optical depth and population density in typical area of East China. Appl. Ecol. Env. Res. 2017, 15, 1041–1056. [Google Scholar] [CrossRef]
Figure 1. (a) Topographic map of the study area, Pakistan, with the elevation represented with the color scale at the bottom left. The areas enclosed in red boxes represent major regions in the study domain, while the locations marked with blue solid circles denote the major cities in the subdomain of the study region. (b) A geographical map in the right panel showing the location of Pakistan in South Asia with its bordering countries, including potential source desert regions, the direction of air masses towards the study region is marked with an arrow. Different colors differentiate the desert area from the normal land cover type.
Figure 1. (a) Topographic map of the study area, Pakistan, with the elevation represented with the color scale at the bottom left. The areas enclosed in red boxes represent major regions in the study domain, while the locations marked with blue solid circles denote the major cities in the subdomain of the study region. (b) A geographical map in the right panel showing the location of Pakistan in South Asia with its bordering countries, including potential source desert regions, the direction of air masses towards the study region is marked with an arrow. Different colors differentiate the desert area from the normal land cover type.
Remotesensing 13 00822 g001
Figure 2. Scatter plot of MERRA-2 versus MODIS-Aqua AOD for (a) DB, (b) DT, and (c) DTB (merged) data product over the study domain Pakistan, for the period (2002–2018). The data points are collocated monthly mean data, indicated by the solid circles. The red line indicates a linear least-square fit the data points. The number of matched data points (N), (R) Pearson’s correlation coefficient, Bias, and root-mean-square error (RMSE) is also indicated in all the panels.
Figure 2. Scatter plot of MERRA-2 versus MODIS-Aqua AOD for (a) DB, (b) DT, and (c) DTB (merged) data product over the study domain Pakistan, for the period (2002–2018). The data points are collocated monthly mean data, indicated by the solid circles. The red line indicates a linear least-square fit the data points. The number of matched data points (N), (R) Pearson’s correlation coefficient, Bias, and root-mean-square error (RMSE) is also indicated in all the panels.
Remotesensing 13 00822 g002
Figure 3. Autocorrelation of annual aerosol optical depth (AOD550) (Modern-era retrospective analysis for research and applications, version 2 (MERRA-2)) time-series as a function of the lag, covering an annual mean data of 20 years period (1980–1999), evaluated for six main study regions (a) North Dry Region (NDR), (b) Katawaz Basin (c) Indus Basin (d) Balochistan Plateau (e) Coastal region and (f) Kharan Desert Region (KDR) covering the study domain of Pakistan.
Figure 3. Autocorrelation of annual aerosol optical depth (AOD550) (Modern-era retrospective analysis for research and applications, version 2 (MERRA-2)) time-series as a function of the lag, covering an annual mean data of 20 years period (1980–1999), evaluated for six main study regions (a) North Dry Region (NDR), (b) Katawaz Basin (c) Indus Basin (d) Balochistan Plateau (e) Coastal region and (f) Kharan Desert Region (KDR) covering the study domain of Pakistan.
Remotesensing 13 00822 g003
Figure 4. Annual and seasonal spatial distributions of AOD averaged for 39 years (1980–2018) retrieved from the MERRA-2 (ai) and Moderate-resolution imaging spectroradiometer (MODIS)-Aqua (DB) (bj; 2002–2018) over the study domain, Pakistan. The blank spaces in the figure are areas for which no MODIS data are available.
Figure 4. Annual and seasonal spatial distributions of AOD averaged for 39 years (1980–2018) retrieved from the MERRA-2 (ai) and Moderate-resolution imaging spectroradiometer (MODIS)-Aqua (DB) (bj; 2002–2018) over the study domain, Pakistan. The blank spaces in the figure are areas for which no MODIS data are available.
Remotesensing 13 00822 g004
Figure 5. Month-to-month variation of the AOD averaged over the years 2002–2018, for 12 cities in Pakistan (a) Gilgit, Muzafarabad (b) Peshawar, D.I.Khan (c) Lahore, Multan (d) Quetta, Sibbi (e) Chhor, Karachi (f) Khuzdar, Panjgur for both MERRA-2 (red-colored histograms) and MODIS-DB (blue colored histograms). Note that the yellow bars are plotted on top of the blue bars. The blank columns in figure(a) are areas for which no MODIS-DB data are available. Note that the scales are different for each figure.
Figure 5. Month-to-month variation of the AOD averaged over the years 2002–2018, for 12 cities in Pakistan (a) Gilgit, Muzafarabad (b) Peshawar, D.I.Khan (c) Lahore, Multan (d) Quetta, Sibbi (e) Chhor, Karachi (f) Khuzdar, Panjgur for both MERRA-2 (red-colored histograms) and MODIS-DB (blue colored histograms). Note that the yellow bars are plotted on top of the blue bars. The blank columns in figure(a) are areas for which no MODIS-DB data are available. Note that the scales are different for each figure.
Remotesensing 13 00822 g005
Figure 6. Temporal variation and linear trends in the annual mean AOD for the cities considered in this study for the years 1980-2018: the blue line shows the monthly mean AOD550, and the solid red line indicates the linear trend fitted to these data (Ya: 1980–2000; Yb: 2001–2018). The intercept and slope of the fit lines are shown at the top right of each figure. City names (abbreviations, see Table 2) are indicated in the top left of each figure.
Figure 6. Temporal variation and linear trends in the annual mean AOD for the cities considered in this study for the years 1980-2018: the blue line shows the monthly mean AOD550, and the solid red line indicates the linear trend fitted to these data (Ya: 1980–2000; Yb: 2001–2018). The intercept and slope of the fit lines are shown at the top right of each figure. City names (abbreviations, see Table 2) are indicated in the top left of each figure.
Remotesensing 13 00822 g006aRemotesensing 13 00822 g006b
Figure 7. Temporal variation and linear trends in monthly mean AOD, in four different seasons during 1980–2018, for the six main study domains in subfigures (a), (g), (m), (s) (NDR), (b), (h), (n), (t) (Katawaz Basin), (c), (i), (o), (u) (Indus Basin) (d), (j), (p), (v) (Balochistan Plateau), (e), (k), (q), (w) (Coastal region), and (f), (l), (r), (x) (KDR). The red lines show the monthly averaged AOD time-series, and the solid blue lines show the linear trend. The fitted trend lines are shown in the top right corner of each plot.
Figure 7. Temporal variation and linear trends in monthly mean AOD, in four different seasons during 1980–2018, for the six main study domains in subfigures (a), (g), (m), (s) (NDR), (b), (h), (n), (t) (Katawaz Basin), (c), (i), (o), (u) (Indus Basin) (d), (j), (p), (v) (Balochistan Plateau), (e), (k), (q), (w) (Coastal region), and (f), (l), (r), (x) (KDR). The red lines show the monthly averaged AOD time-series, and the solid blue lines show the linear trend. The fitted trend lines are shown in the top right corner of each plot.
Remotesensing 13 00822 g007aRemotesensing 13 00822 g007b
Figure 8. Spatial distribution of annual and seasonal mean AOD trend from MODIS-DB (d) 2002–2018 andMERRA-2 during (a) 1980–2018 and for the (b) pre-MODIS (1980–2000) and (c) post-MODIS (2002–2018) periods observed over the study domain, Pakistan and the surrounding region. The shaded color represents the sign and magnitude of trend, the confidence levels (95%) are shown as black dots in each of the subplots. The blank spaces in figure (d) are areas for which no MODIS data are available.
Figure 8. Spatial distribution of annual and seasonal mean AOD trend from MODIS-DB (d) 2002–2018 andMERRA-2 during (a) 1980–2018 and for the (b) pre-MODIS (1980–2000) and (c) post-MODIS (2002–2018) periods observed over the study domain, Pakistan and the surrounding region. The shaded color represents the sign and magnitude of trend, the confidence levels (95%) are shown as black dots in each of the subplots. The blank spaces in figure (d) are areas for which no MODIS data are available.
Remotesensing 13 00822 g008
Figure 9. Hovmoller (time-series versus longitude) diagrams of MERRA-2 retrieved AOD showing the movement of the columnar aerosol load and the associated gradient in time and longitude parameters for different regions (af) for the period 1980–2018.
Figure 9. Hovmoller (time-series versus longitude) diagrams of MERRA-2 retrieved AOD showing the movement of the columnar aerosol load and the associated gradient in time and longitude parameters for different regions (af) for the period 1980–2018.
Remotesensing 13 00822 g009
Figure 10. Variation of AODMERRA-2 with time (2000–2018) for different study regions with different elevations (indicated in the top of each figure) and with population density(unit/km2) on the secondary y-axis (a) NDR, (b) Katawaz Basin (c) Indus Basin (d) Balochistan Plateau (e) Coastal region and (f) KDR. The red data points and line indicate the AOD;green indicates population density.
Figure 10. Variation of AODMERRA-2 with time (2000–2018) for different study regions with different elevations (indicated in the top of each figure) and with population density(unit/km2) on the secondary y-axis (a) NDR, (b) Katawaz Basin (c) Indus Basin (d) Balochistan Plateau (e) Coastal region and (f) KDR. The red data points and line indicate the AOD;green indicates population density.
Remotesensing 13 00822 g010
Figure 11. Probability density distribution functions (PDFs) of monthly mean AOD for different periods (11 years for 1980–1990, 8 years for 2011–2018, 10 years for 1991–2000 and 2001–2010) during the study period (1980–2018) (see legend in Figure (a), for different study regions,(a) NDR, (b)Katawaz Basin (c) Indus Basin (d) Balochistan Plateau (e) coastal region and (f) KDR.
Figure 11. Probability density distribution functions (PDFs) of monthly mean AOD for different periods (11 years for 1980–1990, 8 years for 2011–2018, 10 years for 1991–2000 and 2001–2010) during the study period (1980–2018) (see legend in Figure (a), for different study regions,(a) NDR, (b)Katawaz Basin (c) Indus Basin (d) Balochistan Plateau (e) coastal region and (f) KDR.
Remotesensing 13 00822 g011
Table 1. The trend for the different study regions around the globe, taken from the previous studies in comparison of our results for the study domain of Pakistan. The long-term annual trends of AOD by different data sets are presented for different periods.
Table 1. The trend for the different study regions around the globe, taken from the previous studies in comparison of our results for the study domain of Pakistan. The long-term annual trends of AOD by different data sets are presented for different periods.
RegionTrend (Year−1)PeriodData SetCitation
Pakistan(Indus Basin)0.0081980–2018MERRA-2Present study
Pakistan(Balochistan)0.00521980–2018MERRA-2Present study
Pakistan(Kharan)0.0021980–2018MERRA-2Present study
Pakistan(Karachi)0.012001–2018AERONETKhan et al. [29]
Pakistan(Lahore)−0.0022001–2018AERONETKhan et al. [29]
Pakistan(Multan)0.0022006–2015MODISKumar et al. [9]
Pakistan(Lahore)0.0032006–2015MODISKumar et al. [9]
Pakistan(Karachi)0.0042006–2015MODISKumar et al. [9]
Pakistan(Karachi)−0.00132003–2017MODISAli et al. [14]
Pakistan(Quetta)−0.00252003–2017MODISAli et al. [14]
Pakistan(Lahore)0.00242003–2017MODISAli et al. [14]
Pakistan(Lahore)−0.00062001–2010MODISGupta et al. [30]
Middle East (Iran)0.0052001–2012MODISKlingmülleret al. [27]
Middle East (Iraq)0.0342001–2013MODISKlingmülleret al. [27]
China0.0041960–2009MISRWu et al. [31]
China0.0021980–2017MERRA-2Qin et al. [32]
East Africa0.0022009–2016MODIS-DBBoiyoet al. [10]
China(Yangtze River Delta)0.00651980–2016MERRA-2He et al. [33]
Table 2. Information on selected sites located in the distinct subdomains of Pakistan, and their climate and environment, i.e., geographical location (coordinates: latitude, longitude), elevation, population density (ρ /km2), and precipitation (mm), for the year 2018.
Table 2. Information on selected sites located in the distinct subdomains of Pakistan, and their climate and environment, i.e., geographical location (coordinates: latitude, longitude), elevation, population density (ρ /km2), and precipitation (mm), for the year 2018.
SiteLat (N°)Lon (E°)Elevation (m)ρPrecipitation (mm)Study RegionClimate
(per km2)
Gilgit (Gil)35.9274.33146017159Northern Dry RegionSemiarid region
Muzaffarabad (Mzf)34.3773.488383001457
Peshawar (Psh)34.0271.58328287384Katawaz BasinSemiurban
D.I.Khan (DIK)31.4970.56171220.4268.8
Lahore (Lhr)31.3574.242146300628.8Indus BasinIndustrialized, desert
Multan (Mlt)30.1271.261226500186.8
Quetta (Qt)30.1167.5717195600260.09Balochistan PlateauSuburban, hot desert
Sibbi (Sb)29.3367.5513332144.3
Chhor (Cr)25.3169.42543216.03Coastal RegionUrban, coastal
Karachi (Kr)24.8967.16223900174.7
Khuzdar (Kz)27.566.38123123252.4Kharan Desert RegionSemiurban
Panjgur (Pj)26.5864.196819108.7
Table 3. Annually and seasonally averaged MERRA-2 AOD (±SD) over the years 1980–2018 for thesub-domains in Pakistan.
Table 3. Annually and seasonally averaged MERRA-2 AOD (±SD) over the years 1980–2018 for thesub-domains in Pakistan.
SchemeWinterSpringSummerAutumnAnnual
NDR0.08 ± 0.050.15 ± 0.040.14 ± 0.040.08 ± 0.050.12 ± 0.04
Katawaz Basin0.29 ± 0.170.46 ± 0.120.52 ± 0.160.34 ± 0.170.40 ± 0.15
Indus Basin0.23 ± 0.120.37 ± 0.080.45 ± 0.120.28 ± 0.130.33 ± 0.09
Balochistan Plateau0.14 ± 0.060.30 ± 0.060.35 ± 0.080.19 ± 0.080.24 ± 0.07
Coastal Region0.28 ± 0.160.48 ± 0.110.57 ± 0.150.34 ± 0.170.41 ± 0.13
KDR0.25 ± 0.140.44 ± 0.090.54 ± 0.130.31 ± 0.140.38 ± 0.12
Table 4. Trends (%) in the variations of the seasonally and annually-averaged MERRA-2 AOD during 1980–2018, for selected regions, with the statistical significance (p-value) for the annual trends.
Table 4. Trends (%) in the variations of the seasonally and annually-averaged MERRA-2 AOD during 1980–2018, for selected regions, with the statistical significance (p-value) for the annual trends.
Study RegionsWinter (%)Spring (%)Summer (%)Autumn (%)Annual (%)p-Value
NDR0.090.190.330.540.295210 E-9
Katawaz Basin3.061.610.993.512.290.0001 E-9
Indus Basin2.870.391.033.171.870.001 E-9
Balochistan Plateau3.351.230.653.232.120.0001 E-9
Coastal Region0.580.360.190.820.490.0001 E-9
KDR2.421.090.52.211.560.0001 E-9
Pakistan2.060.810.622.251.430.141 E-0
Table 5. Interdecadal changes in the seasonal mean (±SD) AOD obtained from the MERRA-2 model for the distinct subdomains and entire study domain during 1980–2018.
Table 5. Interdecadal changes in the seasonal mean (±SD) AOD obtained from the MERRA-2 model for the distinct subdomains and entire study domain during 1980–2018.
SeasonPeriodNDRKatawaz BasinIndus
Basin
Balochistan PlateauCoastal RegionKDRPakistan
Winter1980–19900.06 ± 0.040.13 ± 0.030.12 ± 0.040.12 ± 0.030.14 ± 0.030.14 ± 0.040.12 ± 0.03
1991–20000.08 ± 0.080.17 ± 0.080.16 ± 0.080.14 ± 0.070.17 ± 0.070.16 ± 0.060.15 ± 0.07
2001–20100.06 ± 0.130.41 ± 0.090.30 ± 0.070.13 ± 0.020.40 ± 0.080.31 ± 0.070.27 ± 0.05
2011–20180.08 ± 0.020.48 ± 0.090.36 ± 0.100.17 ± 0.090.46 ± 0.110.411 ± 0.140.33 ± 0.09
Spring1980–19900.13 ± 0.040.37 ± 0.130.30 ± 0.080.27 ± 0.050.42 ± 0.120.38 ± 0.100.31 ± 0.08
1991–20000.14 ± 0.040.40 ± 0.130.32 ± 0.090.27 ± 0.040.42 ± 0.110.37 ± 0.080.32 ± 0.07
2001–20100.16 ± 0.020.55 ± 0.120.44 ± 0.090.33 ± 0.040.56 ± 0.130.49 ± 0.110.42 ± 0.08
2011–20180.15 ± 0.040.51 ± 0.120.40 ± 0.090.32 ± 0.070.53 ± 0.130.50 ± 0.100.40 ± 0.08
Summer1980–19900.13 ± 0.050.43 ± 0.120.37 ± 0.070.33 ± 0.080.50 ± 0.120.48 ± 0.110.37 ± 0.08
1991–20000.14 ± 0.050.47 ± 0.120.39 ± 0.090.30 ± 0.070.51 ± 0.110.49 ± 0.100.38 ± 0.08
2001–20100.14 ± 0.030.61 ± 0.170.54 ± 0.100.37 ± 0.060.65 ± 0.160.61 ± 0.120.49 ± 0.09
2011–20180.16 ± 0.030.58 ± 0.130.51 ± 0.090.36 ± 0.070.62 ± 0.130.58 ± 0.110.47 ± 0.08
Autumn1980–19900.07 ± 0.040.19 ± 0.070.17 ± 0.060.17 ± 0.060.21 ± 0.080.21 ± 0.090.17 ± 0.06
1991–20000.09 ± 0.080.25 ± 0.110.21 ± 0.090.19 ± 0.080.24 ± 0.090.24 ± 0.100.20 ± 0.09
2001–20100.07 ± 0.010.46 ± 0.090.36 ± 0.050.18 ± 0.040.45 ± 0.080.38 ± 0.060.32 ± 0.04
2011–20180.10 ± 0.030.51 ± 0.130.41 ± 0.110.24 ± 0.080.49 ± 0.170.46 ± 0.150.37 ± 0.10
Table 6. Trends for the seasonal and annual statistical parameters p, z, and c, over the selected regions during 1980–2018 based on MERRA-2 AOD product. “Bold and italic” numbers represent trends that are statistically significant at 95% (p < 0.05). Here “c” indicates the overall trend (1980–2018) observed in the study region.
Table 6. Trends for the seasonal and annual statistical parameters p, z, and c, over the selected regions during 1980–2018 based on MERRA-2 AOD product. “Bold and italic” numbers represent trends that are statistically significant at 95% (p < 0.05). Here “c” indicates the overall trend (1980–2018) observed in the study region.
Study RegionStatisticsWinterSpringSummerAutumnAnnual
p0.0330.0077.5 E-074.51 E-075.21E-07
NDRz2.122.695.373.55.0185
c0.0070 E-020.026 E-020.048E-020.046 E-020.031E-02
p4.13E-091.48E-076.41E-041.05E-120.00E+00
Katawaz Basinz5.87915.25463.41387.123714.683
c0.00840.00630.00500.01210.0351
p4.22E-158.20E-091.39E-051.33E-071.45E-13
Indus Basinz7.8465.76444.34525.27427.3913
c0.00650.00120.00460.00890.0053
p2.86E-134.05E-060.00166.06E-060.00E+00
Balochistan Plateauz7.34.60893.15744.520
c0.0090.00480.00370.0110.0071
p0.0090.25840.50660.01531.77E-04
Coastal Regionz2.61291.13020.66412.42433.7498
c0.087E-020.0010.064E-020.00160.001
p0.00210.17990.5840.00390.004
KDRz3.07261.3410.54762.88572.8791
c0.00580.00410.002710.00700.0045
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Khan, R.; Kumar, K.R.; Zhao, T.; Ullah, W.; de Leeuw, G. Interdecadal Changes in Aerosol Optical Depth over Pakistan Based on the MERRA-2 Reanalysis Data during 1980–2018. Remote Sens. 2021, 13, 822. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13040822

AMA Style

Khan R, Kumar KR, Zhao T, Ullah W, de Leeuw G. Interdecadal Changes in Aerosol Optical Depth over Pakistan Based on the MERRA-2 Reanalysis Data during 1980–2018. Remote Sensing. 2021; 13(4):822. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13040822

Chicago/Turabian Style

Khan, Rehana, Kanike Raghavendra Kumar, Tianliang Zhao, Waheed Ullah, and Gerrit de Leeuw. 2021. "Interdecadal Changes in Aerosol Optical Depth over Pakistan Based on the MERRA-2 Reanalysis Data during 1980–2018" Remote Sensing 13, no. 4: 822. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13040822

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop