Transboundary sources dominated PM2.5 in Thimphu, Bhutan.

This study estimates the potential source regions contributing to PM2.5 in the capital city of Thimphu, Bhutan, during the years 2018-2020 using the ground-based data, followed by the HYSPLIT back trajectory analysis. The average PM2.5 concentration in the entire study period was 32.47 µg/m3 which is three times of the World Health Organization recommended limit of 10 µg/m3. Less than half of the days in pre-monsoon (43.47%) and post-monsoon (46.41%), and no days in winter were within the 24-h average WHO guideline of 25 μg/m3. During the COVID-19 lockdown imposed from August 11 to September 21 in Bhutan, only a marginal reduction of 4% in the PM2.5 concentrations was observed, indicating that nonlocal emissions dominate the PM2.5 concentrations in Thimphu, Bhutan. Most back trajectories in the analysis period were allocated to south or south-west sector. India was the major contributor (~ 44%), followed by Bangladesh (~ 19%), Bhutan itself (~ 19%) and China (~ 16%). This study confirms that there are significant contributions from transboundary sources to PM2.5 concentrations in Thimphu, Bhutan, and the elevated PM2.5 concentrations need to be tackled with appropriate action plans and interventions. © Islamic Azad University (IAU) 2021.

Introduction

Air pollution has become a serious issue that can pose a threat to human health as well as the environment. In 2016, ambient air pollution caused 7.6% of all deaths (WHO 2016). Due to insufficient control of emissions from rapidly increasing industries, urbanization, and energy consumption, health effects associated with air pollution in developing countries are severe (Guo et al. 2018; Cohen et al. 2005). The situation is more severe in lower–middle-income countries (LMICs), which house around 40% of the global population, when compared with high-income countries (HICs) (Sharma et al. 2020a). It has been estimated that outdoor air pollution leads to 3.3 million premature deaths per year worldwide, with PM2.5 alone causing around 3.15 million deaths per year predominantly in Asia (Lelieveld et al. 2015; Babatola 2018). Fine particulate matter (PM2.5) which are particles with aerodynamic diameter less than 2.5 μm pose great danger to health and are a resultant of different types of combustion activities (motor vehicles, power plants, wood burning, etc.) and certain industrial processes (Sahu et al. 2020; Chen et al. 2016; Pope and Dockery 2006; Morawska and Zhang 2002).

Air quality in Bhutan is believed to be pristine owing to its vast forest cover, but there has been a vague understanding of how the air quality was affected due to the massive development in the country. According to National Environment Commission (NEC) (NEC 2016), Bhutan, air pollution is increasing over time, creating a significant concern to the health of environment and people living in it. NEC, Bhutan, also reported that the air quality of the capital city of Bhutan, Thimphu, is much worst during winter seasons owing to the increasing use of firewood for cooking, forest fires, and the burning of agriculture debris (NEC 2016). Poor air quality observed in 2016 led to the implementation of various national measures like the revision of vehicle emission standards, promotion of electric vehicles, and people were encouraged to use better wood stoves to reduce emission due to firewood for cooking (Zheng et al. 2019). The climate change policy drafted by the NEC is also expected to curb pollution in the future.

Bhutan's concern with pollution was mostly indoor until the last decade because most of the households did not have access to modern forms of energy, which increased their exposure to indoor air pollution and related health problems, such as bronchitis, asthma, and miscarriage (BMCI 2016). Along the southern boundary that the country shares with India and in Bhutan's urban and industrial areas, air quality has been found to continuously degrade (NEC 2016). However, very few studies, if any, were carried out on the air quality of Bhutan with a focus to understand if the pollutants are of local or non-local origin. Also, contribution of air pollution to mortality and morbidity in Bhutan is not well established. Thus, the main aim of this study is to understand the status of PM2.5 pollution and possible source origins in the capital city of Bhutan.

Materials and methods

Study area and data sources

Thimphu, the capital and the largest city of Bhutan, was selected for this study as it is the only city in Bhutan for which historical ground-based PM2.5 data were available. It is situated in the central part of Bhutan at 27.4712°N, 89.6339°E at about 7000 feet (2000 m) above sea level and is the most populated city in Bhutan with over more than 1 lakh people.

The study period of April 2018 to June 2020 was divided into four seasons, viz. pre-monsoon (March, April, and May), monsoon (June, July, and August), post-monsoon (September, October and November), and winter (December, January and February). For this study, archived global reanalysis data from National Oceanic and Atmospheric Administration (NOAA) (ftp://arlftp.arlhq.noaa.gov/archives/reanalysis/) were used as the meteorological input data for back trajectory analysis. The PM2.5 data based on regionally interpolated real-time monitoring station data were downloaded from Berkeley Earth (http://berkeleyearth.lbl.gov/air-quality/local/Bhutan/Thimphu/Thimphu).

Back trajectory

The Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model (Stein et al. 2016) developed by the NOAA was used to calculate air mass trajectories. The model calculation method comprises of a combination of a moving frame of reference that follows air parcels as they move from their initial location for the simulation of advection, diffusion, and deposition (Lagrangian) and a fixed three-dimensional grid as a frame of reference to compute the pollutant air concentrations (Eulerian) (Rolph et al. 2017; Stein et al. 2016). The model uses a puff and particle approach to estimate the horizontal and vertical dispersion, respectively, as explained by Eqs. (1) and (2) where X, X1, and U denote initial, first guest positions, and velocity vector, respectively (Sahu et al. 2019).

Seventy-two hours back trajectory was calculated with the consideration that the residence time of a PM2.5is around three days (WHO 2006b; Jaenicke 1982; Muller 1982). The HYSPLIT model was run in R using the Openair package (Carslaw and Ropkins 2012).

Starting height plays a key role in the computation of back trajectories. Trajectories with very low starting height may hit the ground early and thus would result in information loss. On contrary, very high starting height could result in trajectories not entering the mixing layer. Optimum height for starting a trajectory would be the middle of the planetary boundary layer, which keeps varying (Huang et al. 2013; Sargent et al. 2018). For this study, the starting height was assumed to be equal to the half of the boundary layer height. A similar approach was used by Huang et al. (2013) in which HYSPLIT was run to understand the link between precipitation and high Hg wet deposition events for which back-trajectories were calculated with the starting height as half the planetary boundary layer height.

Back trajectory analysis

Trajectory gridded frequencies analysis was carried out in which the number of unique trajectories in a particular grid square was considered and hexagonal binning was used to gain an idea about trajectory frequencies in a particular grid. The back trajectories obtained for the study period were allocated to eight different sectors, i.e. north (N), north-east (NE), east (E), south-east (SE), south (S), south-west (SW), west (W) and north-west (NW), based on the time a particular trajectory is present in a specific sector. A trajectory is then allocated to a sector if it spends at least 50% (Carslaw and Ropkins 2012) of its time in that sector; otherwise, it is categorized as ‘unallocated’.

Concentration weighted trajectory (CWT)

Seibert et al. (1994) computed concentration fields to identify source areas of pollutants. This approach is sometimes referred to as the CWT or CF (concentration field). It is a method of weighting trajectory residence time with associated pollutant concentrations. This procedure assigns to each grid cell a weighted concentration obtained by averaging sample concentrations associated with trajectories that crossed that grid cell, i.e. each concentration is used as a weighting factor for the residence times of all the trajectories in each grid cell which correspond to that concentration, and then, it is divided by the cumulative residence time from all trajectories. CWT method employs an arbitrary weight function to minimize the inaccuracy caused by the small number of polluted trajectories. In summary, weighted concentration fields show concentration gradients across potential sources and helps to identify the relative significance of potential sources (Reizer and Orza 2018; Stohl 1996).

For each grid cell, the mean (CWT) or logarithmic mean (used in the residence time weighted concentration (RTWC) method) concentration of a pollutant species was calculated using the following equation: where ‘i’ and ‘j’ are the indices of grid, ‘k’ is the index of trajectory, ‘N’ denotes the total number of trajectories used in analysis, c is the pollutant concentration measured upon arrival of trajectory ‘k’, and τ is the residence time of trajectory k in grid cell (i, j). A high value of C means that air parcels passing over cell (i, j) would cause high concentrations at the receptor site.

Identifying the contribution of high concentration back trajectories

An estimate of the contributions of PM2.5 from the regions around Bhutan can help us understand the contributions and dynamics of the regional/local transport. It can also be helpful in determining appropriate strategies in tackling high pollutant concentrations in the capital city of Bhutan.

The PM2.5 contribution to a grid from neighbouring grid is estimated using Eq. (4) as given by (Cheng et al. 2013; Sahu et al. 2019) where Ccs indicates PM2.5 contribution to grid ‘c’ from grid ‘s’; Cs indicates average PM2.5 concentration associated with trajectory endpoints passing through or originating from region s, and T indicates the total number of back trajectory endpoints. (Both Ts and T do not include trajectory endpoint at hour = 0.)

Results and discussion

Air pollution is severe in Thimphu, Bhutan

Change in the monthly PM2.5 concentrations in the study period is shown in Fig. 1. The average concentration in the study period was 32.5 µg/m3, which is three times of the WHO recommended limit of 10 µg/m3 (WHO 2006a). A clear increase in concentrations in the months of September to March, with peaks in January (2020) (182.1 µg/m3), December (2019) (170.4 µg/m3), November (2018) (168.4 µg/m3), and February (2020) (100.6 µg/m3), was observed. The averaged concentrations in pre-monsoon, monsoon, post-monsoon, and winter during the entire study period were 30.4, 15.8, 29.7, and 55.6 µg/m3, respectively. 43.5, 98.2, and 46.4% days in pre-monsoon, monsoon, and post-monsoon were within the 24-h average WHO limit of 25 μg/m3. However, in winters, all days exceeded the 24-h average WHO recommended PM2.5 concentration limit of 25 μg/m3.

Fig. 1

Monthly variation of average PM2.5 (μg/m3) concentration. The red dotted line represents the WHO annual average criteria for PM2.5. The central line in the box represents the interquartile range (IQR) (first (25 percentile, Q1) to third quartile (75 percentile (Q3)), the dark horizontal line in the box shows the median concentrations, the extended horizontal lines at the end of dashed vertical lines represent the maximum and minimum values, and the circles represent the outliers (criteria: minimum = (Q1 − 1.5 × IQR) and maximum = (Q3 + 1.5 × IQR))

Daily averaged PM2.5 concentrations exceeded the WHO guideline of 25 μg/m3 on 38.3, 63.4, 39.4% days in 2018, 2019, and 2020, respectively. Among all the months, the 24-h averaged PM2.5 concentration in August of 2018 and 2019 did not exceed the guideline on any day, unlike January of 2019 and 2020, where the concentrations exceeded the limit on all days.

December had the maximum PM2.5 concentrations in 2018 throughout the day with the highest concentrations observed at 4 am (65.2 µg/m3) and 5 pm (72.4 µg/m3). In 2019, the month of January had the maximum PM2.5 concentration during early mornings at 4 am (76.0 µg/m3) and December had the maximum PM2.5 concentration with peak during the evening at 5 p.m. (90.7 µg/m3). In 2020, the highest concentrations were observed for the month of January with the peaks observed at 5 a.m. (63.7 µg/m3) in the morning and 5 p.m. (92.2 µg/m3) in the evening. Therefore, it can be concluded that in winter season, particularly December and January months, concentrations peaked during mornings and evenings, with the ratio of lowest to highest PM2.5 concentrations being 2.2 for December of 2018, 2.3 for December of 2019, 2.4 for January of 2019, and 2.4 for January of 2020. The 24-h mean values were lowest for the month of August in both 2018 (15.7 µg/m3) and 2019 (12.0 µg/m3) for which the hourly values of 24-h average PM2.5 concentrations were found to be under the WHO recommended PM2.5 concentration limit of 25 μg/m3.

Local emission changes have minor effects

Due to COVID-19 pandemic, lockdowns were enforced in many countries across the globe including countries like India, China, France, Italy, New Zealand, Poland, and the UK, where world’s largest and most restrictive mass quarantines were implemented (Sharma et al. 2020b). Bhutan confirmed its first COVID-19 case on 6 March 2020. The first nationwide lockdown in Bhutan was enforced on 11 August 2020 which ended on 21 September 2020. During the lockdown, travel and goods transport between districts were banned. All schools, nonessential government facilities, and commercial establishments were closed.

To assess the effect of these restricted emissions, PM2.5 concentrations during the lockdown phase of August 11 to September 21, for the years 2018, 2019, and 2020, were analysed (Fig. 2). It was observed that the average concentration for the lockdown phase for 2018 and 2019 combined was 15.1 µg/m3, and in 2020 was 13.6 µg/m3, indicating a marginal ~ 4% reduction in the PM2.5 concentrations due to lockdown in 2020. However, from August 11 to 30, a significant reduction of 18.7% in the PM2.5 concentrations was observed. This analysis shows that local emission sources do not affect the total PM2.5 concentrations much. This suggests that regional transport of PM2.5 is mainly responsible for the high PM2.5 concentrations in Thimphu, Bhutan.

Fig. 2

Percentage reduction in the PM2.5 concentrations during the lockdown period in 2020 when compared to the average concentrations observed in 2018 and 2019. The dashed red line shows the average percentage reduction (3.71%) observed during the lockdown period. The positive values show reduction in the PM2.5 concentrations in the year 2020, and the negative values indicate an increase in the PM2.5 concentrations

Transboundary sources dominated PM2.5

Trajectories allocated to wind sectors

The trajectories obtained by running HYSPLIT for the study period were allocated to eight wind sectors as shown in Fig. 3. It was observed that in the year 2018, on an average 49 μg/m3 of PM2.5 was contributed majorly from the trajectories allocated to west sector (17% of the trajectories) and 42.4 μg/m3 was allocated to trajectories allocated to the northwest (15% of the trajectories) sector. 56.3% of the trajectories for the year 2018 were not allocated to any sector. In the year 2019, 50.5 μg/m3 of PM2.5 was contributed from trajectories allocated to northwest sector, followed by the southwest sector. For the year 2020, trajectories allocated to southeast and southwest sectors were observed to be the dominant sectors and contributed 63.2 μg/m3 and 53.4 μg/m3 of PM2.5, respectively. However, from the results it can be concluded that for the year 2018, south was the dominant wind directions with majority of the trajectories originating from southern sector. For the year 2019, south and west, for 2020, south-west were observed to be the dominant sector with most of the trajectories originating from these sectors. Thus, from the allocation analysis, it is evident that majority of the trajectories reaching Thimphu, Bhutan, originate from the southern neighbouring regions of Bhutan.

Fig. 3

Figure showing the proportion of time in a year, a back trajectory resides in a specific sector and the mean concentration associated with the back trajectories associated with the specific sector

Back Trajectory Analysis

The number of trajectories originating from different regions in various seasons is plotted in Fig. 4. It was observed that for the pre-monsoon months, the trajectories originated from Afghanistan, Bangladesh, Bhutan, China, India, Kazakhstan, Kyrgyzstan, Nepal, Pakistan, Tajikistan, Turkmenistan, Uzbekistan, and Indian Ocean. However, maximum number (~98%) of trajectories originated from Indo-Gangetic plains, northern parts of Bangladesh, and north-eastern states of India. During monsoon, the trajectories originated from Bangladesh, Bhutan, Myanmar, China, India (Indian states of Orissa, Bihar, West Bengal, and the north-eastern states of Assam, Meghalaya, Manipur, Nagaland, Tripura), Nepal, and Indian Ocean. Majority of trajectories (~89%) originated from Bhutan, northern part of Bangladesh and Indian states of Assam, Meghalaya, and Sikkim. During post-monsoon, the trajectories originated from Afghanistan, Bangladesh, Bhutan, Burma, China (Tibet and Xinjiang Province), India, Iran, Kazakhstan, Nepal, Pakistan, Tajikistan, Turkmenistan, Uzbekistan, and Bay of Bengal. Nearly 81% of the trajectories were observed to originate from north-east Indian states, eastern parts of Bihar, northern parts of Bengal, and eastern parts of Nepal and northern parts of Bangladesh. From the back trajectory analysis for winter season, it was observed that the trajectories originated from Afghanistan, Armenia, Azerbaijan, Bangladesh, Bhutan, Burma, China (Tibet and Xinjiang Province), India (North Indian states and Indo-Gangetic Plain), Iran, Iraq, Kazakhstan, Kyrgyzstan, Nepal, Tajikistan, Turkey, Turkmenistan, Uzbekistan, and Bay of Bengal. However, ~96.5% of the trajectories originated from north-east Indian states, northern regions of West Bengal, Uttar Pradesh and Bihar, northern parts of Bangladesh, and eastern Nepal. Overall, from the back trajectory analysis, it can be concluded that nearly 94% of the wind trajectories originate from Indian states of Sikkim, Assam, Meghalaya, northern West Bengal, Northern parts of Indo-Gangetic Plain; northern Bangladesh and Bhutan itself. The general movement of the air parcels during the pre-monsoon is from the west. During monsoon and the post-monsoon, the air parcels are observed to move in from the south. However, during the winter, the air parcels move in from the south-west.

Fig. 4

Gridded back trajectory frequencies with hexagonal binning

Concentration weighted trajectories

The gridded smoothed back trajectory PM2.5 concentrations estimated using the CWT approach are shown in Fig. 5.

Fig. 5

Gridded and smoothed back trajectory concentrations showing mean PM2.5 concentrations using concentration weighted trajectory approach

In all the three years studied, during the pre-monsoon season, the grid cells in Nepal, Indo-Gangetic plains and from some regions of China (Tibet and Xinjiang) were associated with high CWT values up to 40 μg/m3, showing that these regions may be the probable major source regions associate with the PM2.5 concentrations in Bhutan during pre-monsoon. Grid cells in Bhutan during the pre-monsoon were found to have moderate CWT values of up to 30 μg/m3. The high CWT values from southern parts of Nepal and the Indo-Gangetic plains with shorter trajectories indicate stable situation with reduced removal of PM2.5. During monsoon, the grid cells in north-eastern states and West Bengal in India along with Bangladesh, western parts of Myanmar and Bay of Bengal were associated with moderate CWT values of up to 20 μg/m3 of PM2.5. However, grid cells in Bhutan were found to have moderate CWT values of ~25 μg/m3. This shows that majority of the PM2.5 during monsoons majorly contributed to Thimphu from other areas of Bhutan. For post-monsoon, grid cells in eastern Nepal, parts of China (Southern Tibet Province) and north-eastern parts of Uttar Pradesh, Bihar and Uttarakhand and Punjab were associated with high CWT values up to 40 μg/m3 of PM2.5. Grid cells in Bhutan during the post-monsoon season were associated with moderate CWT values of up to 30 μg/m3. In post-monsoons eastern Nepal, parts of China (Southern Tibet Province) and north-eastern parts of Uttar Pradesh, Bihar and Uttarakhand and Punjab can be the major potential source regions for PM2.5 in Thimphu, Bhutan.

During winter, very high CWT values up to 60 μg/m3 were found to be associated with the grid cells in Nepal, north Indian states of Himachal Pradesh, Uttarakhand, Punjab, Haryana, and Delhi, along with northern parts of West Bengal, Uttar Pradesh and Bihar and parts of China (mainly southern Tibet). High CWT values of ~40 μg/m3 were associated with the grid cells in Indo-Gangetic Plain, southern China (Tibet), northern Bangladesh and north-eastern Indian states and West Bengal. Bhutan itself was also found to have grid cells with very high CWT values of up to 60 μg/m3. This goes to show that during winters apart from transboundary source regions, the contribution from local sources was also high.

CWT analysis for 2018 revealed that grid cells lying in Nepal, Southern Tibet and Assam and Arunachal Pradesh, along with grid cells in Bhutan, were associated with very high CWT values up to ~60 μg/m3 and above in winter; grid cells in Bangladesh, Assam, and Sikkim were associated with high CWT values up to 35 μg/m3 in monsoon with moderate CWT values (30 μg/m3) associated with grid cells in Bhutan; grid cells in West Bengal, northern Bangladesh, Arunachal Pradesh, and regions of Tibet adjoining Bhutan, and Bhutan itself were associated with moderate CWT values of a little over 25 μg/m3 in pre-monsoon; and the grids cells in southern part of Tibet, eastern Nepal, north eastern parts of Uttar Pradesh and Bihar in India are associated with high CWT values of ~45 μg/m3 in post-monsoon, while the cells in Bhutan were found to have high CWT values up to 40 μg/m3 form the north-western regions and moderate values up to 30 μg/m3 from eastern regions.

CWT analysis for 2019 showed that grids lying in Nepal, southern parts of Tibet, and northern region of Indian state West Bengal as well as western Bhutan were associated with concentration of PM2.5 as high as 70 μg/m3 in winters; grid cells in Nepal, northern parts of Indian state Uttar Pradesh and Bihar were associated with 45 μg/m3 in pre-monsoon, with grid cells in north-eastern parts of Bhutan associated with low CWT values up to 20 μg/m3; a little over 16 μg/m3 was associated with grid cells located in Bangladesh, Bay of Bengal, eastern parts of West Bengal in India and Meghalaya and Assam along with south western Bhutan for the monsoon season; grids cells located in Indian states of Bihar, West Bengal were associated with high CWT values of 70 μg/m3 PM2.5, and Southern Tibet, Nepal and Indian sates Uttarakhand and Punjab were associated with concentrations up to 40 μg/m3 during post-monsoon. For post-monsoon, grid cells in western Bhutan were found to have high CWT values up to 60 μg/m3, while the eastern parts of Bhutan had CWT values up to 45 μg/m3.

For 2020, for the pre-monsoon season, the CWT analysis revealed that Nepal, southern Tibet, and Indian states including Sikkim, Northern Uttar Pradesh, Bihar, and West Bengal were the major source regions and were associated with high CWT values of up to 40 μg/m3. Grid cells in Northern Bangladesh, parts of western Assam and Meghalaya and central Tibet were associated with moderate CWT values up to 30 μg/m3. However, the grid cells in eastern Bhutan were associated with low CWT values of up to 20 μg/m3, while the cells in the western parts were found to have high CWT values up to 35 μg/m3.

From the CWT analysis, it can be concluded that the regions outside Bhutan with Tibet, north-east India, Indo-Gangetic Plain and Nepal are responsible for significant contributions to the PM2.5 concentrations in Thimphu and can be considered as the major source regions.

Contribution from neighbouring states

The percentage contributions from the neighbouring states of Bhutan were estimated using Eq. (4). The contributions from different countries and the Indian states are shown in Figs. 6 and 7. For simplicity, the contributions below 10% are not represented in this table. In pre-monsoon, the seasonal average of PM2.5 observed in Thimphu was 29.11 μg/m3. India was the highest contributor to the PM2.5 with around ~45% (13.9 μg/m3) with Indian states of Bihar (~12%, 3.8 μg/m3) and Uttar Pradesh (~10%, 3.1 μg/m3) being the major contributors. The next major contributors were Nepal (~17%, 5.5 μg/m3), Bangladesh (~13%, 4.2 μg/m3), and China (Tibet Province) (~12%, 3.7 μg/m3). During monsoon, the average PM2.5 concentration observed during the study period was 15.9 μg/m3. India was the major contributor with contribution of around 43% (7.3 μg/m3) of the PM2.5 concentrations with the states Assam (~12%, ~2 μg/m3), Meghalaya (~14%, 2.2 μg/m3), and West Bengal (~10%, 1.6 μg/m3) being the major contributors. Bangladesh contributed about 31.04% (5.2 μg/m3), and regions of Bhutan around Thimphu contributed around 21% (4.1 μg/m3) of the PM2.5 concentrations. This also shows that there are significant contributions of the regional transport to the total PM2.5 concentrations. During the post-monsoon, the average PM2.5 concentration observed during the study period was 30.4 μg/m3. India contributed around ~42% (12.8 μg/m3) to the total PM2.5 during this season and was the major contributor. In India, the states of West Bengal (~13%, 3.9 μg/m3), Bihar (~10%, 3.1 μg/m3), and Assam (~8%, 2.5 μg/m3) were the major contributors. Apart from India, contribution of Bangladesh (~17%), China (~6%), regions of Bhutan also contributed ~19% of the PM2.5 concentrations in Thimphu. The highest seasonal PM2.5 concentration average of 50 μg/m3 was observed in the winter season. During winter, India was again the major contributor of PM2.5, with contribution of ~45% (22.6 μg/m3). West Bengal and Bihar contributed to ~10% (5.4 μg/m3) and were the major contributing states of India. China contributed ~20% (10.3 μg/m3), Nepal ~13% (6.6 μg/m3), and regions of Bhutan contributed to around ~23% (11.5 μg/m3). From the analysis, it is evident that India is one of the major contributors to the total PM2.5 in Bhutan. The analysis is also in line with the CWT analysis which also suggests that Indian states are mostly associated with the high PM2.5 concentrations. From the analysis, it can be concluded that Thimphu experiences the highest PM2.5 concentrations during winter. Regions of Bhutan also contribute significantly to the total PM2.5 concentrations observed in Thimphu.

Fig. 6

PM2.5 contributions from different countries to total PM2.5 in Thimphu for different seasons

Fig. 7

PM2.5 contributions from Indian states to total PM2.5 in Thimphu for different seasons

Conclusion

From the study, it was concluded that PM2.5 concentrations in Thimphu, Bhutan, exceeded the WHO prescribed limits of 25 μg/m3 on 38.3, 63.4, 39.4% days in 2018 (March–December), 2019 (January–December), and 2020 (January–September), respectively. It was observed that even during the national lockdown, imposed due to the spread of COVID-19, from 11 August 2020 to 21 September 2020, only a marginal reduction of 4% was observed in the PM2.5 concentration in 2020 when compared to the average concentration in 2018 and 2019 combined. This confirms that transboundary sources dominate the contributions to total PM2.5 in Thimphu, Bhutan. The back trajectories allocation shows that for the years 2018 and 2019 south was the dominant sector, and south-west was for 2020. Significant PM2.5 concentrations were attributed to west sector; northwest sector; and southeast and southwest sectors for the years 2018, 2019, and 2020, respectively. From the back trajectory analysis, it was concluded that nearly 94% of the wind trajectories originated from Indian states of Sikkim, Assam, Meghalaya, northern West Bengal, Northern parts of Indo-Gangetic Plain; northern Bangladesh, and Bhutan itself. For all the seasons, most of the trajectories were found to originate from the north-eastern states of India, northern parts of Bangladesh, and eastern parts of Nepal. From the CWT analysis, it was concluded that Tibet, north-east India, Indo-Gangetic Plain, and Nepal are the major contributors of PM2.5 to Thimphu. Also, India was the highest contributor to the PM2.5 contributing around ~ 45% of the average PM2.5 with the mean PM2.5 contribution from being ~ 45% (13.9 μg/m3) in pre-monsoon, 43% (7.3 μg/m3) during monsoon, ~ 42% (12.8 μg/m3) during post-monsoon, and ~ 45% (22.6 μg/m3) during winter season. As hourly ground-based pollutant concentration data are still not available for Bhutan, more ground-based continuous monitoring stations are required for monitoring different pollutant concentrations. Also, daily hospital admission data logging is required to enable estimation of risks associated with PM2.5 in Bhutan. Also, the study’s scope was limited due to the data availability; however, more detailed modelling studies using chemical transport models or dispersion models are required to explore the nature of PM2.5 in Thimphu, Bhutan.