First Measurement of Atmospheric Mercury Species in Qomolangma Nature Preserve, Tibetan Plateau, and Evidence of Transboundary Pollutant Invasion

Located in the world’s ‘Third Pole’ and a remote region connecting the Indian Ocean plate and the Eurasian plate, Qomolangma National Nature Preserve (QNNP) is an ideal region to study the long-range transport of atmospheric pollutants. In this study, gaseous elemental mercury (GEM), gaseous oxidized mercury (GOM) and particle-bound mercury (PBM) were continuously measured during the Indian monsoon transition period in QNNP. A slight increase in GEM concentration was observed from the preceding the Indian Summer Monsoon period (1.31±0.42 ng m) to the Indian Summer Monsoon period (1.44±0.36 ng m), while significant decreases were observed in GOM and PBM concentrations, with concentrations decreasing from 35.2±18.6 to 19.3±10.9 pg m (p<0.001) for GOM and from 30.5±12.5 to 24.9±19.8 pg m (p<0.001) for PBM. A unique daily pattern of GEM concentration in QNNP was observed, with a peak value before sunrise and a low value at noon. Relative to the low GEM concentrations, GOM concentrations (with a mean value of 21.4±13.4 pg m, n=1239) in this region were relatively high compared with the measured values in some other regions of China. A cluster analysis indicated that the air masses transported to QNNP changed significantly at different stages of the monsoon, and the major potential Hg sources shifted from north India and west Nepal to east Nepal and Bangladesh. Because there is a large area covered in glaciers in QNNP, local glacier winds could increase transboundary transport of pollutants and transport polluted air masses to the Tibetan Plateau. The atmospheric Hg concentration in QNNP in the Indian Summer Monsoon period was influenced by transboundary Hg flows. This sets forth the need for a more specific identification of Hg sources impacting QNNP and underscores the importance of international cooperation for global Hg controls.


Introduction
Understanding atmospheric mercury (Hg) concentrations in remote regions is vital to understand the global atmospheric Hg cycling processes (Zhang et al., 2016a;Angot et al., 2016;AMAP/UNEP, 2013). Generally, atmospheric Hg can be divided into three major types: gaseous elemental Hg (GEM), gaseous oxidized Hg (GOM) and particlebound Hg (PBM) (Selin, 2009). Over 95% of atmospheric Hg exists in the form of GEM (Ebinghaus et al., 2002;Huang et al., 2014). Due to its stable chemical properties and long lifetime in the atmosphere (approximately 0.3 to 1 year), GEM can be transported over long distances (Horowitz et al., 2017;Travnikov et al., 2017;Selin, 2009). In contrast, GOM and PBM could deposit quickly from the atmosphere, exposing local environments to significant impacts (Lindberg and Stratton, 1998;Seigneur et al., 2006;Lynam et al., 2014). To understand the global and regional cycling of atmospheric Hg, different Hg monitoring networks and sites have been established in recent decades, such as the Atmospheric Mercury Network (AMNet) (Gay et al., 2013) and Global Mercury Observation System (GMOS), which contains over 40 ground-based monitoring stations distributed in the world (Sprovieri et al., 2016). Generally, atmospheric Hg background concentrations range between 1.5 to 1.7 in the northern hemisphere and 1.1 to 1.3 ng m -3 in the southern hemisphere (Lindberg et al., 2007;Slemr et al., 2015;Venter et al., 2015;Sprovieri et al., 2016). However, existing studies are still far from sufficient to obtain a full understanding of long-range Hg transport because of insufficient monitoring data in remote and less-populated regions (Zhang et al., 2015a;Fu et al., 2012a).
The trans-boundary and long-range transport of pollutants have attracted considerable attentions in the northeastern and southeastern regions of the Tibetan Plateau (Yang et al., 2018;Li et al., 2016;Zhang et al., 2015b;Pokhrel et al., 2016). The transboundary flows of atmospheric pollutants to the Tibetan Plateau have been identified for pollutants such as persistent organic pollutants and black carbon (Yang et al., 2018;Li et al., 2016;Zhang et al., 2015b;Pokhrel et al., 2016). It was reported that smoke from biomass burning in the Indian subcontinent could pass over natural barrier of the Himalaya (Wang et al., 2015;Pokhrel et al., 2016). HCHs, DDTs and PCBs were all found to have their highest concentrations in the southeast Tibetan Plateau during the monsoon season . Similar conditions have also occurred for black carbon . However, studies of the trans-boundary transport of Hg on the Tibetan Plateau are still limited. The existing Hg monitoring data is affected to varying extents by local emission sources (Fu et al., 2012a;Zhang et al., 2015;. Fu et al. (2012a) report that air masses with high Hg concentrations passed over the urban and industrial areas in Western China and Northern India, and influenced the atmospheric Hg concentrations in Waliguan on the northeastern edge of the Tibetan Plateau. At Shangri-La, located on the southeastern edge of the Tibetan Plateau, the atmospheric Hg sources were reported to be Southeast Asia, India and mainland China (Zhang et al., 2015a). Nevertheless, studies are still lacking on trans-boundary transport of Hg in the Qomolangma National Nature Preserve (QNNP), which directly connects the Indian Subcontinent and Eurasia. The detailed pollutant transport pathways and seasonal or daily patterns of atmospheric Hg concentrations in this region are still not clear.
QNNP, located on the southern edge of the Tibetan Plateau, is considered one of the world's cleanest regions (Qiu, 2008). With an average altitude of ~4,500 m a.s.l., QNNP is a remote region with sparse human population and rare industries (Qiu, 2008;Yao et al., 2012b;Li et al., 2016). However, it is surrounded by two large potential pollution sources: the populated and developed eastern China region, which has experienced about 30 years of rapid industrial development, and South Asian developing countries (e.g., India, Nepal, and Bangladesh), which have also been developing rapidly in recent years (Streets et al., 2011;Zhang et al., 2015b;Yang et al., 2018). China and India are reported as the largest coal consumers in the world (BP Statistical Review of World Energy 2018), and coal combustion is the largest source of atmospheric Hg emissions globally, accounting for ~86% of Hg emissions (Chen et al., 2016a). China is predicted to become the largest economy in the world in the next 20-50 years, and India is predicted to catch up with the Euro area before 2030 (Pacyna et al., 2016). The rapidly growing economies have led to rapid increases in energy demands and hence increasing domestic Hg emissions (Pacyna et al., 2016). With the implementation of control strategies, the atmospheric Hg emissions is forecasted to be about 242 tonnes in China in 2020 . However, atmospheric Hg emissions in India are expected to increase to about 540 tonnes Hg by 2020 (Burger Chakraborty et al., 2013). Because QNNP is located on the pathway of air mass transport due to the Indian Summer Monsoon (ISM) , meteorological conditions in QNNP vary significantly during the monsoon transition period (Wang et al., 2001). The monthly average precipitation can range from less than 50 mm in the non-ISM period to 950 mm in the ISM period (Panthi et al., 2015). In addition to the monsoon, the glacial coverage in QNNP is approximately 2,710 km 2 (Nie et al., 2010). Glacier winds could therefore have direct effects on the local pollutant transport because downslope glacier winds can transport polluted air from the upper levels to the land surface . The atmosphere in QNNP is therefore vulnerable to surrounding pollution sources Xu et al., 2009).
To the best of our knowledge, the present work is the first study regarding Hg monitoring and source identification in the QNNP covering both the period preceding the Indian Summer Monsoon (PISM) and during the Indian Summer Monsoon (ISM).
We performed continuous measurements of GEM, GOM and PBM concentrations for 2 weeks during the onset of the monsoon and for 3.5 months during the monsoon itself.
To identify the detailed sources, we combined the real-time Hg monitoring data with a backward trajectory analysis, clustering analysis and potential source contribution function (PSCF) analysis. We further discuss the effects of local glacier winds, caused by the large spatial extent of QNNP glaciers, on the trans-boundary transport of pollutants. This combined monitoring and modeling study could help researchers and government managers to better understand the global Hg cycling processes and potential impacts from the rapidly developing countries in South Asia on the atmospheric Hg concentrations in QNNP.

Atmospheric Hg monitoring site
Atmospheric Hg monitoring was conducted at the "Atmospheric and Environmental Comprehensive Observation and Research Station, Chinese Academy of Sciences on Mt. Qomolangma" (latitude: 28°21'54" N, longitude: 86°56'53" E) in QNNP, at an altitude of 4,276 m a.s.l. (Figure 1). In QNNP, Mt. Qomolangma spreads from east to the west along the border between the Indian subcontinent and the Tibetan Plateau ( Figure 1). Due to its high altitude, QNNP is naturally isolated from the populated regions, and only rare local Hg emission sources have been observed (AMAP/UNEP, 2013). The most populated region near this monitoring site is Tingri County (with a population density of 4 persons per km 2 ), located ~40 km to the southwest of the monitoring site. The average annual temperature is 2.1 °C and the total annual rainfall is 270.5 mm in QNNP (Chen et al., 2016b). QNNP is located along the air mass transport pathway of the ISM , and the meteorological conditions in QNNP have significant variations between the PISM and ISM periods (Wang et al., 2001). During the transition period, the temperature in the Tibetan Plateau and South Asia changes from "southern warm -northern cool" to "northern warm -southern cool" (Wang et al., 2001). This reversal leads to a significant increase of diabatic heating over South Asia and the southern slope of the Tibetan Plateau (Ge et al., 2017), which further affects the wind directions and speeds. Local glacier winds could also affect the transport of air masses in QNNP. Glaciers cover ~2,710 km 2 in QNNP (Nie et al., 2010), and most of the glaciers are located on the northern slope of the mountain (Figure 1) (Bolch et al., 2012). The glacier wind is a continuous downslope wind blowing from glacier surfaces down to the foothills of the mountain throughout the day. Hence, the transport of air masses in this region is a combination of atmospheric circulation (monsoon) and local weather conditions (glacier winds). The structure of the boundary layer over QNNP is also significantly affected by glaciers (Li et al., 2006). The height of the atmospheric boundary layer follows a diurnal profile ranging from ~350 m above ground level during the night to ~2000 m during the day (Li et al., 2006).

GEM, GOM and PBM monitoring
To describe the changes of atmospheric Hg concentrations during the PISM and ISM periods, real-time continuous measurements of GEM, GOM and PBM concentrations were carried out using the Tekran 2537B, 1130 and 1135 instruments (Tekran Inc., Toronto, Canada) from 15 April, 2016 to 14 August, 2016. During the operation of the Tekran instruments, ambient air was introduced into the instrument for 60 minutes through an impactor, a KCL-coated annular denuder, and a Quartz Fiber Filter (QFF).
All the Hg species were converted into Hg(0) and then measured by cold vapor atomic fluorescence spectroscopy (CVAFS). The collected PBM and GOM were desorbed in succession to Hg(0) at the temperature of 800 °C and 500 °C, respectively. Hg-free air was used to flush the 1130 and 1135 systems to introduce the desorbed PBM and GOM into model 2537B for analysis. The GEM was collected at 5-minutes intervals. The sampling inlet was set at ~1.5 m above the instrument platform (shown in Figure S1).
To mitigate the impacts of low atmospheric pressures on the pump's train, a low air sampling rate of 7 L min -1 for the pump model and 0.75 L min -1 (at a standard pressure of 1013 hPa and temperature of 273.14 K) for model 2537B was applied Zhang et al., 2015a;Zhang et al., 2016a). The Tekran 2537B analyzer was calibrated automatically using the internal Hg permeation source inside the instrument every 23 h, and the internal source was calibrated before and after the monitoring by an external Hg source using a syringe. The Tekran ambient Hg analyzer has been described in more details in the previous publications (Landis et al., 2002;Rutter et al., 2008;de Foy et al., 2016). Recent studies have suggested that there may be a low bias of GOM and PBM concentrations for small sample loads of Hg(e.g. less than 10 pg) (Slemr et al., 2016;Ambrose, 2017). Hence, the monitoring data with GOM or PBM concentrations below 23.8 pg m -3 were recalculated by the method of Slemr et al. (2016). The updated GOM concentrations increased slightly from 21.3±13.5 pg m -3 to 21.4±13.4 pg m -3 and from 25.5±19.2 pg m -3 to 25.6±19.1 pg m -3 for PBM.

Meteorological data
Throughout the sampling period, the meteorological information was recorded using the Vantage Pro2 weather station (Davis Instruments, USA) with a 5-minute resolution.
The monitored parameters included the temperature (with a precision of 0.1°C), relative humidity (with a precision of 1%), wind speed (with a precision of 0.1 m s -1 ), wind direction (with a precision of 1°), air pressure (with a precision of 0.1 hPa), solar radiation (with a precision of 1 W m -2 ) and UV index (with a precision of 0.1 MEDs).
The snow cover data was obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument on board the Terra and Aqua satellites (MOD10A1, Hall et al., 2010) with a daily 0.05° resolution.

Backward trajectory simulation
To identify the atmospheric Hg sources, the Hybrid Single-Particle Lagrangian Integrated Trajectory (HYSPLIT) model was applied to perform a backward trajectory simulation (Stein et al., 2015;Chai et al., 2016;Chai et al., 2017;Hurst and Davis, 2017).
The HYSPLIT model, known as a complete and mature system for modeling air parcel trajectories of complex pollutant dispersion and deposition, was developed by the US National Oceanic and Atmospheric Administration (NOAA). Global Data Assimilation System (GDAS) data with 1°×1° latitude and longitude horizontal spatial resolution and 23 vertical levels at 6-hour intervals was used for the backward trajectory simulation. All the trajectory arrival heights were set to 1000 m above ground level.
Every backward trajectory was simulated for 72 hours in 6-hour intervals, and the air mass transport regions covered China, Nepal, India, Pakistan and majority of west Asia.
Backward trajectories during the whole monitoring period were calculated, and cluster analysis was carried out to identify the Hg transport pathways. The cluster statistics summarize the percentage of back trajectories in each cluster, and the average GEM concentrations are linked with each cluster. The clustering algorithm utilized in this study is based on Ward's hierarchical method (Ward Jr, 1963), which minimizes angular distances between corresponding coordinates of the individual trajectories. By averaging similar or identical pathways from existing air mass pathways to the receptor site, clusters can help identify the mean transport pathways of air masses and provide the primary directions of pollutants transported to the measurement site.
The Potential Source Contribution Function (PSCF) model is a hybrid receptor model using the calculated backward trajectories to estimate the contributions of different emission sources in upwind regions and has been applied in many previous studies (Kaiser et al., 2007;Fu et al., 2012b;Kim et al., 2005;Zhang et al., 2013). The PSCF calculation is made based on counting the trajectory segments that terminate within each cell to determine the values for the grid cells in the study domain (Ashbaugh et al., 1985). In this study, the PSCF model was used to identify the possible sources of atmospheric GEM. The study domain was separated as i × j cells. Then, the PSCF value for the ij th cell is defined as follows: where N ij is the total number of endpoints that fall into ij th cell during the whole simulation period, and M ij is the number of endpoints for the same cell that correspond to GEM concentrations higher than a set criterion. In this study, PSCF values were calculated based on the average GEM concentration during the whole sampling campaign. The PSCF value stands for the conditional probability that the GEM concentration at the measurement site is larger than the average GEM concentration if the parcel passes through the ij th cell before it reaches the measurement site.
To account for and reduce the uncertainty due to low values of N ij , the PSCF values were scaled by an arbitrary weighting function W ij (Polissar et al., 1999). While the total number of the endpoints in a cell (N ij ) is less than ~three times the average value of the end points for each cell, the weighting function will decrease the PSCF values.
In this study, W ij was set using the following piecewise function: , %& =

Comparisons of atmospheric Hg concentrations between PISM and ISM
The GEM, GOM and PBM concentrations at the sampling site were 1.42±0.37 ng m -3 (n=15180), 21.4±13.4 pg m -3 (n=1239) and 25.6±19.1 pg m -3 (n=1237), respectively, during the whole study period ( Figure 2 and Table 1). GEM accounted for over 95% of all the atmospheric Hg species. Figure S2 shows    In QNNP, with the wide distribution of glaciers, glacier winds could bring the upper air masses to the land surface layer , which could further strengthen the subsidence movement. Low wet deposition rate of GOM caused by the rare precipitation in QNNP (~270mm) (Chen et al., 2016c) could be another reason for the high GOM concentrations (Prestbo and Gay, 2009).
The increases of GEM concentrations during the ISM period could indicate the impacts of trans-boundary transport, which has been confirmed by previous studies (Fu et al., 2012a;Zhang et al., 2016a). The deposition of GEM from the atmosphere to the land surface is difficult, and GEM has a much longer residence time than the other Hg species (Horowitz et al., 2017;Travnikov et al., 2017;Selin, 2009). At Ailaoshan in Yunnan province (Zhang et al., 2016a), a higher TGM concentration during the ISM period (2.22±0.58 ng m -3 ) than the PISM period (1.99±0.66 ng m -3 ) was also observed.
The TGM concentration during the ISM period (2.00±0.77 ng m -3 ) was also higher than that during the PISM period (1.83±0.78 ng m -3 ) at Waliguan station in the northeastern Tibetan Plateau (Fu et al., 2012a). In contrast to GEM, the GOM and PBM levels during the ISM period were lower than the monitored values during the PISM period ( Figure   S2 and Table 2). In previous studies, the PBM concentration in the Kathmandu Valley was lower during the monsoon period (with a value of 120.5±105.9 pg m -3 ) than the pre-monsoon (with a value of 1855.4±780.8 pg m -3 ) and post-monsoon period (with a value of 237.6±199.4 pg m -3 ) (Guo et al., 2017). In India, PBM concentrations during the monsoon period (with a value of 158±34 pg m -3 ) were lower than those in the nonmonsoon season (with a value of 231±51 pg m -3 ) (Das et al., 2016). This fact could be possibly attributed to precipitation increases brought by the monsoon, which further causes wet depositions of PBM from atmosphere. During the ISM period, the precipitation could increase by up to 25% in the South Asia and Tibetan Plateau (Ji et al., 2011).

Diurnal variation of atmospheric Hg species in QNNP
During the PISM period, all the atmospheric Hg species showed clear diurnal patterns (Figure 3 and Figure S3). For GEM, the minimum concentrations usually  (Han et al., 2009;Tie et al., 2007;Quan et al., 2013). With a large glacier coverage (~2,710 km 2 ), the structure of the boundary layer over QNNP was significantly affected by glacier winds (Li et al., 2006). The local PBL may be subject to impacts from the glacier-covered environment and have a significant diurnal variation. The height of the atmospheric boundary layer could vary significantly from ~350 m above ground level to ~2000 m in one day (Li et al., 2006). Following sunrise, with the strengthening of the glacier wind, a strong convection current starts to grow in the troposphere, and the stock of GEM in the nearground atmosphere is depleted quickly, leading to the quick decrease in concentrations.
In contrast, after sunset, with the weakening of the glacier wind, the nocturnal stable boundary layer takes a dominate position controlling the surface layer, and its height is relatively low (Li et al., 2006), which could lead to increases in GEM concentrations.
Comparing the diurnal variations between the PISM and ISM period, the atmospheric and approximately 60% of the Hg deposited to snow cover could eventually be reemitted to the air. A shorter reservoir lifetime for deposited Hg in snowpack was also reported when temperature rises (Faïn et al., 2007). With the increase of ambient temperature and radiation from April to August, the reemission of GEM from the glaciers could increase as well. As the snow coverage in the QNNP decreased significantly from the PISM to the ISM period ( Figure S4), some of the released Hg may become a source of new GEM from the initial ISM to the final stage of the ISM period. More Hg(0) could be released due to the higher temperature and stronger radiation in the afternoon. However, some other factors such as changes in the PBL heights and in wind directions could also be partly responsible for the diurnal variations of GEM concentrations (Horowitz et al., 2017;Travnikov et al., 2017;Selin, 2009;Li et al., 2006). Relatively low GEM concentrations (<1.5 ng m -3 ) were observed in most of the samples (80.0%) of air masses in the predominant Hg-transport direction (from southwest to west) during the PISM period, which is due to the control of westerlies.

Wind direction dependence of Hg concentrations
With high wind speeds (Table 1) and coming from Central Asia, the westerlies are the predominant wind containing low pollutant levels that spread in the QNNP during the PISM period (Kotlia et al., 2015). Relatively high GEM concentrations (>1.5 ng m -3 ) were found in 92.4% of the samples for the predominant Hg direction during the ISM period under the control of the monsoon (Kotlia et al., 2015), which might indicate that the transported air masses are coming from polluted regions. GOM and PBM had similar patterns under the control of the westerlies and monsoon during the PISM and ISM period, respectively.

Air mass back trajectories analysis
To further quantify the contributions of different sources to GEM concentrations, an air mass back trajectory simulation and trajectory cluster analyses were applied in this study. Figure (Xiao et al., 2015;Chen et al., 2015;Rhode et al., 2007;Huang et al., 2016).
Cluster 1 (17%) and cluster 3 (60%) represent the pollutant coming from Nepal, and the trajectory is relatively short. During the ISM2 period, all the clusters originated from or passed through central Asia, northern India and northwestern Nepal (Figure 5c).
The clusters were similar to most of the clusters during PISM period; however, the GEM concentrations in these clusters were higher than those during the PISM period, which might be caused by the large Hg emissions from frequent fires in the source region during ISM2 ) ( Figure S5). During the ISM3 period ( Figure   5d) India are related to the rapidly growing economy and increasing usage of fossil fuels (Sharma, 2003). Considering the heavy air pollutions in Nepal (Rupakheti et al., 2017;Forouzanfar et al., 2015) and in Bangladesh (Islam et al., 2015;Rahman et al., 2018;Rana et al., 2016;Mondol et al., 2014), Nepal and Bangladesh might be underestimated Hg source regions in the modeling and should be taken into consideration in further study.
Under the control of the ISM during the ISM2 period, the high PBM concentration may be related to the biomass burning in the source region. According to the PSCF analysis, northern India and Nepal are the potential source regions during the ISM2 period. The source identification by back trajectory simulation and trajectory cluster analyses also indicated that northern India and Nepal are in the air mass transport trajectory that would transport Hg to QNNP. Finley et al. (2009) reported that PBM concentrations could be associated with Hg emissions from wildfire events. One possible cause of the observed high PBM concentration is the frequent fire events that occurred during the ISM2 period in the air masses trajectory. Figure S4 shows the fire hotspots observed by MODIS from April to August 2016. During the ISM2 period, frequent fire hotspots were identified in the source region, and large amounts of PBM may have been released into the atmosphere from biomass burning ).
The transport of those air masses with enriched PBM was controlled by the ISM and intensified by glacier winds. The transport of polluted air to QNNP resulted in the outburst of PBM concentration during the ISM2 period. During the PISM period, although the number of fire hotspots was much higher, most of the fire hotspots locations were not in the potential source region (Figure 6a and Figure S4), resulting in the low PBM concentration observed.

Implications from this study
At a high altitude and located in the deep southern Tibetan Plateau, QNNP is isolated from anthropogenic perturbations and industrial activities, and this area was thought to be shielded from pollutant inputs from South Asia. However, our results show that the Hg concentration in this region is not as low as previously expected. During the whole monitoring period, the highest GEM concentration reached 3.74 ng m -3 (with trajectories passing through the north of India), ~2.5 times higher than the average concentration in the Northern Hemisphere (~1.5 to 1.7 ng m -3 ) (Lindberg et al., 2007;Slemr et al., 2015;Venter et al., 2015). The average GEM concentration in the middle stage of the ISM was 1.56 ng m -3 , which is inside the average range of observed We now recognize that trans-boundary transportation is an important mechanism that can influence Hg distribution in this region. In particular, the air masses transported to QNNP might be primary under the control of mesoscale ISM drivers and intensified by regional glacier winds (Figure 7). From the PISM to ISM periods, the warm center gradually shifts northwestward from low latitudes to the QNNP (Wang et al., 2001;Ge et al., 2017), and the South Asian High moves onto the Tibetan Plateau and maintains a strong upper-level divergence and upward motion. The upward motion makes the air masses cross the high-altitude Himalayan Mountains and move to mainland China (Xu et al., 2009;Bonasoni et al., 2010). During the ISM period, the transboundary transport of atmospheric Hg is strengthened by both monsoon and glacial winds. However, this effect seems to be weaker during the PISM period. The transboundary-transported air masses can be pumped down right after crossing Mt. Qomolangma due to the control of the regionally unique wind transportation mode, the glacier wind. Hence, in addition to the monsoon, the trans-boundary transport of Hg could also be intensified by regional glacier winds, leading to the increases of atmospheric Hg in this region. As showed in other studies in the northern or eastern Tibetan Plateau, the glacier wind can pump down air masses from upper level to the surface in QNNP . The pump movement is remarkably efficient at transporting air masses (Zhu et al., 2006), and could bring significant amount of pollutants to QNNP.
In 2013