Impacts of Aerosol Loading in the Hindu Kush Himalayan Region Based on MERRA-2 Reanalysis Data

Abstract

The impacts of climate change have severely affected geosphere, biosphere and cryosphere ecosystems in the Hindu Kush Himalayan (HKH) region. The impact has been accelerating further during the last few decades due to rapid increase in anthropogenic activities such as modernization, industrialization and urbanization, along with energy demands. In view of this, the present work attempts to examine aerosol optical depth (AOD) over the HKH region using the long-term homogeneous MERRA-2 reanalysis data from January, 1980 to December, 2020. The AOD trends are examined statistically with student’s t-test (t). Due to a vast landmass, fragile topography and harsh climatic conditions, we categorized the HKH region into three sub-regions, namely, the northwestern and Karakoram (HKH1), the Central (HKH2) and the southeastern Himalaya and Tibetan Plateau (HKH3). Among the sub-regions, the significant enhancement of AOD is observed at several potential sites in the HKH2 region, namely, Pokhara, Nainital, Shimla and Dehradun by 55.75 × 10${}^{-4}$ ± 3.76 × 10${}^{-4}$, 53.15 × 10${}^{-4}$ ± 3.94 × 10${}^{-4}$, 51.53 × 10${}^{-4}$ ± 4.99 × 10${}^{-4}$ and 39.16 × 10${}^{-4}$ ± 4.08 × 10${}^{-4}$ AOD year${}^{-1}$ (550 nm), respectively, with correlation coefficients (Rs) of 0.86 to 0.93. However, at a sub-regional scale, HKH1, HKH2 and HKH3 exhibit 23.33 × 10${}^{-4}$ ± 2.28 × 10${}^{-4}$, 32.20 × 10${}^{-4}$ ± 2.58 × 10${}^{-4}$ and 9.48 × 10${}^{-4}$ ± 1.21 × 10${}^{-4}$ AOD year${}^{-1}$, respectively. The estimated trends are statistically significant (t > 7.0) with R from 0.81 to 0.91. Seasonally, the present study also shows strong positive AOD trends at several potential sites located in the HKH2 region, such as Pokhara, Nainital, Shimla and Dehradun, with minimum 19.81 × 10${}^{-4}$ ± 3.38 × 10${}^{-4}$ to maximum 72.95 × 10${}^{-4}$ ± 4.89 × 10${}^{-4}$ AOD year${}^{-1}$ with statistical significance. In addition, there are also increasing AOD trends at all the high-altitude background sites in all seasons.

1. Introduction

The Hindu Kush Himalayan (HKH) region is highly sensitive to global warming and its impacts [1,2]. Due to rapidly increasing surface temperatures over the HKH (also well known as the World’s “Third Pole”) region, the melting of glaciers is occurring much faster than in other parts of the Earth and is approaching an alarming rate [1,3,4]. The enhancement of dust particles with carbonaceous aerosols such as black carbon (BC) over the HKH region has been one of the main concerns during the last few decades [5,6,7,8]. Due to the consequences of their impacts, a great threat is emerging for 240 million people living in the HKH region along with several endangered species of wild fauna and flora as reported in the Hindu Kush Himalaya Assessment report [9]. The enhancement of BC aerosol could be responsible for the melting of glaciers as it causes a decrease in snow surface albedo by depositing over snow and darkening its surface [3,10,11,12]. Further, Reference [13] reported BC and other light-absorbing impurities in snow in the Chilean Andes during Austral winters 2015 and 2016 at 22 sites between latitudes 18$°$ S and 41$°$ S. They noticed that light absorption by impurities in snow was dominated by dust rather than BC. In case of the Himalayan region, such study is required to examine the different characteristics of BC and ligh-absorbing impurities in detail. The BC aerosol absorbs the solar radiation and causes warming in the atmosphere [14] and also affects the distribution, lifetime and micro-physical properties of clouds through indirect and semi-direct effects [15]. Since the entire HKH region is situated in the close vicinity of Indo-Gangetic Plains (IGP), there is a significant contribution from anthropogenic aerosols undergoing long-range transportation through the prevailing winds over the high altitude pristine region (above 4000 m above mean sea level (amsl)) in the HKH region [16,17,18,19]. Sometimes, these transported aerosols can reach the Arabian Sea and central south India [20] as well as remote areas of the northwestern and central Himalayas [17,21,22,23,24].

The robust estimation of long-term climate variable parameters in the HKH region is inadequate due to sparse and discontinuous observations [1,25]. With the increasing emissions of anthropogenic greenhouse gases, the cryosphere processes coupled with the hydrological regimes are attributed to warming over the region during the last few decades [1]. Besides global warming, other regional atmospheric forcing parameters such as anthropogenic aerosols have influenced the South Asian summer monsoon [26]. Intensive monsoonal rainfall over the northern part of India and western Nepal during 2013 caused severe rainfall with landslides over the region. This caused the worst floods in HKH region and has been linked with the increasing greenhouse gases and atmospheric aerosols [27]. Apart from the increasing levels of greenhouse gases, a trend of increasing aerosols over the Indian sub-continent was reported by several studies in the past [28,29,30,31]. However, most of these studies are reported only from the urban or plain sites, and there is limited information about the high-altitude aged background sites located in the HKH region in particular. Due to vast landmass, different fragile topography and climatic conditions, the HKH entire region is divided into three sub-regions, namely, the northwestern Himalaya and Karakoram (HKH1: 71–79° E, 32–39$°$ N), the central Himalaya (HKH2: 76–93$°$ E, 27–32$°$ N) and the southeast Himalaya and Tibetan Plateau (HKH3: 93–103$°$ E, 28–36$°$ N), as reported by [4]. Among these sub-regions, the HKH-1 region is one of the most climate-sensitive regions in the sub-continent. Due to the complex topography and coarse resolution of global climate models, the common consensus among the models for the HKH region is very weak, as reported by [1]. Further, the availability of data has large inconsistencies as well as non-homogeneity, apparently due to lack of ground-based observation over the harsh climatic condition and mountain terrain [1,9,25]. Recently, [32] have reported multidecade trend of aerosol-radiative properties at 52 observatories worldwide, including high-altitude mountain site Mukteshwar (29.44$°$ N, 79.62$°$ E, 2180 m amsl), located in the HKH2 region under the Global Atmospheric Watch (GAW) program. However, the HKH2 region has more aerosol measuring stations than the HKH1 and HKH3 sub-regions. Despite the lack of ground-based observation and other challenges such as harsh climatic condition and mountain terrain, [33] found increasing Aerosol Optical depth (AOD) trends over the high-altitude remote background sites located in the HKH region using AErosol RObotic NETwork (AERONET; [34]) and Sky radiometer network (SKYNET; [35]) data.

With this background knowledge, in the current work, we examined the long-term aerosol optical depth (AOD) trends over the HKH region using 41 years (1980–2020) of homogeneous reanalysis data from MERRA-2 (Modern-Era Retrospective Analysis for Research and Applications, Version 2; [36]) as well as ground-based data, preferably from AERONET, SKYNET, SONET (Sun–Sky Radiometer Observation Network; [37]), etc., which have minimum 4–5 years of observation. The purpose of selecting such ground-based data (or stations) in the present work is to validate or compare with reanalysis data available from the early 1980s across the globe. The reanalysis data have more advantages than the satellite data as far as the long-term data are concerned. The satellite era data are available from the early 2000s. Such study of long-term continuous data for AOD trend analysis may be the first of this kind in the region. Further, the present study has many important aspects for identifying various climate forcing parameters in the HKH region, where the satellite-based measurements are inaccurate and biased due to extremely low AOD [38,39,40,41]. This paper is organized as follows: Section 2 describes observational sites with data sources; Section 3 describes methodology and data analysis; Section 4 presents results and discussion for yearly AOD trend, seasonal long-term AOD trends and climatological AOD in the HKH sub-regions; and Section 5 summarizes and discusses the implications of the present work.

2. Observational Sites

In order to study long-term AOD trends, it is desirable to use homogeneous and continuous long-term data within a range of more than 20–30 years. In the recent past, satellite data have become available from the early 2000s. However, such satellite data have several issues over the background mountainous sites when the retrieved parameters are of the order of uncertainties of the satellite data, as reported by several authors [38,39,40,41]. Due to limitations of such satellite data, the use of reanalysis data has several advantages due to availability of sufficiently long-term data (i.e, available from the early 1980s) without any discontinuity. Moreover, such reanalysis data are evaluated using observation (ground and space) and model data. In the present work, we used MERRA-2 data AOD at 550 nm [36] for the period 1980–2020. Further, the MERRA-2 data have had the bias corrected using several space-based remote-sensing platforms such as Moderate Resolution Imaging Spectroradiometer (MODIS), Advanced Very High-Resolution Radiometer (AVHRR) instruments [42], Multiangle Imaging Spectro Radiometer (MISR) [43], and ground-based direct measurements from AERONET [34,44]. Several further improvements have been made in MERRA-2 data products, including the hydrological cycle, the stratosphere, ozone and cryospheric process, as reported therein. Therefore, MERRA-2 data used in the present work are very comprehensive, including augmentation products from the previous version that consider all the relevant aspects.

To study the AOD trend over the entire HKH region, we also selected 16 stations due to different fragile topography/orography and climatic conditions within 1 × 1$°$ grid areas in addition to the wider areas of HKH1 (71–79$°$ E, 32–39$°$ N), HKH2 (76–93$°$ E, 27–32$°$ N) and HKH3 (93–103$°$ E, 28–36$°$ N) sub-regions. The main criteria for selection of these sites are as follows: (i) they must either be a major tourist site or be becoming a potential sources of anthropogenic emission, and (ii) they must be an aerosol measuring station (AERONET/SKYNET/SONET etc.) with a minimum of 4–5 years of continuous data. The purpose of selecting such ground-based data is to examine the MERRA-2 AOD data since satellite data have several uncertainties in high-altitude mountain regions in particular, as reported earlier. Moreover, most of the comparative studies between MERRA-2 and ground-based data are from plain topographic areas [45,46]. Among the 16 selected sites, there are 5 from AERONET, 1 from SKYNET (representing only one site from Hanle and Merak; see the details in Section 3) and 2 from SONET that have more than 4 years of AOD data. The details of the sites such as latitude, longitude and station altitude in metres (m) above mean sea level (amsl) are presented in

Table 1. Stations located in the HKH region as well as the three sub-regions of HKH considered in the present study for AOD trend estimation.

Site Name	Type of Sites	Lat.	Long.	Elevation (m)
HKH1	background & non-background	32–39	71–79	∼1000–5000
Hanle	background	32.78	78.95	4500
Merak	background	33.80	78.62	4310
Leh	background	34.15	77.58	3500
Issyk	non-background	42.60	76.98	1650
Srinagar	non-background	34.01	74.80	1590
Kashi	non-background	39.50	75.93	1320
HKH2	background & non-background	27–32	76–93	∼300–5000
EVK2	background	27.96	86.80	5080
Nam Co	background	30.77	90.96	4750
QOMS	background	28.37	86.95	4280
Lhasa	background	29.60	91.20	3680
Shimla	non-background	31.11	77.17	2280
Nainital	non-background	29.36	79.46	1940
Pokhara	non-background	28.19	83.97	800
Dehradun	non-background	30.32	78.03	430
HKH3	mostly background	28–36	93–103	∼1000–5000
Litang	background	29.98	100.26	3930
Chamdo	background	31.13	97.17	3240

. The classification of these sites in the HKH region is made roughly as background (or remote) and non-background (urban, semi-urban or significant anthropogenic/dust sources) from the observed yearly MERRA-2 AOD at 550 nm < 0.10 and >0.10, respectively. The non-background classification of the sites may have significant sources from anthropogenic (urban or semi-urban) or desert-dust origin. It is found that most of the sites selected in HKH1 and HKH3 sub-regions are remote/mountainous regions and far from the most polluted IGP regions unlike the HKH2 region.

The AERONET data used in the present work are Level 2.0, which refer cloud-screened and quality-controlled data products [34,44]. For SKYNET, we used AOD data for Hanle and Merak sites, located in the Ladakh region. In the data process, aerosol parameters (such as AOD, single scattering albedo, etc.) at five wavelengths are estimated using an inversion software package (SKYRAD.pack version 4.2) developed by [35]. Further, the estimated aerosol parameters are a performing cloud screening code [47] to remove the cloud-contaminated data. Therefore, AOD used from the SKYNET are cloud-screened and quality-controlled data products. The details of the instruments and their uncertainties on the estimated aerosol parameters are described by [48,49]. In the case of SONET data, to ensure the data quality, SONET established a routine calibration scheme, consisting of laboratory and field calibration approaches. The calibration of the instruments, multiwavelength-polarized sun–sky radiometer CE318-DP manufactured by Cimel Electronique in France, is performed every year in the field (Ling Mountain site, ∼1600 m amsl) via intercomparison with a master instrument as reported by [37]. These instruments are also routinely maintained by site managers along with the AERONET instruments in the region. The average AOD difference (0.002 ± 0.001) between SONET and AERONET is much smaller than the AERONET AOD uncertainty as reported by [37]. Further details of the data processing and quality assessment as well as specification of the instrument are reported therein. At present, the SONET has more than 20 sites, which are located in different typical regions of China, including urban, rural, desert, coastal, basin, mountain and plateau areas. Among these SONET sites, only two sites, Kashi and Lhasa, fall under the periphery of the HKH region. Therefore, we selected only these two sites from the network.

3. Data Analysis

The validation of AOD data from MERRA-2 are reported by several authors using ground-based measurement. For example, [45] have reported about the validation of ECMWF’s CAMS and MERRA-2 AOD at 550 nm using 15 years of AERONET observations from its 793 stations across the globe. In their analysis, MERRA-2 performs better than CAMS with root-mean-square error (RMSE) in the range 0.031–0.268 for CAMS and 0.017–0.232 for MERRA-2. [46] reported about the MERRA-2 AOD data with Sky radiometer observation at 550 nm from Hefei (31.90$°$ N, 117.17$°$ E, 30 m amsl), China, for the different seasons during March 2007 to February 2010. The correlation coefficients (R) were 0.88, 0.83, 0.88, and 0.80 in the spring, summer, autumn, and winter seasons, respectively. Further, Reference [50] has reported good agreement between MERRA-2 and AERONET AOD with correlation coefficients between 0.59 and 0.94 at nine sites in Australia. In the present work, a comparative study of MERRA-2 AOD at 550 nm was made with the available AERONET/SKYNET/SONET AOD data at 500 nm. In order to make a uniform comparison with MERRA-2 AOD at 550 nm, the AERONET/SKYNET/SONET AOD at 500 nm was converted into AOD at 550 nm using the relation:

$$AO{D}_{550}=AO{D}_{500}\ast {(500/550)}^{\alpha }$$

(1)

where $\alpha$ is the Angstrom Exponent, estimated from the slope of the spectral AOD plots in logarithmic scales. Figure 1 shows the scatter plot of monthly AOD from ground-based observation (AERONET/SKYNET/SONET) versus MERRA-2 data at 550 nm. Despite the discontinuity, the starting and end points of the data period with the available month data (N) are shown inside the Figure panel. Since AOD data from Hanle and Merak are very close with a difference in the 3rd decimal places [16], Merak AOD data during 2011 to 2015 can be added to missing data during this period at Hanle. Therefore, in the scatter plot, we used AOD data from 2008 to 2020 without any break. The same AOD at the two sites are also noticed in

Table 2. Yearly AOD (550 nm) trend (AOD year${}^{-1}$) at different sub-regions in the Hindu Kush Himalayan region in 1980–2020 using MERRA-2 reanalysis data. The results also show a Pearson correlation coefficient (R) and confidence level (t).

Stations	Mean AOD	Trend ( × 10${}^{-4}$)	R	t
HKH1	0.18 ± 0.05	23.33 ± 2.28	0.87	10.22
Hanle	0.06 ± 0.02	2.66 ± 0.93	0.45	2.85
Merak	0.06 ± 0.02	2.57 ± 0.96	0.43	2.67
Leh	0.08 ± 0.03	3.59 ± 1.04	0.52	3.47
Issyk	0.16 ± 0.04	2.85 ± 1.72	0.28	1.65
Srinagar	0.16 ± 0.04	15.85 ± 2.21	0.78	7.16
Kashi	0.29 ± 0.09	34.69 ± 3.88	0.85	8.95
HKH2	0.17 ± 0.04	32.20 ± 2.58	0.91	12.47
EVK2	0.06 ± 0.03	6.76 ± 1.00	0.77	6.74
Nam Co	0.06 ± 0.02	3.40 ± 0.92	0.55	3.71
QOMS	0.05 ± 0.02	4.28 ± 0.87	0.66	4.91
Lhasa	0.06 ± 0.03	4.76 ± 0.95	0.66	4.99
Shimla	0.29 ± 0.10	51.53 ± 4.99	0.88	10.33
Nainital	0.26 ± 0.09	53.15 ± 3.94	0.92	13.50
Pokhara	0.25 ± 0.08	55.75 ± 3.76	0.93	14.85
Dehradun	0.23 ± 0.08	39.16 ± 4.08	0.86	9.61
HKH3	0.09 ± 0.03	9.48 ± 1.21	0.81	7.82
Litang	0.06 ± 0.03	5.70 ± 0.90	0.75	6.34
Chamdo	0.06 ± 0.03	5.60 ± 1.01	0.70	5.53

. It may be mentioned that aerosol measurements at Hanle and Merak are performed by a same instrument (POM-01, Prede) by shifting the instrument at the two places during different time periods. Due to similar AOD at Merak and Hanle, we used the long-term data from 2008–2020 and represented as Hanle station in seasonal and climatological analysis. The duration of the AOD data from the available eight sites has minimum of 53 months (EVK2) to amaximum 140 (Hanle/Merak). In the scatter plot shown in Figure 1, the MERRA-2 AOD data from Hanle are found to be slightly higher than the ground measurement, i.e., 1.149 × SKYNET (slope) with offset 0.017 (intercept). Further, several sites also show lots of scatter data points between the two data products. For example, scatter data points at Kashi, Lhasa, Pokhara, etc. are seen in Figure 1. The MERRA-2 data at Pokhara have lesser than 0.206 × AERONET (or slope) with offset 0.133. Further, the MERRA-2 data at Lhasa and Kashi have 0.572 and 0.181 × SONET (or slope) with an offset of 0.017 and 0.234, respectively. The estimated correlation between the ground (AERONET/SKYNET/SONET) and MERRA-2 data is shown to be from minimum 0.31 (Kashi) to a maximum 0.76 (QOMS) from the monthly data as seen in Figure 1. Such large scatter data may attributed to the improper assumption of atmospheric parameters (such as surface albedo) over the non-uniform topography of mountain terrain used in the re-analysis data. Further, due to poor spatial resolution of the MERRA-2 data (0.5$°$ × 0.625$°$), it is likely unable to capture AOD within a small selected site as observed for other moisture parameters [51]. Similar to such discrepancies, the validation of satellite data over the high-altitude Tibetan Plateau shows low accuracy, as the observed data sometimes fall within the retrieval errors due to very lower AOD [41].

Due to a lack of consistency (or more scatter data) among the selected sites, we did not implement the ground-based offset to MERRA-2 data. Therefore, using a single data product from MERRA-2 provides uniform and consistent results in the long-term trend estimation. In addition, using such homogeneous long-term data will also provide a uniform pattern while estimating AOD trend over the Indian sub-continent. In order to test the statistical confidence level of the estimated trends, we estimated the Student’s t-test (t). In the process, t is calculated by:

$$t=R\times \frac{\sqrt{n-2}}{\sqrt{1-R^2}}.$$

(2)

Here, n is the number of the data set (e.g., 41 years) and R is the Pearson correlation coefficient in the linear trend analysis. From the estimated trend, if the absolute value of t is 1.965 or 1.648, then the trend is statistically significant at the 95% or 90% confidence level, respectively [52]. In the analysis, we removed high AOD data in 1982–1984 and 1991–1994. Such a high AOD was due to high impacts of two major volcanic eruptions, which happened in 1982 at El Chich$\acute{o}$n and 1991 at Mt. Pinatubo. The impacts of the two major volcanic period have lower impacts on lower sites (less than 2000 m amsl). However, at high-altitude sites (above 3000 m amls), the volcanic impacts were prolonged for a few more years than the lower sites, as seen in Supplementary Figures S1–S3.

4. Results and Discussion

4.1. AOD Trend in the Indian Sub-Continent

Figure 2 shows the AOD trend over the Indian sub-continent during the study period. The highest AOD trends in the sub-continent are observed in the IGP region (∼130.00 × 10${}^{-4}$ AOD year${}^{-1}$) and the estimated trends are statistically significant. There are high-AOD trends (above 100.00 × 10${}^{-4}$ AOD year${}^{-1}$) in HKH1 and HKH2 which are located in the closed peripheral area of the IGP region as seen in Figure 2. Due to longer study time period, the most of the estimated trends are statistically significant (t > 2.0) with high correlation coefficients of 0.95 (results of R and t are shown in Supplementary Figures S4 and S5. These results of AOD trends are consistent and agree with the work reported by several authors in the recent past. For example, [28] have reported increasing AOD trends from 88.00 × 10${}^{-4}$ to 266.00 × 10${}^{-4}$ AOD year${}^{-1}$ over the Indian sub-continent using ∼25 years of AOD data from several ARFINET (Aerosol Radiative Forcing India Network [6,28]) sites in the sub-continent. In addition, there are several papers that have reported about the increasing AOD trend over the Indian sub-continent using ground-based as well as satellite data. However, there are very few papers reported in the HKH1, HKH2, and HKH3 regions in particular due to lack of long-term continuous ground-based data. Further, the estimated trends vary depending upon the duration of the study period as well as continuity of the data as noticed in the present. Further, the impacts of the two major volcanic eruptions have been prolonged in the atmosphere for several years, and that has effect on high-altitude sites in particular in the trend estimation (see supplementary Figures S1–S3). Among the few papers available, the estimated AOD trend is statistically insignificant at the HKH1 region due to discontinuity or lack of long-term data [33]. The details of the AOD trends in yearly and seasonal trends among the HKH sub-regions are shown in the following sub-sections.

4.2. AOD Trend in the HKH Region

Due to a vast area with different topographic features, the estimated AOD trends among the sub-regions of HKH have different magnitudes. Among the sub-regions, HKH2 is located in the closed vicinity of the IGP region. However, HKH1 and HKH3 are located not far from the IGP region with a higher elevation as seen in Table 1. The estimated trends at few selected sites in the HKH region varied from the minimum 2.57 × 10${}^{-4}$ ± 0.96 × 10${}^{-4}$ AOD year${}^{-1}$ at Merak in the HKH1 to maximum 55.75 × 10${}^{-4}$ ± 3.76 × 10${}^{-4}$ AOD year${}^{-1}$ at Pokhara at HKH2 sub-regions of HKH during 1980–2020. The estimated trends located at high-altitude (above 3000 m amsl) or far from IGP such as Hanle, Merak, and Issyk exhibit weak AOD trends. However, the estimated trends near the IGP region are strong and statistically significant (t > 4.0). Figure 3, Figure 4 and Figure 5 show the increasing AOD trends at a few selected sites that have higher and lower trends in 1980–2020. The estimated trends, along with t and correlation coefficient (R), are shown in the graphs of each site. The HKH, HKH1, HKH2 and HKH3 sub-regions exhibit 23.33 × 10${}^{-4}$ ± 2.28 × 10${}^{-4}$, 32.20 × 10${}^{-4}$ ± 2.58 × 10${}^{-4}$, and 9.48 × 10${}^{-4}$ ± 1.21 × 10${}^{-4}$ AOD year${}^{-1}$, respectively. The HKH3 region has a very low number of ground-observing stations despite an increasing AOD trend, as seen in Figure 2. It is seen that the AOD trend starts decreasing after 2018 onward at most of the sites (except Issyk and Kashi) shown in Figure 3 and Figure 5 and Supplementary Figures S1–S3. Further, there were significantly lower AOD values in 2020, which is apparently due to continuous national lockdown from the COVID-19 pandemic situation. Due to the national lockdown, air quality across the Indian sub-continent improved significantly in 2020, as reported in several papers [53,54,55]. This reduction in AOD has been significantly influential in the close vicinity of the IGP region of the HKH. However, the lockdown has a lower impact on sites that are far from the IGP region and that may be one of the reasons for increasing AOD at Issyk and Kashi even after the lockdown period. However, the rising AOD at these two sites may be due to local effects from the surrounding desert, Taklamakan. It may be mentioned that the Taklamakan desert is one of the largest sandy deserts in the world. Apart from the MERRA-2 AOD data, an increasing AOD in 2020 at Kashi is also observed from the SONET ground-based data (not shown here). The details of the AOD trends in each sub-section of HKH are described in the following sub-sections.

4.2.1. HKH1

There are six potential sites in the HKH1 region, namely Hanle, Merak, Issyk, Leh, Srinagar and Kashi, although Issyk and Kashi fall near the periphery of the HKH1 region. Among these sites, Kashi site, which is very close to the Taklimakan desert, is often affected by desert dust, as reported by [37]. The estimated trends at HKH1 varied from a minimum 2.57 × 10${}^{-4}$ ± 0.96 × 10${}^{-4}$ AOD year${}^{-1}$ at Merak to a maximum 34.69 × 10${}^{-4}$ ± 3.88 × 10${}^{-4}$ AOD year${}^{-1}$ at Kashi, as seen in Table 2. The estimated AOD trends at Leh, Hanle and Srinagar are 3.59 × 10${}^{-4}$ ± 1.04 × 10${}^{-4}$, 2.66 × 10${}^{-4}$ ± 0.93 × 10${}^{-4}$ and 15.85 × 10${}^{-4}$ ± 2.21 × 10${}^{-4}$ AOD year${}^{-1}$, respectively, with statistical significance (t > 2). However, there are high AOD trends in the southwestern part of HKH1, as seen in Figure 2. The high AOD trends in this sub-region are mainly due to major contribution from this southwestern region. Among the selected sites, Srinagar is very close to this sub-region. Apparently, the present results of AOD trends at Srinagar and Leh towns using such a long-term data are the first reported results so far from the region. Moreover, the two towns have slightly more AOD trends than Hanle, as expected, due to the significant contribution of anthropogenic pollution from the major tourist activities as well as transported from the IGP region. Among the sites, Issyk has the lowest AOD trend 2.85 × 10${}^{-4}$ ± 1.72 × 10${}^{-4}$ AOD year${}^{-1}$, despite the site is falling beyond the HKH1 region. Due to the weak trend, the estimated trend is statistically insignificant, as seen in Table 2. Although both Issyk and Kashi are located above 1000 m amsl, Kashi is located in a closer vicinity to Taklamakan desert than Issyk, and consequently, Kashi is often affected by desert dust, where AOD is obviously higher (0.56, from SONET data) than the surrounding desert [37]. Since yearly AOD (MERRA-2) at Srinagar, Issyk and Kashi are more than 0.10 (see Table 1), we consider these sites as non-background.

4.2.2. HKH2

In the HKH2 region, we selected eight potential sites, and several sites are under AERONET. However, there are only sites that have more than 4–5 years of continuous aerosol measurement. Further, many of these sites are located in the close vicinity of the IGP region. Among these sites, Nam Co and Lhasa are located far from the IGP region, as seen in Figure 2. Nam Co is an AERONET station, while Lhasa falls under SONET. Due to low AOD (>0.1 from MERRA-2), EVK2, Nam Co, QOMS and Lhasa in this sub-region are considered as background sites, as seen in Table 1. Here, it may be mentioned that IGP region is located in one of the most densely populated regions in the northern India and is also one of the most polluted regions in the Indian sub-continent [22,24,56,57]. The region has the highest estimated trends (55.75 × 10${}^{-4}$ ± 3.76 × 10${}^{-4}$ AOD year${}^{-1}$, at Pokhara) among the sub-regions of HKH. This indicates that there is increasing significant contribution of aerosols, apparently due to the increase in anthropogenic sources from local as well as long-range transportation from the IGP region. Among the selected eight potential sites, the lowest AOD trend 3.40 × 10${}^{-4}$ ± 0.92 × 10${}^{-4}$ AOD year${}^{-1}$ is observed at Nam Co, followed by 4.28 × 10${}^{-4}$ ± 0.87 × 10${}^{-4}$, 4.76 × 10${}^{-4}$ ± 0.95 × 10${}^{-4}$, and 6.76 × 10${}^{-4}$ ± 1.00 × 10${}^{-4}$ at QOMS, Lhasa and EVK2, respectively. These estimated trends are statistically significant, as seen in Table 2. In an earlier work, [33] have reported 10.00 × 10${}^{-4}$ AOD year${}^{-1}$ at EVK2, and 11.00 × 10${}^{-4}$ AOD year${}^{-1}$ at Nam Co and QOMS using AERONET data. Due to lack of long-term data, the estimated trends were weak and statistically insignificant. However, the present study has improved the estimated trends. The remaining sites, which show significant enhancement of AOD in the present study, are Dehradun, Shimla and Nainital with 39.92 × 10${}^{-4}$ ± 4.08 × 10${}^{-4}$, 51.53 × 10${}^{-4}$ ± 4.99 × 10${}^{-4}$ and 53.15 × 10${}^{-4}$ ± 3.94 × 10${}^{-4}$, respectively. These estimated trends are significantly strong and statistically significant with R of 0.86 to 0.93 due to the presence of the homogeneous long-term data set, as seen in Table 2.

4.2.3. HKH3

We selected two potential sites in the HKH3 region, namely Litang and Chamdo in the Tibet Autonomous region in China. Chamdo is Tibet’s third largest city after Lhasa and Shigatse. There are also limited sites from AERONET or other networks, unlike the HKH2 region, while Litang, located in the eastern edge of the Tibetan Plateau (TP), is an AERONET site under the Chinese Academy of Sciences. However, the site has few days of observation and has been kept only for temporary testing (after communication from the AERONET site owner). Due to low AOD (>0.1 from MERRA-2), both of these sites are considered as background sites (see Table 1). The estimated trends are 5.60 × 10${}^{-4}$ ± 1.01 × 10${}^{-4}$ and 5.70 × 10${}^{-4}$ ± 0.90 × 10${}^{-4}$ AOD year${}^{-1}$ at Chamdo and Litang, respectively, while HKH3 region exhibits 9.48 × 10${}^{-4}$ ± 1.21 × 10${}^{-4}$ AOD year${}^{-1}$. Due to continuous long-term and homogeneous data, the estimated trends are strong and statistically significant. With relevance to the present increasing AOD trends, [58] have also reported increasing AOD trends (440 nm) in the TP by 10.00 × 10${}^{-4}$ ± 30.00 × 10${}^{-4}$ at Lhasa, 130.00 × 10${}^{-4}$ ± 30.00 × 10${}^{-4}$ at Mt_WLG and 20.00 × 10${}^{-4}$ ± 20.00 × 10${}^{-4}$ at Nam CO using ground-based data during the past decade. The positive trends in most TP areas may be due to the rapid increase in human activities, such as the expansion of tourism to the TP and biomass burning in South Asia, as reported therein.

4.2.4. Seasonal Long-Term AOD Trend

Table 3. Seasonal (Winter and Autumn) AOD trend ( AOD year${}^{-1}$), R and t at different sub-regions in the Hindu Kush Himalayan region using MERRA-2 data om 1980–2020 in the present study.

	Winter			Autumn
Stations	Trend ( × 10${}^{-4}$)	R	t	Trend ( × 10${}^{-4}$)	R	t
HKH1	21.99 ± 2.04	0.89	10.78	29.79 ± 2.96	0.87	10.97
Hanle	1.40 ± 0.69	0.34	2.02	0.90 ± 0.71	0.22	1.29
Merak	0.24 ± 0.49	0.09	0.48	0.44 ± 0.72	0.11	0.61
Leh	2.17 ± 0.78	0.44	2.79	2.36 ± 0.76	0.48	3.09
Issyk	1.88 ± 2.29	0.15	0.82	5.74 ± 1.47	0.57	3.91
Srinagar	23.68 ± 1.92	0.91	12.32	20.57 ± 3.16	0.76	6.52
Kashi	29.03 ± 4.88	0.72	5.94	37.75 ± 5.79	0.76	6.52
HKH2	41.40 ± 2.47	0.95	16.74	38.47 ± 3.27	0.90	11.78
EVK2	0.36 ± 0.85	0.08	0.43	3.47 ± 0.92	0.55	3.77
Nam Co	–0.36 ± 0.58	–0.11	–0.61	2.41 ± 0.97	0.40	2.48
QOMS	–0.56 ± 0.68	–0.15	–0.83	1.90 ± 0.54	0.53	3.51
Shimla	55.91 ± 4.02	0.93	13.90	51.30 ± 5.35	0.86	9.59
Nainital	50.80 ± 3.50	0.93	14.48	54.27 ± 7.97	0.77	6.81
Pokhara	63.56 ± 5.53	0.90	11.51	48.97 ± 4.33	0.90	11.32
Dehradun	40.65 ± 3.54	0.89	11.47	36.61 ± 4.16	0.84	8.79
HKH3	7.79 ± 1.22	0.74	6.38	3.66 ± 0.67	0.69	5.46
Litang	3.80 ± 0.71	0.69	5.35	0.97 ± 0.68	0.24	1.43
Lhasa	0.18 ± 0.70	0.05	0.26	2.90 ± 1.03	0.45	2.82
Chamdo	4.28 ± 1.06	0.58	4.03	2.98 ± 0.75	0.58	3.99

and

Table 4. Seasonal (Spring and Summer) AOD trend (AOD year${}^{-1}$), t and R values at different sub-regions in the Hindu Kush Himalayan region during 1980–2020 in the present study.

	Spring			Summer
Stations	Trend ( × 10${}^{-4}$)	R	t	Trend ( × 10${}^{-4}$)	R	t
HKH1	10.57 ± 1.71	0.73	6.19	34.68 ± 4.10	0.82	8.45
Hanle	–0.35 ± 1.25	–0.05	–0.28	0.30 ± 1.77	0.03	0.17
Merak	1.76 ± 1.62	0.19	1.08	5.82 ± 1.99	0.47	2.93
Leh	1.80 ± 1.89	0.17	0.95	7.13 ± 0.85	0.83	8.37
Issyk	9.63 ± 4.25	0.37	2.27	–1.90 ± 2.97	–0.11	–0.64
Srinagar	1.64 ± 2.59	0.11	0.63	22.18 ± 2.54	0.84	8.72
Kashi	32.19 ± 4.10	0.81	7.85	23.56 ± 6.04	0.57	3.90
HKH2	26.10 ± 3.62	0.78	7.22	33.82 ± 2.05	0.94	16.52
EVK2	14.09 ± 2.06	0.77	6.83	9.05 ± 1.14	0.81	7.95
Nam Co	6.37 ± 1.62	0.57	3.93	5.92 ± 0.81	0.80	7.31
QOMS	9.14 ± 1.82	0.67	5.04	7.65 ± 0.90	0.83	8.46
Shimla	27.89 ± 3.91	0.78	7.14	70.56 ± 6.33	0.89	11.13
Nainital	34.95 ± 4.39	0.82	7.97	72.95 ± 4.88	0.94	14.93
Pokhara	65.80 ± 7.45	0.84	8.84	41.45 ± 5.14	0.82	8.06
Dehradun	19.82 ± 3.38	0.72	5.86	53.31 ± 5.87	0.85	9.09
HKH3	19.62 ± 1.66	0.90	11.83	7.75 ± 1.51	0.68	5.13
Litang	11.67 ± 1.33	0.83	8.73	8.17 ± 1.04	0.81	7.86
Lhasa	6.54 ± 1.34	0.65	4.89	5.98 ± 0.90	0.76	6.63
Chamdo	11.26 ± 1.16	0.86	9.71	2.75 ± 1.12	0.40	2.45

show seasonal AOD trends (AOD year${}^{-1}$) at 550 nm, t and R over few selected sites located in the sub-regions of HKH during 1980–2020. It was found that all the high-AOD groups located in the close vicinity of IGP region exhibited a significant increasing trend with statistical significance in all the seasonal data. Among the high-AOD groups, Nainital shows the strongest positive trends with statistical significance, that is, a maximum 72.95 × 10${}^{-4}$ ± 4.88 × 10${}^{-4}$ AOD year${}^{-1}$ in summer (June–July–August) and minimum 34.95 × 10${}^{-4}$ ± 4.39 × 10${}^{-4}$ AOD year${}^{-1}$ in spring (March–April–May), as seen Table 3 and Table 4. The minimum positive trend in spring was relatively lower than winter and autumn, and the reason for this phenomenon needs to be explored. The increasing AOD trends over the region may be attributed to rapid economic development and large-scale urbanization since the late twentieth century, in addition to enhancement from long-range transportation process from the Middle East, the Arabian Sea and the Thar Desert [16,17,59,60].

The present study found that most of the high-altitude sites showed weakly increasing AOD trends, except for a few negative trends, during all the seasons. The results are mixed, with some statistically significant and some and non-statistically significant. For example, there are positive trends of 1.80 × 10${}^{-4}$ ± 1.89 × 10${}^{-4}$ in spring, 7.13 × 10${}^{-4}$ ± 0.85 × 10${}^{-4}$ in summer, 2.36 × 10${}^{-4}$ ± 0.76 × 10${}^{-4}$ in autumn and 2.17 × 10${}^{-4}$ ± 0.78 × 10${}^{-4}$ in winter at Leh. These trends, except for the spring one, are statistically significant. Further, there is a negative AOD trend at Hanle, including –0.35 × 10${}^{-4}$ ± 1.25 × 10${}^{-4}$ in spring, and increasing AOD trends, including 0.30 × 10${}^{-4}$ ± 1.77 × 10${}^{-4}$ in summer, 0.90 × 10${}^{-4}$ ± 0.70 × 10${}^{-4}$ in autumn and 1.40 × 10${}^{-4}$ ± 0.69 × 10${}^{-4}$ in winter. These trends, except for the winter one, are statistically non-significant. Figure 6 shows seasonal AOD trends at each representative site of the sub-HKH three regions. Such weak trends in the high-altitude sites could show slight improvement in air quality during the national lockdown period coupled with the decreasing AOD trend in the previous few years. The present study found slightly decreasing AOD trends during the last 5–6-year period at these high-altitude sites in particular, as seen in Figure 3 to Figure 5. However, such variation is not significant at the sites located in the close vicinity of the IGP region. Further, the impacts of two major volcanic eruptions happening during 1982 and 1991 at El Chich$\acute{o}$n and Mt. Pinatubo, respectively, weakened in the trend estimation, particularly at high-altitude sites. Similarly to in this present work, there are weak and statistically insignificant trends with positive as well as negative trends over the high-altitude mountainous regions in the previous studies [33]. Further, [28] have also reported negative trends of AOD near the foothills of the Himalayas at Gandhi College, Ballia, in India using ground-based data in the recent past. The present study results of all positive AOD trends in all the foothills of the Himalayan region indicate a strong enhancement of AOD during the last decade due to rapid economic development, energy demands and large-scale urbanization/industrialization.

4.3. Climatological AOD

Figure 7 shows a violin plot of climatological AOD at 16 selected sites, as well as the three sub-regions of HKH using MERRA-2 data during the years 1980–2020. The central solid circle in white indicates the average AOD, while top and bottom of solid circles in black represent maximum and minimum values, respectively. Due to large variation in AOD magnitude, the Y-axis shown in Figure 7 cannot make a uniform upper limit. These figures clearly show a distinct seasonal variation with peak AOD values (average) during spring (March–May), as well as summer (June–July), at different sub-regions of HKH. These include a few selected high-AOD groups, including 0.48, 0.44 and 0.42 at Shimla, Nainital and Pokhara, respectively, in June–July, and at a few selected low AOD-groups, including 0.08 and 0.06 at Leh and Hanle, respectively, in April–May. The minimum AODs are observed during winter (December–January–February) at 0.05–0.07 and 0.02 for the few selected high and low AOD groups, respectively. However, on a wider scale, the magnitude of the maximum AOD values were 0.25, 0.25 and 0.15 in HKH1, HKH2 and HKH3, respectively. The minimum AOD observed at the three sub-regions of HKH varied from 0.03 to 0.05. The yearly average AOD values at the present study sites are listed in Table 2. The observed AOD values reported at Nainital are consistent with ground observational data reported in the past [59,61]. In the case of Hanle/Merak, the seasonal AOD peaked as ∼ 0.05–0.06 in March–May (spring) according to ground observation [16,18,62], and such results from ground observation are consistent with the present reanalysis data. The high AOD during spring/summer at the background sites Hanle/Merak is attributed to long-range transported aerosols [16,18]. The prevailing wind at Hanle/Merak was shown to be southwesterly. However, during the Indian monsoon months (June–August), there are two types of air mass arrivals at the study region, one from the northeast and another from the north-west, which bring anthropogenic aerosols [62].

5. Conclusions

The present study found significant increasing AOD trends in the IGP region in the Indian sub-continent (∼130 × 10–04 AOD year${}^{-1}$), and the trends are statistically significant, with correlation coefficient above 90%. However, the increasing AOD trends over the background site and climate-sensitive region of the HKH are one of the main concerns in the present study. Due to a lack of such ground-based data, expanding the ground-based network is needed in order to challenge severe threats from global and regional warming in the context of climate change. For the first time, we present the long-term (1980–2020) AOD trends over the HKH region using homogeneous and continuous long-term MERRA-2 reanalysis data. The main findings of the AOD trend over the HKH region are as follows:

• The estimated AOD${}_{550nm}$ trends over the HKH1 varied from a minimum 2.57 × 10${}^{-4}$ ± 0.96 × 10${}^{-4}$ AOD year${}^{-1}$ at Merak (Ladakh) to a maximum 34.69 × 10${}^{-4}$ ± 3.88 × 10${}^{-4}$ AOD year${}^{-1}$ at Kashi. Among the HKH2 sites, significant enhancement of AOD is observed at many stations, namely, Dehradun, Shimla, Nainital and Pokhera by 39.16 × 10${}^{-4}$ ± 4.08 × 10${}^{-4}$, 51.53 × 10${}^{-4}$ ± 4.99 × 10${}^{-4}$, 53.15 × 10${}^{-4}$ ± 3.94 × 10${}^{-4}$ and 55.75 × 10${}^{-4}$ ± 3.76 × 10${}^{-4}$ AOD year${}^{-1}$, respectively. However, the estimated AOD trends at HKH3 varied from minimum 4.76 × 10${}^{-4}$ ± 0.95 × 10${}^{-4}$ at Lhasa to maximum 5.70 × 10${}^{-4}$ ± 0.90 × 10${}^{-4}$ AOD year${}^{-1}$ at Litang with statistical significance.
• At the sub-regional scale, HKH1, HKH2 and HKH3 exhibit 23.33 × 10${}^{-4}$ ± 2.28 × 10${}^{-4}$, 32.20 × 10${}^{-4}$ ± 2.58 × 10${}^{-4}$ and 9.48 × 10${}^{-4}$ ± 1.21 × 10${}^{-4}$. AOD year${}^{-1}$, respectively, with correlation (R) of 0.81 to 0.91. The estimated trends are statistically significant, and the high AOD trends in the HKH1 and HKH2 region are due to long-range transported anthropogenic sources from the IGP region.
• On a seasonal basis, Pokhara, Nainital, Shimla and Dehradun show strong positive AOD trends with minimum 19.82 × 10${}^{-4}$ ± 3.38 × 10${}^{-4}$ to maximum 72.95 × 10${}^{-4}$ ± 4.88 × 10${}^{-4}$ AOD year${}^{-1}$, during all the seasons, with statistical significance. In addition, there are also increasing AOD trends at all the high-altitude background sites. However, the estimated trends are weak and statistically insignificant in several cases. Such weak trends in the high-altitude sites could be due to slight improvement of air quality during the national lockdown period coupled with the decreasing AOD trend in the previous few years. However, such variation is not significant at the sites located in the close vicinity of the IGP region. Further, the impacts of two major volcanic eruptions happening during 1982 and 1991 at El Chich$\acute{o}$n and Mt. Pinatubo, respectively, shows weakening of the trend estimation, particularly at high-altitude sites.
• The weakening AOD trends at the high-altitude sites in particular, which are attributed to a slight reduction in AOD during the last 5–6 years, may need to be re-examined with a longer and homogeneous data sets in future.