Emission-dominated gas exchange of elemental mercury vapor over natural surfaces in China 1

Abstract. Mercury (Hg) emission from natural surfaces plays an important role in global Hg cycling. The present estimate of global natural emission has large uncertainty and remains unverified against field data, particularly for terrestrial surfaces. In this study, a mechanistic model is developed for estimating the emission of elemental mercury vapor (Hg0) from natural surfaces in China. The development implements recent advancements in the understanding of air–soil and air–foliage exchange of Hg0 and redox chemistry in soil and on surfaces, incorporates the effects of soil characteristics and land use changes by agricultural activities, and is examined through a systematic set of sensitivity simulations. Using the model, the net exchange of Hg0 between the atmosphere and natural surfaces of mainland China is estimated to be 465.1 Mg yr−1, including 565.5 Mg yr−1 from soil surfaces, 9.0 Mg yr−1 from water bodies, and −100.4 Mg yr−1 from vegetation. The air–surface exchange is strongly dependent on the land use and meteorology, with 9 % of net emission from forest ecosystems; 50 % from shrubland, savanna, and grassland; 33 % from cropland; and 8 % from other land uses. Given the large agricultural land area in China, farming activities play an important role on the air–surface exchange over farmland. Particularly, rice field shift from a net sink (3.3 Mg uptake) during April–October (rice planting) to a net source when the farmland is not flooded (November–March). Summing up the emission from each land use, more than half of the total emission occurs in summer (51 %), followed by spring (28 %), autumn (13 %), and winter (8 %). Model verification is accomplished using observational data of air–soil/air–water fluxes and Hg deposition through litterfall for forest ecosystems in China and Monte Carlo simulations. In contrast to the earlier estimate by Shetty et al. (2008) that reported large emission from vegetative surfaces using an evapotranspiration approach, the estimate in this study shows natural emissions are primarily from grassland and dry cropland. Such an emission pattern may alter the current understanding of Hg emission outflow from China as reported by Lin et al. (2010b) because a substantial natural Hg emission occurs in West China.


Introduction
Accurate inventories of mercury (Hg) emission are the foundation for assessing Hg global biogeochemical cycling (Selin, 2009;Streets et al., 2009Streets et al., , 2011)).The inventories of Hg emission include the emission from anthropogenic sources and the so-called "natural" emission that includes the primary natural release (i.e., from geogenic activities) and the re-emission of legacy Hg stored in the terrestrial and water surfaces.Hg emission from anthropogenic sources has been quantified and updated with reasonable consistency since the 1990s (Streets et al., 2009(Streets et al., , 2011;;Zhang et al., 2015;Zhang et al., 2016).In particular, the inclusion of the release from commercial products and modifications of Hg emission speciation profiles corresponding to the implementation and upgrade of air pollution control technologies have substantially reduced the uncertainty of anthropogenic Hg emission esti-Published by Copernicus Publications on behalf of the European Geosciences Union.
mates (Horowitz et al., 2014;Zhang et al., 2016).In contrast, estimates of natural Hg emission are poorly constrained and have large uncertainties (±2000 Mg yr −1 ), limiting the understanding of global and regional Hg cycling budgets (Pirrone et al., 2010;X. Wang et al., 2014;Song et al., 2015).In light of the control of anthropogenic Hg emission by the legally binding Minamata Convention (Pacyna et al., 2016), a better quantification of natural Hg emission is critical in evaluating the effectiveness of policy actions (Selin, 2009;Pirrone et al., 2010;Song et al., 2015).
One of the greatest challenges in predicting natural Hg emissions is the limited understanding in the air-surface Hg 0 exchange processes among multiple environmental compartments, such as those in a terrestrial vegetative ecosystem.Estimates from regression-based models derived from the correlations between Hg flux and environmental parameters (temperature, solar radiation, etc.) may not be representative (Xu et al., 1999;Bash et al., 2004;Lin et al., 2005;Gbor et al., 2006;Shetty et al., 2008;Selin et al., 2008;Smith-Downey et al., 2010) because the relationships between measured fluxed and environmental factors are based on limited field data that are site specific, which cannot account for the heterogeneity in soil properties and vegetative coverages.Although the bidirectional resistance schemes describing Hg 0 gas exchange may be appropriate (Bash, 2010;X. Wang et al., 2014;Wright and Zhang, 2015), they are limited by the availability of required soil property data and other physicochemical parameters such as Hg II reduction kinetics and characteristics of interfacial exchanges (Bash, 2010;X. Wang et al., 2014), leading to inconsistences between simulated and measured fluxes.Other challenges, including acquiring and integrating the meteorology, land use, and soil property data (Hg content and other characteristics) in a model domain covering China, call for further model development to estimate natural Hg emission (X.Wang et al., 2014).
Recent advances in the understanding of Hg II reduction provide new opportunities to build a more physically robust air-soil exchange scheme.These include constrained 10 −11 to 10 −10 s −1 pseudo-first-order rate constant of Hg II reduction in soil (Scholtz et al., 2003;Qureshi et al., 2011) and 0.2-1.0h −1 rate constant in natural water (O'Driscoll et al., 2006;Qureshi et al., 2010).In these reactions, the UV band of actinic light has been shown to be the primary driver for Hg II photoreduction in soils and water bodies (Moore and Carpi, 2005;Si and Ariya, 2011), and the role of functional substructures (e.g., -C(O)OH, -SH, -OH) of DOM (dissolved organic matter) in the reduction has been determined by kinetic studies using model compounds (He et al., 2012;Si andAriya, 2011, 2015).The presence of dissolved O 2 has been shown to inhibit most aqueous Hg II reduction but not the photoreduction of Hg II bound to R-SH (Si and Ariya, 2011).In dry soil, the first-order rate constants of Hg II photoreduction are 0.007-0.028h −1 for HgCl 2 coated over sand and 0.003-0.006h −1 for Hg II in a natural soil (Quinones and Carpi, 2011).In the absence of light, Hg II reduction in soil is also observed at a rate of 0.001-0.002h −1 at 293 K (Pannu, 2012).
Intercontinental transport of Hg emission in China has been suggested to enhance Hg deposition in North America (Jaffe et al., 2005;Strode et al., 2008;Lin et al., 2010b;Chen et al., 2014).However, the natural Hg emission inventory used in earlier modeling studies may not be representative.S. Wang et al. (2014) suggested the Hg emissions in China from natural and anthropogenic sources are comparable.Shetty et al. (2008) estimated the natural Hg emission in China to be 462 Mg yr −1 using an outdated model scheme, which was applied for assessing regional Hg budgets in Est Asia (Lin et al., 2010b;Zhu et al., 2015).The large uncertainty associated with the earlier estimate warrants a reassessment of these earlier modeling efforts.In addition, the data of soil Hg concentration used in Shetty et al. (2008) are more than 3 decades old and require updates to appropriately represent spatial distribution of soil Hg contents that have been modified due to the rapid industrial development and urbanization occurring in China since 1980s.This data deficiency has been addressed by the National Multi-Purpose Regional Geochemical Survey (NMPRGS) of China, which was completed in 2014 (Li et al., 2014).This database provides soil Hg content for the agricultural and industrialized regions at a resolution of 4 km, which may substantially reduce uncertainty caused by data deficiency.In addition, the datasets of terrestrial flux in mainland China recently reported in the literature allow verification of model results and optimization of model schemes.These field data of Hg 0 air-surface exchange in China document the flux characteristics over different land uses including urban-rural-remote differences and effects of crop rotation over agricultural lands (Fu et al., 2008(Fu et al., , 2010(Fu et al., , 2012(Fu et al., , 2013a;;Zhu et al., 2011;Zhu et al., 2013;Zhu et al., 2016;Sommar et al., 2013aSommar et al., , b, 2016a)).
Given the scientific advancements and new data availability discussed above, this work develops a state-of-the-science mechanistic model for estimating the natural Hg emission in China.For the first time, the simulated natural emission flux is verified with field measurements over different land surfaces in a modeling effort.The spatial, temporal, and seasonal characteristics of the model-estimated emissions over soil, vegetative surface, and water are presented and compared to the earlier estimates performed by Shetty et al. (2008).The implications of the new estimate are discussed in terms of chemical transport of Hg emission in China and the need for a re-assessment of mercury emission outflow from China.

Model description
Compared to the earlier mechanistic schemes (X.Wang et al., 2014;Bash, 2010;Scholtz et al., 2003;Zhang et al., 2012b), this model (1) builds a new scheme for estimating the airsoil flux based on the reduction pathways of reactive Hg in soil identified in the literature, (2) develops a scheme for the Hg flux exchange over rice paddy, which is an important land use feature in China, and (3) updates the scheme for the airsnow interface and chemical parameters for air-foliage flux (Table 1).

Reduction of Hg II in soil
Based on Hg II reduction mechanisms proposed in peerreviewed literature (Moore and Carpi, 2005;Quinones and Carpi, 2011;Si and Ariya, 2011;Pannu, 2012), a new model describing Hg 0 formation from Hg II reduction in soil is developed using three reaction pathways: (1) photoreduction of Hg II in soil pore water (Hg 0 1 ), (2) photoreduction of Hg II associated with soil particles (Hg 0 2 ), and (3) Hg II reduction through nonphotochemical pathways (Hg 0 3 ).The production of Hg 0 in surface soil is calculated as where K 1 is the photoreduction rate constant of Hg II in soil pore water (a comprehensive parameter list with units is shown in Table 1), K 2 is the photoreduction rate constant of Hg II associated with soil particles, K 3 is the Hg II reduction rate constant in soil through nonphotochemical pathways, Hg s,m is the Hg II pool in soil pore water, Hg p,m is the labile Hg II pool available for reduction on soil particles, and Hg t,m is the total reducible Hg II pool in soil.Based on the Arrhenius equation, K 1 and K 2 are parameterized as a function of solar irradiance and soil temperature, and K 3 is parameterized as a function of soil temperature and soil moisture: where k 1 and k 2 are the photoreduction rate constants at the reference soil temperature (T f ; Table 1).k 3 is the darkreduction rate constant at the reference soil temperature and soil moisture (θ f ; Table 1).R and R i are total solar irradiance in the soil profile and under the canopy, respectively.γ is the ratio of UV over total solar irradiance.An empirical rule suggests that a 10 • C temperature increase doubles the rates for chemical reactions near room temperature (Kissinger, 1957;Hood et al., 1975), which has been shown to be applicable to Hg II reduction in boreal soil (Pannu et al., 2014).In addition, Hg emission flux from soil substrate doubles at ∼ 25 % increase of soil moisture content (Lin et al., 2010a).Based on these observations, Eqs. ( 4)-( 6) can be simplified as R and R i in Eqs. ( 7)-( 8) are calculated based on the Beer-Lambert law: where R 0 is solar irradiance above the canopy, K is the canopy light extinction coefficient, LAI is leaf area index, k r is the light extinction coefficient in soil, and l is the depth of surface soil.Hg s,m , Hg p,m , and Hg t,m are calculated based on Lyon et al. (1997): where Hg t is the total Hg II concentration in soil, BD is the soil bulk density, θ is the soil moisture, V is the soil volume, and ϕ is the ratio of reducible Hg over total Hg in soil.kd is the soil-water partition coefficient and calculated following Lee et al. (2001) and Sauve et al. (2000): where f oc is the fraction of organic carbon in surface soil.The values r, s, and t are regression coefficients.Following Obrist et al. (2014), we assume that the Hg 0 emission from soil is controlled by diffusion after Hg II reduction.Basing on the Fick's first law, the observed air-soil flux exchange can be calculated as where D soil and D 0 are the diffusivities of Hg 0 vapor in soil and ambient air, GEM is the atmospheric Hg 0 concentration, Photoreduction rate constant in soil solution 6.0 Reference soil temperature 32 • C (Eq. 8), 20 and σ is the soil porosity.Hence, during a time period t, the soil Hg 0 vapor compensation point used in the bidirectional resistance model can be derived as

Updates for air-terrestrial exchanges
Extending from the two categories (vegetated canopy and bare land) from our earlier work (X.Wang et al., 2014, the terrestrial system is divided into four categories: vegetated surface with unsaturated soil moisture (forest, grassland, shrubland, etc.), vegetated surface with saturated soil (i.e., rice paddy), barren or sparsely vegetated land, and snow/ice surface.The governing equation for calculating Hg 0 airsurface exchange over vegetated surfaces is where R a is the aerodynamic resistance, R b is the quasilaminar sublayer resistance, C atm is the atmospheric Hg concentration, and χ cnp is the overall compensation point.For the canopy biomes with unsaturated soil, χ cnp is parameterized as in X. Wang et al. (2014): where χ c is the cuticular compensation point, χ s is the stomatal compensation point, R c is the cuticular resistance, R s is the stomatal resistance, R g is the soil diffusion resistance, and R ac is the in-canopy aerodynamic resistance (Table 1 in detail).While for vegetated surface with saturate soil, χ cnp is parameterized as the following: where χ w is the air-water compensation point and R w is the diffusional resistance on water surface.The governing equation for air-surface exchange in barren or sparsely vegetated land, and over snow/ice surface is For bare land, R g is calculated following Zhang and Lindberg (1999): For snow/ice surface, R g is evaluated following Zhang et al. (2012b): where R g(SO 2 ) and R g(O 3 ) are the diffusion resistances of SO 2 and O 3 , α Hg 0 is the Hg 0 scaling factor based on SO 2 , and β Hg 0 is the Hg 0 scaling factor based on O 3 .The formulation of R g(SO 2 ) and R g(O 3 ) has been described previously (Zhang et al., 2003;X. Wang et al., 2014).The χ g for the air-snow interface is assumed to be 3 ng m −3 based on field measurements at air-snow interface (Mann et al., 2015;Lalonde et al., 2003;Fain et al., 2007;Maxwell et al., 2013).

Model configuration and data
The modeling domain is in Lambert conformal projection, with 223 × 149 grid cells at a 36 km spatial resolution.The modeling period is 1 year (2013).Hourly meteorological data are prepared using the Weather Research and Forecasting (WRF) model version 3.7.Sensitivity analysis in X. Wang et al. (2014) showed that accurate representation of environmental parameters (temperature, solar irradiance, etc.) greatly improves the flux estimate.To obtain the best physics and dynamics options of WRF for the China domain, a L 25 (5 6 ) orthogonal design of experiments is utilized (Supplement, Table S1).The best combination of meteorological physics options is selected based on model evaluation metrics R (correlation coefficient) and RMSE (root mean square error) between simulated outputs of each combination of options and observed values in 750 meteorological stations.The selected physics options are Thompson (microphysics options), Betts-Miller-Janjic (cumulus parameterization options), RRTMG (radiation physics options), and BouLac (PBL physics options) based on the results of meteorological model performance evaluation (Fig. S1 in the Supplement).
The datasets for surface soil properties (1 km spatial resolution) containing organic matter contents, pH, bulk density, and porosity are adopted from Shangguan et al. (2013).The land cover data (1 km spatial resolution) are obtained from Ran et al. (2012).The land cover map represents the best available datasets and follows the IGBP (International Geosphere-Biosphere Programme) classification system (Fig. S2).The ratio of rice planting fields in China during each month are classified following the method used in Liu et al. (2013).The rice planting seasons are April to October in South China (including double rice planting), and May to October in Northeast China (single rice planting).The LAI data, also with a 1 km spatial resolution, are adopted from Yuan et al. (2011).The soil Hg content information utilized by Shetty et al. (2008) is updated and greatly expanded with the comprehensive NMPRGS dataset (Li et al., 2014).These high resolution datasets were regridded into the domain specification for each land use using the spatial tools in ArcGIS 10.1.The soil Hg content in 0-20 cm surface soil varies with land use types, containing mean concentrations of 119-211, 61-197, 80-82, 80-82, and 31-162 ng g −1 of Hg for for- est ecosystems, shrubland, savanna/grassland, cropland, and other land uses, respectively (Table 1 and Fig. 1).Though the mean Hg concentration in 0-20 cm soil profile could somewhat underestimate Hg concentration in the top soil layer, the dataset is the best available one describing the soil Hg concentration in China.Datasets of Hg concentration in the top soil layer (e.g., 0-5 cm depth) are recommended when they become available.
In the simulation, the validated Hg 0 concentration retrieved from the output of the Hg extension of Community Multi-scale Air Quality modeling system (CMAQ-Hg) version 4.7 for the same modeling period is applied to represent the ambient air concentration of Hg 0 (Lin et al., 2010b).The simulation does not incorporate the feedback of the airsurface exchange to the air concentration because the feedback of the air-surface exchange to the air concentration does not significantly modify the atmospheric Hg 0 concentration in each grid cell, and the typical variation range of ambient Hg 0 concentration is not a sensitivity parameter for flux change (X.Wang et al., 2014)

Evaluation for soil Hg reduction scheme
Values of all model parameters used in the simulation are showed in Table 1.The value of k 1 is assumed to be 6 × 10 −9 m 2 W −1 s −1 based on the relationship between irradiance intensity and apparent photoreduction rate constant in aerated solution observed by Si and Ariya (2011).Considering the 2 mm maximum photolysis penetration depth in soil (Hebert and Miller, 1990), the measured rate constant in soil particles (depth = 2.07 mm) from Quinones and Carpi (2011) is 2 × 10 −3 m 2 W −1 h −1 (k 2 ) with respect to the pool of labile Hg II available for reduction.The value of k 3 is assumed to be 1.0 × 10 −3 h −1 based on the average rate constants for dark (thermal) reduction (Pannu, 2012).The mean ratio of reducible Hg in soil is assumed to be 0.03 for the soil with vegetation based on measurements from Pannu (2012).No data are available for the bare soil.Data from Lindberg et al. (1999) imply that Hg-enriched desert soil (1400-5000 ng g −1 total Hg) produces a nominal Hg 0 efflux in the range from 40 to 60 ng m −2 h −1 .Derived from back calculation taking pore diffusion into account, the fraction of reducible Hg is predicted at least 10 times lower (≤ 0.003) than that in the soil with vegetation.
Sensitivity analyses using a box model for a typical forest soil are performed to gauge the selected rate coefficients, and the results are shown in Fig. 2 the Hg 0 vapor concentration in soil porous media is estimated to be 4.5 ng m −3 , comparable to the measured concentration (4.1 ± 2.0 ng m −3 ) in the surface forest floor (Moore and Castro, 2012).This suggests that the selected empirical constants appropriately represent typical environmental conditions.Generally, the range of Hg 0 vapor in all simulations is 1.5-6.7 ng m −3 .Less than 0.1 % Hg 0 vapor is from photoreduction in soil solution as the Hg pool in soil solution is small (≤ 0.1 % of total Hg concentration).A ∼ 16 % fraction of the pore Hg 0 concentration derives from thermal Hg II reduction, contributing to 0.5-1 ng m −3 of Hg 0 vapor present in soil gas.Hg 0 concentrations in pore gas are typically lower than the 1-2 ng m −3 atmospheric Hg concentration in background forest at night (Carpi and Lindberg, 1998;Ericksen and Gustin, 2006;Kuiken et al., 2008a, b;Obrist et al., 2014;Fu et al., 2015), suggesting that the forest floor acts as a Hg 0 sink during nighttime.This is consistent with the sign of nocturnal fluxes observed over forest floor (Carpi and Lindberg, 1998;Ericksen et al., 2006;Kuiken et al., 2008a, b).Moore and Carpi (2005) reported that the Hg flux under sun-lit condition is 3-5 times higher than the value observed in the dark.The developed model is capable of producing results consistent with the observation that photoreduction on soil particles dominates the formation of Hg 0 vapor.
Figure 3 illustrates the model response to the model variables at the two experimental levels in Table 2.The twolevel factorial design of experiments is meant to gauge the extreme variation of flux caused by the possible range of all parameters.This method is statistically robust, and therefore the synergistic and antagonistic interactions among model parameters can be estimated with indications of statistical significance.On average, increasing soil bulk density from 0.1 to 1.5 g cm −3 , Hg content from 10 to 400 ng g −1 , soil temperature from 0 to 30 • C, and solar radiation from 0 to 1000 W m −2 will significantly enhance the flux by 20-30 ng m −2 h −1 .Additional 18-20 ng m −2 h −1 synergistic ef-  fects from the combination of above parameters are also predicted.Filed measurements suggest the combined effects of soil Hg content (from 60 to 590 ng g −1 ) and soil temperature (from 5 to 30 • C) enhance the flux by ∼ 40 ng m −2 h −1 (Fu et al., 2012(Fu et al., , 2008)).On the other hand, increasing leaf area index (LAI) from 0 to 7 m 2 m −2 reduces the flux by 19 ng m −2 h −1 .Furthermore, LAI could offset the positive effects from bulk density, soil Hg concentration, and solar radiation above canopy, leading to an additional −19 to −16 ng m −2 h −1 decrease, indicating that the canopy shading substantially constrains soil Hg evasion, consistent with the data that the shading could decrease 70-90 % fluxes compared to nonshaded soils in filed measurements (Carpi and Lindberg, 1998;Zhang et al., 2001;Choi and Holsen, 2009).
In reality, since the actual variation of the parameters is much smaller than the possible range, the flux change will also be much milder.To illustrate this, we run the model using the center values of selected parameters (i.e., showing the model results by running the model at half of the experimental level).Using the center values of soil Hg content, LAI, soil bulk density, solar radiation, and soil temperature in Table 2 (close to the environmental parameters in a typical forest ecosystem), the air-soil flux is 4.5 ng m −2 h −1 , similar to the measured fluxes (0.5-9.3 ng m −2 h −1 ) in forest ecosystems of China (Fu et al., 2012(Fu et al., , 2008).
In the new scheme, the soil organic matter is not incorporated into either K 2 or K 3 , in accordance with the findings of Pannu (2012).While Hg 0 evasion from substrates coated with HgCl 2 and humic matter is inversely correlated with humic matter content both in the dark and under irradiation, the inhibitory effect from humic matter is not linear to its content (Mauclair et al., 2008).For instance, relatively small differences are observed at humic matter content > 1 % (Mauclair et al., 2008).In addition, the effect of soil organic matter type has not been comprehensively investigated (Zhang and Lindberg, 1999;Bash et al., 2007).Further studies that quantify the reduction rate constants associated with different types of soil organic matters (or species) under light, as well as field flux data that relate the observed flux intensity to a given type Table 3. Mean annual air-surface fluxes and annual total Hg emissions from individual land use.SHg is the Hg content in surface soil (0-10 cm), FF is the Hg 0 flux over foliage, and FS is the Hg 0 flux over soil.

Diurnal variation of natural Hg 0 emissions in China
Table 3 shows the annual mean air-surface fluxes for different land use types.Annual mean air-foliage fluxes range from −0.2 to −4.5 ng m −2 h −1 , with the highest value over the woody savannas, and the lowest over deciduous forests (Table 3).The diurnal variation for air-foliage flux is displayed in Fig. 4. Higher deposition occurs during midmorning (08:00-10:00 UTC + 8) and late afternoon (16:00-17:00 UTC + 8) due to the suitable air temperature and solar irradiance that induces Hg uptake by stomata.The rates of Hg uptake during midday are comparatively weaker due to the stronger irradiance and higher temperature.This bimodal pattern is consistent with field observations (Lindberg et al., 2002;Poissant et al., 2008;Fritsche et al., 2008;Sommar et al., 2016a), suggesting that the model is capable of simulating the diurnal pattern of air-foliage exchange of Hg 0 .Such a diurnal pattern is caused by re-emission of the deposited Hg on the surface foliage through photoreduction under the strong solar radiation during noontime, and the emissions from underlying soil surfaces.Except for urban lands, the strength of diurnal deposition for the other land uses is controlled by LAI, solar radiation, and air temperature.The elevated atmospheric Hg concentration is an important parameter to induce Hg uptake by growing foliage in urban lands.Simulated mean air-soil fluxes range from 0.1 to 23.3 ng m −2 h −1 , with the lowest flux over barren vegetated lands and the highest over urban lands (Table 3).This suggests that the simulated air-soil fluxes greatly vary over different land uses.There are distinct diurnal variations in terrestrial ecosystems (Fig. 5).The diurnal pattern is caused by the variation of solar irradiance, close to zero at night and peaking from 13:00 to 15:00 (UTC + 8).Similar diurnal patterns have been observed during filed measurements for forest, grassland, and cropland in China (Feng et al., 2005;Fu et al., 2008Fu et al., , 2012;;Zhu et al., 2013).The degree of diurnal variability for each land use is highly related to the LAI.Higher LAI gives a more intensive canopy shading and largely inhibits Hg evasion from soil under canopy.This is also the main reason for relative weaker diurnal variation over forest soils compared to shrubland, grassland, and cropland (Fig. 5).The synergistic interactions between low vegetation cover and high soil concentration (Mean = 162 ± 83 ng g −1 ) results in the strongest diurnal variation for urban land types.
The simulated annual mean of air-water flux is 3.4 ng m −2 h −1 .The diurnal variability for air-water flux is weaker since wind speed is a more influential driver than irradiance (X.Wang et al., 2014), consistent with the diurnal variation observed in field studies that meteorology and photochemical process are the primary factors (Feng et al., 2002(Feng et al., , 2003(Feng et al., , 2008;;Wang et al., 2006;Fu et al., 2010Fu et al., , 2013a, b), b).
Overall, the annual net natural emission in China is 465.1 Mg Hg (Table 3), including 565.5 Mg yr −1 from soil, 9.0 Mg yr −1 from water bodies, and −100.4Mg yr −1 from the vegetation.The annual quantity of emission from soil is comparable to the estimate (528 Mg yr −1 ) based on the scale-up calculation using measured air-soil fluxes (Fu et al., 2015) that suggest emissions from cropland and grassland are the most important contributor.Of the total Hg 0 emission estimated by the model, 50 % is from shrubland, savanna, and grassland (C6-C11, 38 % total land use); 33 % is from cropland (C12-C13, 22 % total land use); 9 % is from forest (C1-C5, 14 % total land use); and 8 % is from other land use types.Forest contributes to 28 % of Hg uptake by foliage; shrubland, savanna, and grassland contribute to 38 %; cropland contributes to 33 %; and other land use types contribute to 1 %.
Although soil Hg contents in forest ecosystems are 2-4 times higher than that in grassland and cropland, total annual fluxes above the canopy (soil plus foliage) of forest ecosystems are 1-6 times lower than the values in other two types of land uses (Table 3).This highlights the importance of canopy cover in natural emission process of Hg 0 .It is noteworthy that the land use data are based on the survey in 2000 (Ran et al., 2012).From 2000 to 2013, the forested area in China increased from 14.0 to 21.6 % (FAO, 2014), benefiting from implementation of the governmental Grain for Green Project and stricter natural forest protection actions.Assuming that annual mean air-surface fluxes are at the same level as in this study, the total quantity of natural Hg emission in 2013 is approximately 5 % smaller than this estimate because of the increasing forest coverage.Given the forest coverage is projected to be 24-26 % during 2030-2050(FAO, 2014)), the quantity of natural Hg emission in China would decrease by 9-10 % compared to the estimated level of 2013.

Spatial distribution of natural Hg emission in China
The spatial distribution of annual air-foliage flux can be divided by the well-known geodemographic demarcation line, the Heihe-Tengchong Line (Fig. 6a).The vegetation on the east side of the line is much denser than on the west side of the line because of the higher annual precipitation (≥ 800 mm; Fig. S2) that leads to stronger Hg 0 uptake by vegetation (> 90 % of the grid cells have a flux below −1.0 ng m −2 h −1 on the east side, compared to > 90 % of the grid cells have a flux above −0.5 ng m −2 h −1 on the west side).There is an enhanced Hg deposition in South China (22-27 • N, 105-113 • E; Fig. 6a), with the fluxes ranging −3.8 to −19.1 ng m −2 h −1 .This can be explained by the elevated atmospheric Hg concentrations (2-10 ng m −3 ) and dense vegetation (i.e., high LAI; Fig. S3) that enhance Hg uptake (Fu et al., 2015;Zhu, 2014).Specifically, evergreen broadleaf forest has the highest LAI compared to other type of forests (Liu et al., 2012) and shows enhanced Hg uptake (up to −4.5 ng m −2 h −1 mean flux).Although the direct measurement of Hg deposition flux through vegetative uptake is still not presently feasible, the measured Hg input through litterfall (Fu et al., 2015) suggested the rate of Hg uptake by foliage could be in the range of 4-12 ng m −2 h −1 , comparable to the simulation results in this study.
Figure 6b shows the spatial distribution of annual airsoil fluxes.There are three high flux regions (mean flux ≥ 10 ng m −2 h −1 ): cropland/grassland in South China and Southwest China (mainly in the Guangdong, Guangxi, Guizhou, Yunnan, Chongqing, and Sichuan provinces), cropland in North China (the Heibei, Henan, and Shangdong provinces), and grassland in North China (Inner Mongolia, Shaanxi, and Xinjiang provinces).Such elevated fluxes in the first two regions have been confirmed in field studies www.atmos-chem-phys.net/16/11125/2016/(Feng et al., 2005;Wang et al., 2005;Wang et al., 2006;Fu et al., 2008;Sommar et al., 2016b).The high fluxes in South China and Southwest China are attributed to the elevated Hg concentration in soil (85 % of grid cells have a soil Hg content ≥ 100 ng g −1 ; Fig. 1).Interestingly, soil Hg content is not the primary factor causing the high fluxes in the other two regions (70 % of grid cells have a soil Hg content ≤ 50 ng g −1 ).Dry deposition of PBM and/or GOM plausibly supply the reducible Hg in soil for gradual reduction and volatilization as Hg 0 (Sommar et al., 2016b).The relatively low LAI (Fig. S3), strong solar irradiance, and high soil temperature (Figs.S4-S5) during summer/autumn contribute to the high simulated emissions.The lower simulated fluxes in desert regions compared to fluxes over grassland (Fig. 6b) are caused by the lower fraction of reducible Hg in soils.
Since the soil Hg 0 flux is the primary source of natural Hg emission, the spatial distribution of the natural Hg emission is strongly influenced by air-soil flux (Fig. 7a).There is a distinct seasonal variation in the emission quantity: 8 % in winter, 28 % in spring, 51 % in summer, and 13% in autumn (Fig. 7b-e).Elevated fluxes mainly cluster in South China and Southwest China during winter because of higher soil Hg content (Fig. 1), and relatively higher temperature and stronger irradiance.Highest correlation coefficients are found between the flux and soil Hg concentration and soil bulk density (Table 4), suggesting that the soil Hg 0 pool is a major factor influencing Hg emission in winter.From the cold to warm season, fluxes gradually increase from low latitude to high latitude with the seasonal change of temperature and solar radiation (Fig. 7b-d).Under the strong irradiance and temperatures in summer, > 65 % of the grid cells in the domain have a flux above 10 ng m −2 h −1 and the effect of soil Hg content is relatively weaker (Table 4).In autumn, higher flux occurs over the cropland of Central and North China, and over the regions with high soil Hg content (Fig. 7f) because of the decreasing temperature and solar irradiance in the south and the influence of soil Hg content in the north.Overall, 72 % of natural Hg emission occurs from May to September, with higher emission over grassland and cropland in North China in these months.
It is worth noting that parts of regions in South China (23-31 • N, 110-120 • E; mainly in the Fujian, Jiangxi, Hunan, Hubei, and Anhui provinces) and Northeast China (39-51 • N, 130-134 • E; mainly in the Liaoning, Jilin, and Heilongjiang provinces) have relatively lower fluxes (−6.9-9.0 ng m −2 h −1 ) during summer and autumn (Fig. 7e-d).In  addition to the impact from the intensive canopy cover in forested areas (Fig. S2), agricultural activities in these regions also contribute to the smaller fluxes.Based on Liu et al. (2013), 60 % croplands in these regions are flooded for rice planting in summer and autumn.Field-scale flux measurement using micrometeorological methods (i.e., aerodynamic gradient method) suggest that a typical oilseedrice rotated cropland in Southwest China is a significant source during oilseed planting seasons with fluxes of 10.1-89.4ng m −2 h −1 and a mild sink during rice planting seasons with fluxes of −3.4 to −15.8 ng m −2 h −1 (Zhu, 2014).The model also successfully simulates such a pattern, with simulated fluxes at 1.1-101.5 ng m −2 h −1 (Fig. 7b-c) during winter and earlier spring when croplands are not flooded and at −3.5 to 1.5 ng m −2 h −1 during the rice growing season (Fig. S6).Overall, 3.3 Mg Hg 0 is predicted to deposit into rice paddies during the rice growing season, with 56 % of the deposition occurring in summer, 41 % in autumn, and 3 % in late spring.

Verification of model estimates
For the first time, the simulated natural Hg emission in China is verified against field observational data in this study.The dataset of Hg deposition through litterfall is utilized for verifying the simulated air-foliage fluxes because of two rea-sons: (1) it has been shown that Hg deposition through litterfall dominates dry deposition (≥ 70 %) in the forests of China (Fu et al., 2015), and the annual Hg deposition through litterfall has been used as a surrogate to constrain air-foliage fluxes in forest ecosystems (Risch et al., 2012;Zhang et al., 2012a); and (2) the litterfall data in China comprehensively include different forest types compared to the limited locations where air-foliage flux data are available.For verifying the exchange fluxes over water and soil surfaces, the flux measurements over forest soil, grassland, cropland, and water bodies in China (Table S2) are utilized.
To estimate the annual Hg deposition through litterfall in the study domain, a Monte Carlo simulation (described in details in the Supplement) is applied for constructing the probability distribution of the Hg deposition through litterfall using litter biomass production and litterfall Hg concentration in China reported in peer-reviewed literature (Fig. 8).The sampling locations include 20 sites in Tibetan Plateau, 27 sites for evergreen forests, and 12 sites for deciduous forests.The quality-assured data of litter biomass production (number of replicates ≥ 3, collector size = 1 m 2 ) are obtained from the China National Knowledge Infrastructure (CNKI).This dataset includes the measurements at 5 sites in Tibetan Plateau, 277 sites for evergreen forests, 74 sites for deciduous forests, and 61 sites for mixed forests.
The measured air-soil flux (Table S2) ranges from −1.4 to 20.7 ng m −2 h −1 over forest soil (n = 19; mean = 6.1 ± 5.1 ng m −2 h −1 ), from −18.7 to 114 ng m −2 h −1 over grassland (n = 14; mean = 26 ± 36 ng m −2 h −1 ), from −4.1 to 135 ng m −2 h −1 over cropland (n = 33; mean = 21.3 ± 36.7 ng m −2 h −1 ).For water body (n = 51), the flux range is 0-43.8ng m −2 h −1 with a mean of 4.6 ± 6.6 ng m −2 h −1 (Table S2).The mean fluxes from May to October are substantially higher than those from November to April: 3.3 times for water surface (p = 0.004), 3.2 times higher for forest soil (p = 0.08), and 1.4 times for cropland (p = 0.50).A reverse  S2), (c) comparison between simulated exchange and measured exchange over soil under canopy, and (d) comparison between simulated exchange and measured exchange over grasslands, cropland, and water surface.The mean and median of Fig. 8c and b are based on the filed data from peer-reviewed literature (n = 19 for forests; n = 12 for grasslands; n = 42 for croplands; n = 51 for water bodies).Note that the exchange over deciduous needleleaf forests in (a) is small because of the small forest area.
trend is found for grassland, which has higher mean flux in cold seasons (50 % higher, p = 0.36).More measurements in grassland where few data exist will greatly improve the accuracy of the current estimate.
Figures S8, S9, and 8a compare the model estimates to the mean and variability level predicted by a Monte Carlo simulation using field data.The annual Hg uptake simulated by the bidirectional exchange model is not significantly different from the field observations (p > 0.05, t test) and the spatial patters are similar (Fig. S8) in coniferous forest ecosystems, demonstrating the model capability for simulating the air-foliage flux.However, the bidirectional exchange model did not capture the spatial distribution of air-foliage flux in broadleaf forest ecosystems (particularly in evergreen broadleaf forest; Fig. S9).One possible explanation is that the resistance terms obtained from temperate/boreal forests (Zhang et al., 2012b) may not appropriately represent the value in evergreen broadleaf forests.Filed measurements suggests that the leaf stomatal conductance of broadleaf is usually higher than the value of needleleaf (Wang et al., 2015;Ishida et al., 2006;Sobrado, 1991;Eamus et al., 1999), leading to a larger air-foliage Hg 0 exchange (Graydon et al., 2006).Further studies on the Hg transport and chemical reactions at the air-foliage interface in evergreen broadleaf forests will help constrain the model.
Figure 8b shows the scatter plot of the measured and model-predicted fluxed over soil and water (R 2 = 0.73).Modeling results for over water surfaces and soil under forest canopy also agree with filed data (Fig. 8c-d).
The model results somewhat underestimate the high fluxes (≥ 30 ng m −2 h −1 ; Fig. S9.2) measured over grassland and cropland (Fig. 8d), which can be attributed to several possible reasons.One reason is the bias caused by the comparatively coarser spatial resolution (36 km) of meteorological parameters and soil properties that limit the reproduction of the instantaneously measured fluxes at observational sites.In addition, the limited mechanistic understanding of the reemission process after Hg dry deposition (Gustin et al., 2015(Gustin et al., , 2008b;;Lindberg et al., 2007;Ariya et al., 2015) complicates model parameterization.Finally, the uncertainties caused by flux quantification methodology (Lin et al., 2012;Zhu et al., 2015a, b) and the typically short campaign periods (mostly ranging from several days to a couple of weeks) could bias the measurement data (Feng et al., 2005;Fu et al., 2008Fu et al., , 2012Fu et al., , 2015;;Zhu et al., 2015a, b).Improvement of flux methods and extended campaign periods at more study sites for cropland/grassland will help constrain the model estimates.

Comparison with earlier estimates and implications on Hg emission outflow in China
Figures S10 and S11 show the gridded natural Hg emission in the east Asian domain reported by Shetty et al. (2008) and S. Wang et al. (2014), which have two distinct differences compared to the model estimate in this study.One is the role of vegetation in natural Hg emission, the other is spatial distribution of the emission.Vegetation is clearly an important sink of Hg 0 with the mechanistic model algorithms implemented in this study and the shading of vegetation suppresses Hg evasion from soil under canopy.In contrast, vegetation is considered a major source, accounting for 76 % of total emissions in Shetty et al. (2008) because the model treats Hg evasion as an evapotranspiration process that transports Hg from root zone through vascular tissues in foliage (Gbor et al., 2007;Shetty et al., 2008).Recent experimental evidence using stable Hg isotope tracers points to exclusion of this pathway for cereal plants (Cui et al., 2014).In addition, Hg isotopic signatures observed in air and leaf samples (Demers et al., 2013;Yin et al., 2013) and during air-foliage exchange process (Graydon et al., 2006;Gustin et al., 2008a) indicate uptake of atmospheric Hg by foliage, pointing to vegetation as an Hg 0 sink.In contrast to the spatial distribution of the emission in this study, the earlier Hg 0 emission estimates occur mainly in the regions on the east side of the Heihe-Tengchong Line (Shetty et al., 2008;S. Wang et al., 2014), due to the spatial distribution of vegetation (Shetty et al., 2008), or soil Hg content (S.Wang et al., 2014).This study advances upon the earlier estimates (Shetty et al., 2008;S. Wang et al., 2014) in three areas.Firstly, the recent soil survey data, including soil Hg content and other soil characteristics, are a major advantage in this study.The soil Hg data applied in Shetty et al. (2008) are outdated with a coarse spatial resolution, while the data in S. Wang et al. (2014) are based on the output of the global terrestrial Hg model in GEOS-Chem, calculated from Hg / C ratios.Secondly, the mechanistic model scheme better describes the air-surface exchange process compared to the regression and evapotranspiration in the earlier studies.Finally, the model estimates are verified against the field flux data with generally good agreement, which has not been attempted in earlier works.
Although the total quantity of annual natural emission estimated in this study is comparable to earlier estimates (400-600 Mg yr −1 ) by Shetty et al. (2008) and S. Wang et al. (2014), the distinct spatial distribution of natural emissions simulated in this study may alter the current understanding of Hg emission outflow from China assessed by Lin et al. (2010b).The outflow of Hg emissions in China is mainly driven by the prevailing west-wind drift (Lin et al., 2010b;Chen et al., 2014).The predominant natural Hg emission in the west side of model domain results in a longer residence time of evaded Hg, which can be more readily oxidized and deposited within the domain.Furthermore, the dense vegetation in the east side of the domain can also help capture the nature Hg 0 emission, potentially leading to substantially larger domestic deposition and smaller quantity of outflow compared to the estimates by Lin et al. (2010b).We are presently reassessing the emission outflow using a regional chemical transport model (e.g., CMAQ-Hg) and a similar mass balance approach by Lin et al. (2010b) and will report the model results in a future paper.

Conclusions
Using a mechanistic model incorporating the present state of understanding in Hg transformation in soils and on foliage surface with up-to-date datasets of soil characteristics and land use changes, the natural emission of Hg 0 vapor in China is estimated to be 465.1 Mg yr −1 , including 565.5 Mg yr −1 of emission from soils, 9.0 Mg yr −1 of emission from water bodies, and −100.4Mg yr −1 deposition (uptake) by vegetation.The air-surface exchange is strongly dependent on land use and meteorology, with 9 % of net emission from forest ecosystems; 50 % from shrubland, savanna, and grassland; 33 % from cropland; and 8 % from other land uses.Given the large agricultural land area in China, farming activities play an important role on the air-surface exchange.Particularly, rice fields shift from a net sink (3.3 Mg uptake) during the growing season in rice paddy to a net source during the season when the farmland is not flooded.The estimated natural Hg 0 emission in this study yields similar Hg 0 evasion quantity but exhibits contrasting spatial distribution compared to the estimate by Shetty et al. (2008).The difference in the spatial patterns may alter the current understanding of Hg emission outflow from China as reported by Lin et al. (2010b) because of a substantial amount of natural Hg 0 emission occurs in West China, which requires further assessment.
For future model improvement, studies focusing on fundamental understanding of Hg II reduction in soil (especially the role of soil organic matter, contribution of photochemical and nonchemical pathways, and radiation transfer in soil) and air-vegetation exchange mechanisms are needed.Continuous updates on the data of soil characteristics and Hg content are also essential.More field data for model performance evaluation are also important for constraining the model.In particular, data of air-foliage flux, and air-soil flux over cropland and grassland in the remote regions of North China are valuable for model calibration.
The Supplement related to this article is available online at doi:10.5194/acp-16-11125-2016-supplement.

Figure 1 .
Figure 1.Updated Hg concentrations (ng g −1 ) in surface soil (0-20 cm) of China.Sampling areas in NMPRGS cover most agriculturally and industrially developed regions of East and Central China, and are presented in more detail in Li et al. (2014).
. The model algorithms are coded in FORTRAN 90 and Network Common Data Form (NetCDF) version 4.3.The gridded model results are visualized by the Visualization Environmental for Rich Data Interpretation (VERDI) version 1.5.

Figure 3 .
Figure 3. Sensitivity analysis on model parameters for air-soil exchange using a two-level factorial design after prescreening the model variables shown in Table 2 for the identified significant factors.The effects shown in the figure are based on a significance level of 95 % (i.e., p < 0.05).

Figure 4 .
Figure 4. Diurnal variation of mean simulated exchange fluxes of Hg 0 over canopy in the model domain (UTC + 8).

Figure 5 .
Figure 5.Diurnal variation of mean simulated exchange fluxes of Hg 0 over soil and water surfaces in the model domain (UTC + 8).

Figure 6 .
Figure 6.Simulated results of (a) mean annual air-foliage flux and (b) mean annual air-soil flux in the study domain.

Figure 7 .
Figure 7. Model estimates of (a) annual mean Hg 0 fluxes in the model domain, (b) seasonal mean Hg 0 fluxes in winter, (c) seasonal mean Hg 0 fluxes in spring, (d) seasonal mean Hg 0 fluxes in summer, (e) seasonal mean Hg 0 fluxes in autumn, and (f) monthly Hg 0 fluxes in the grid cells (box and whisker chart showing maximum, 75th percentile, mean, median, 25th percentile, and minimum).

Figure 8 .
Figure 8. Model verification: (a) model estimates of Hg 0 uptake by foliage (which include the uptake by stoma and the re-emission from cuticle) and by stoma, compared to the estimate (mean and 95 % confidence interval) of Hg 0 uptake using a Monte Carlo (M-C) simulation of the observational data; (b) scatter plot of the observed fluxes vs. simulated fluxes for different land uses (the flux observations are described in detailed in TableS2), (c) comparison between simulated exchange and measured exchange over soil under canopy, and (d) comparison between simulated exchange and measured exchange over grasslands, cropland, and water surface.The mean and median of Fig.8c and bare based on the filed data from peer-reviewed literature (n = 19 for forests; n = 12 for grasslands; n = 42 for croplands; n = 51 for water bodies).Note that the exchange over deciduous needleleaf forests in (a) is small because of the small forest area.

Table 1 .
Model variables, constants, and rate coefficients used in the model simulation.

Table 2 .
Examined model variables and the experimental levels of factorial design for air-soil exchange.

Table 4 .
Pearson correlations between mean total fluxes and major controlling environmental parameters in each season.
a means p < 0.01 and b means p < 0.05.