Mercury vapor air – surface exchange measured by collocated micrometeorological and enclosure methods – Part I : Data comparability and method characteristics

Reliable quantification of air–biosphere exchange flux of elemental mercury vapor (Hg) is crucial for understanding the global biogeochemical cycle of mercury. However, there has not been a standard analytical protocol for flux quantification, and little attention has been devoted to characterize the temporal variability and comparability of fluxes measured by different methods. In this study, we deployed a collocated set of micrometeorological (MM) and dynamic flux chamber (DFC) measurement systems to quantify Hg flux over bare soil and low standing crop in an agricultural field. The techniques include relaxed eddy accumulation (REA), modified Bowen ratio (MBR), aerodynamic gradient (AGM) as well as dynamic flux chambers of traditional (TDFC) and novel (NDFC) designs. The five systems and their measured fluxes were cross-examined with respect to magnitude, temporal trend and correlation with environmental variables. Fluxes measured by the MM and DFC methods showed distinct temporal trends. The former exhibited a highly dynamic temporal variability while the latter had much more gradual temporal features. The diurnal characteristics reflected the difference in the fundamental processes driving the measurements. The correlations between NDFC and TDFC fluxes and between MBR and AGM fluxes were significant (R>0.8, p<0.05), but the correlation between DFC and MM fluxes were from weak to moderate (R = 0.1–0.5). Statistical analysis indicated that the median of turbulent fluxes estimated by the three independent MM techniques were not significantly different. Cumulative flux measured by TDFC is considerably lower (42 % of AGM and 31 % of MBR fluxes) while those measured by NDFC, AGM and MBR were similar (<10 % difference). This suggests that incorporating an atmospheric turbulence property such as friction velocity for correcting the DFC-measured flux effectively bridged the gap between the Hg fluxes measured by enclosure and MM techniques. Cumulated flux measured by REA was ∼ 60 % higher than the gradient-based fluxes. Environmental factors have different degrees of impacts on the fluxes observed by different techniques, possibly caused by the underlying assumptions specific to each individual method. Recommendations regarding the application of flux quantification methods were made based on the data obtained in this study.

fluxes estimated by the three independent MM techniques were not significantly different.Cumulative flux measured by TDFC is considerably lower (42 % of AGM and 31 % of MBR fluxes) while those measured by NDFC, AGM and MBR were similar (<10 % difference).This suggests that incorporating an atmospheric turbulence property such as friction velocity for correcting the DFC-measured flux effectively bridged the gap between the Hg 0 fluxes measured by enclosure and MM techniques.Cumulated flux measured by REA was ∼ 60 % higher than the gradient-based fluxes.Environmental factors have different degrees of impacts on the fluxes observed by different techniques, possibly caused by the underlying assumptions specific to each individual method.Recommendations regarding the application of flux quantification methods were made based on the data obtained in this study.

Introduction
Mercury (Hg) is a ubiquitously distributed neurotoxin in the environment (Lindqvist et al., 1991).The bulk of atmospheric Hg is made up of gaseous elemental Hg (Hg 0 , > 95 % of the total mass) with minor contribution from the analytically defined fractions of gaseous oxidized Hg (GOM) and particulate bounded Hg (PBM) (Gustin, 2011).Being chemically inactive and partitioning less favorably into aqueous Published by Copernicus Publications on behalf of the European Geosciences Union.

W. Zhu et al.: Mercury flux data comparability and method characteristics
phase, Hg 0 is prone to undergo hemispherical-scale tropospheric transport (Durnford et al., 2010).Hg 0 is subject to bi-directional exchange between atmosphere and natural surfaces through complex and yet not well understood processes (Bash, 2010;Gustin and Jaffe, 2010).Recent estimation indicates that annual natural emission accounts for two-thirds of global release of atmospheric Hg (Pirrone et al., 2010).However, current estimates of natural exchange quantity remain highly uncertain due to the limitations in accuracy and representativeness of measurement techniques (Gustin and Jaffe, 2010;Pirrone et al., 2010).
There exist multiple experimental approaches to gauge Hg 0 air-surface exchange, which can be grouped into enclosure and micrometeorological (MM) methods (Sommar et al., 2013a).Dynamic flux chambers (DFCs) representing the smallest scale, as the areas covered are typically in the order of 0.1 m 2 , are the most extensively applied method for quantifying Hg 0 evasion from and deposition to soil (Poissant and Casimir, 1998;Stamenkovic and Gustin, 2007;Xiao et al., 1991;Carpi and Lindberg, 1998).For measuring Hg 0 fluxes on larger landscape scales, MM techniques represent an attractive alternative to DFCs.They allow spatially averaged measurements over a large area without disturbing ambient environmental conditions.For trace gases such as CO 2 , CH 4 , O 3 , NH 3 , HNO 3 , and selected volatile organic compounds (VOCs), eddy covariance (EC) is the preferred MM technique for quantifying air-landscape gas exchange (Aubinet et al., 2012;Farmer et al., 2006;Park et al., 2013;Whitehead et al., 2008).However, due to the lack of a sufficiently fast and sensitive sensor for the ultra-trace levels of Hg 0 in air, true EC measurement of background Hg 0 flux has not yet been accomplished.MM techniques applied in Hg 0 flux (also called turbulent flux) quantification include the relaxed eddy accumulation method (REA, also known as conditional sampling, CS) (Bash and Miller, 2008;Cobos and Baker, 2002;Olofsson et al., 2005;Sommar et al., 2013b), the aerodynamic gradient methods (AGMs) (Baya and Van Heyst, 2010;Cobbett and Van Heyst, 2007;Converse et al., 2010;Edwards et al., 2005;Fritsche et al., 2008a;Fritsche et al., 2008b;Marsik et al., 2005), and the modified Bowen ratio method (MBR) (Converse et al., 2010;Fritsche et al., 2008a;Fritsche et al., 2008b;Lindberg et al., 1995;Poissant et al., 2004).MM methods estimate turbulent transport with the assumptions of fetch homogeneity and the measurements are made within the constant flux layer (Wesely and Hicks, 2000).For example, REA-derived flux relies on accurate measurement of the concentration difference between upward and downward moving air parcels while gradientderived flux is estimated from the vertical concentration gradient and the associated turbulent exchange parameters.For the traditional DFC (TDFC) methods, flux is derived from a steady-state mass balance over the chamber.More recently, we have designed and deployed a DFC of novel design (NDFC) based on surface wind shear condition (friction ve-locity) rather than on artificial fixed flow to account for natural shear conditions (Lin et al., 2012).
Limited efforts have been devoted to Hg 0 flux measurement comparison.In the Nevada STORMS campaign (4-day duration), TDFCs and MM gradient methods were deployed to measure Hg 0 flux over a heterogeneously Hg-enriched fetch.The TDFC-and MM-derived fluxes differed by one order of magnitude (Gustin et al., 1999;Gustin and Lindberg, 2000;Poissant et al., 1999;Wallschläger et al., 1999).Subsequent investigations have suggested that TDFCs of different sizes, shapes and operation flow rates yield different fluxes (Eckley et al., 2010;Lin et al., 2012;Zhang et al., 2002;Wallschläger et al., 1999).Gradient methods were deployed to measure seasonal Hg 0 fluxes over grasslands in the Alps (Fritsche et al., 2008b) and over a meadow in the Appalachians (Converse et al., 2010), the observed flux means varied by up to one order of magnitude.Collocated flux measurement using both MM and DFCs techniques for method evaluation and data synthesis remains scarce (Gustin, 2011).This limits a thorough comparison of flux data obtained by different techniques.
Measured fluxes are estimates of unknown quantities of air-surface exchange under field conditions and a reference technique for validating the estimates does not exist.Each available technique has its specific advantages and drawbacks and its applicability to obtain representative fluxes is limited under particular atmospheric conditions and site characteristics.It is therefore essential to compare and review uncertainties of the major techniques deployed for measuring air-ecosystem Hg 0 exchange.The objective of this study is to investigate the method characteristics, data comparability and measurement uncertainty of Hg 0 exchange fluxes as measured by five collocated MM and DFC methods including REA, MBR, AGM, TDFC and NDFC.We improved a number of measurement platforms (Lin et al., 2012;Sommar et al., 2013b) and performed two intensive field campaigns over both bare and vegetated landscapes.The results of this integrated assessment are presented in part by two companion papers.In Part I, we evaluate the technical merits of the examined flux quantification methods, assess the flux variability and data comparability, and address the method applicability under a given set of environmental conditions.In Part II, we quantify the bias and uncertainty of the examined flux measurement methods.

Dynamic flux chamber techniques
In this study, chambers of traditional and the new design described in Lin et al. (2012) (Feng et al., 2005;Fu et al., 2008Fu et al., , 2010Fu et al., , 2012;;Li et al., 2010;Wang et al., 2005Wang et al., , 2007;;Zhu et al., 2013a).The NDFC was fabricated of thin polycarbonate sections and enclosed a soil surface of 0.09 m 2 (for details, see Lin et al., 2012).The NDFC internal flow condition was precisely controlled to relate to the applied flushing flow rate to the atmospheric boundary shear condition (therefore wind shear condition) and the calculated flux was re-scaled to boundary shear condition (Eq. 2 below).Both DFCs were operated at a relatively high flushing flow rate of 15 L min −1 , corresponding to turnover times (TOTs) of 0.32 min and 0.47 min for TDFC and NDFC, respectively.The flux from TDFC and NDFC were calculated following Eq.( 1) and (2), respectively (Xiao et al., 1991;Lin et al., 2012): , where F TDFC Hg 0 is Hg 0 flux measured from the TDFC method, F NDFC Hg 0 is Hg 0 flux from the NDFC method, Q is applied flow rate (0.9 m 3 h −1 ), A is footprint (0.06 m 2 for TDFC, 0.09 m 2 for NDFC), C o and C i are the DFC outlet and inlet air Hg 0 concentration, and k mass(a) and k mass(m) are the overall mass transfer coefficient (m s −1 ) in the near-surface boundary layer and in the internal layer within NDFC, respectively.A c is the NDFC flow cross-sectional area (0.009 m 2 ), l is the distance measured from the starting point of the measurement zone (0.15 m), h is the height of NDFC (0.03 m), u * is the atmospheric boundary layer friction velocity, and z 0 is surface roughness height (m).D H and D are the NDFC hydraulic radius (0.0545 m) and diffusivity of Hg 0 (1.194 × 10 −5 m 2 s −1 ), respectively.

Relaxed eddy accumulation (REA) method
A REA system of whole-air type was deployed with the design and operation parameters described elsewhere (Sommar et al., 2013b;Zhu et al., 2013b).The REA apparatus constitutes of open path EC (OPEC) and conditional gas sampling system.The OPEC part included a 3-D fast-response anemometer, an open path CO 2 /H 2 O analyzer, and a micrologger with processing and control capabilities.MM data collected at 10 Hz are acquired and processed by the latter, which also control the execution of conditional sampling valves from its 12 V terminal following the implemented dynamic wind dead-band algorithm to accurately isolate up-and down-drafts present in sampled turbulent air parcels.Turbulent REA flux was computed according to where σ w (m s −1 ) is the standard deviation of vertical wind speed (m s −1 ) and C ↑/↓ is the concentration of Hg 0 (at standard temperature and pressure) for the up-and down-moving eddies corrected for dilution of zero air injection, respectively (ng m −3 ).The operational form of Eq. ( 3) is given in the right-hand side, in which, for sample i, m ↑/↓ i is the mass of Hg 0 derived for the up-or down-draft channels (pg), t i is the total duration (min), Q ↑/↓ i is the continuous flow rate through the up-or down-draft channels (L dry air min −1 ), and α ↑/↓ i is the fraction of time the up-or down-draft conditional sample valves are activated.β s is a dimensionless relaxation coefficient (calculated from scalar s) which for each averaging period (20 min) was calculated on-line from suitable scalar s those fluxes (F EC s = ρ d • w χ s ) can be measured by the OPEC system (in addition to CO 2 flux, buoyancy flux C P • w T s and for latent heat flux λ • w q , symbol definitions see appendix in Sommar et al., 2013b) as well as by REA according to where χ ↑/↓ s is the mixing ratio of the specific scalar quantity for the up-and downdraft (kg kg −1 ).

Aerodynamic gradient micrometeorological (AGM) method
The AGM method is based on an analogy application of Fick's first law stating that turbulent bi-directional flux of a scalar from surface (F AGM s ) is proportional to its local vertical concentration gradient (∂C/∂z) and eddy diffusivity of sensible heat (K H ), which is a function of friction velocity (u * ) and the dimensionless stability parameter ς m = (z m − d) /L (z m is the sampling height above ground, d is the zero plane displacement height and L is the Monin-Obukhov length (Monin and Obukhov, 1954).Assuming that measurements are made within a vertical layer of constant flux that forms over homogeneous terrain, after integration between W. Zhu et al.: Mercury flux data comparability and method characteristics two heights, the flux can be expressed as where κ is von Kármán constant (∼ 0.41), u * is the friction velocity (m s −1 ), the υ tr term is the transfer velocity (m s −1 ), z 2 and z 1 are the heights of the upper and lower sampling inlet (m), and ψ H is the integrated universal function for sensible heat to correct for deviations from the ideal logarithmic profile.ψ H is parameterized as a function of ς m (ς 1 and ς 2 represent the parameter at z 2 and z 1 respectively), and furthermore C Z 2 and C Z 1 are the Hg 0 concentration (ng m −3 ) at z 2 and z 1 , respectively.

Modified Bowen ratio (MBR) method
The MBR method assumes that the flux of a trace gas can be related to that of a surrogate scalar determined from OPEC measurements (e.g., sensible and latent heat, CO 2 flux, and H 2 O flux) (Converse et al., 2010;Lindberg et al., 1995).In this study, temperature was used as the proxy scalar, which was monitored at the heights coinciding with measurement of Hg 0 concentration.The Hg 0 flux is calculated following Walker et al. ( 2006): where F MBR Hg 0 is the Hg 0 flux (ng m −2 h −1 ) measured with the MBR method, w T is kinematic heat flux (K m s −1 ) measured by EC, while C and T are the vertical gradients of Hg 0 concentration (ng m −3 ) and air temperature (K), respectively.The ratio w T / T is known as the eddy diffusivity for heat.

Site description and sampling
The flux measurement experiments were conducted at Yucheng Comprehensive Experimental Station, Chinese Academy of Sciences (36 • 57 N, 116 • 36 E), which is a semi-rural agricultural station located in the North China Plain approximately 50 km from Jinan, Shandong Province.Within a radius of ∼ 5 km the planting system is winter wheat (Triticum aestivum Linn., November-May) or summer maize (Zea mays, June-October) for a rotation in a year.The surface soil texture in this area is silty loam consisting of 12 % sand, 22 % clay and 66 % silt with moderate salinity and alkalinity (pH = 8.6) (Hou et al., 2012).The agricultural fields adjacent to the sampling site are relatively flat (level differences < 1.5 m within 1 km) and the total Hg content in surface soil is spatial homogeneously distributed (45 ± 3.9 µg kg −1 , n = 27) (Zhu et al., 2015a).Two intensive field campaigns were performed: one in late autumn 2012 (IC#1, 4-24 November, DOY (day of year) 309-329) and the other in spring 2013 (IC#2,(16)(17)(18)(19)(20)(21)(22)(23)(24)(25).IC#1 was carried out over the ploughed bare soil surface using AGM, MBR, TDFC, and NDFC.IC#2 was carried out over wheat canopy (average height ∼ 0.36 m, leaf area index of 3.4) using REA, AGM and MBR.Given the tight row spacing of the grain field, the deployment of DFCs was not permissible during IC#2.

Instrumentation
A 6.5 m MM flux tower was installed at the same location for both campaigns (Fig. 1).The instrumentation system consists of the tower-based MM systems and ground-based DFCs.
The OPEC system consisted of a Campbell CSAT-3 sonic anemometer-thermometer, Licor LI-COR 7500A open-path CO 2 /H 2 O analyzer and HMP155A humidity-temperature sensors, a standard instrumentation combination used in long-term ecosystem instrumentation networks (Mauder et al., 2013).REA sampling inlet was positioned at 2.96 m above ground.By using a set of 2/3-way automated magnetic switching unit (Tekran ® 1110) coupled with an automated Tekran ® 2537B Hg vapor analyzer operated at a flow rate of 0.75 L min −1 , up-and down-draft conditional samples were sequentially routed into the analyzer at 10 min intervals (two 5 min samples).For gradient measurements, the temperature and relative humidity sensor (HMP155A, Vaisala Oy, Finland) housed in radiation shields and corresponding Hg 0 intake was assembled at two heights of 2.96 m and 0.76 m.The two-level Hg 0 vertical gradient profiling system consisted of two separate inlet lines (PFA Teflon), each with an inlet filter (0.2 µm PFA Teflon), were routed to another sampling manifold (Model 1110).Another Hg 0 gas analyzer (Model 2537B) is connected to the outlet of the manifold and the profile inlets are opened one at a time synchronized with 2537B's sampling cycles.The manifold was configured to allow the inlet not in use to be continually flushed by a bypass pump.Both the pump and 2537B are operated at a flow rate of 1.0 L min −1 .An estimate of the vertical Hg 0 concentration gradient was derived every 20 min from measurements of the two heights sequentially, 5 min integrated samples.
The TDFC and NDFC were operated in tandem using one 2537B analyzer (sampling flow rate 1.0 L min −1 ).A fourport automated magnetic dual switching unit (Tekran ® 1115) was utilized to sequentially sample the two DFCs inlet and outlet twice at 2.5 min intervals in the sequential order: inlet of TDFC, outlet of TDFC, inlet of NDFC and outlet of NDFC, thereby retrieving 2.5 L samples for Hg 0 analysis.20 min Hg 0 flux was calculated using Eqs.( 1) and ( 2) for TDFC and NDFC.Prior to sampling, the internal clocks of all instrumentation were synchronized (UTC + 8 h) and therefore the reported fluxes resembled identical 20 min integration periods.

Quality assurance/control (QA/QC), data evaluation and EC flux corrections
The three Tekran ® 2537B analyzers (Fig. 1) were operated and maintained following the standard operation procedures of NADP (2011).The analyzers were regularly calibrated in the laboratory by manual injections of known amount of Hg 0 .The yielded recovery was 98-101 %.In the field, instruments were calibrated every 48 h using the internal Hg 0 permeation source.A soda-lime trap and a 0.2 µm Teflon membrane filter were located upstream of the inlet of all analyzers.The analyzers are sensitive to insufficient power and were therefore always supplied with grid power passing a 10 kW voltage stabilizer to ensure proper operation in the field.All the tubing and system valve blanks were checked before and after the campaigns by flushing with zero air obtained from a zero-air generator (Tekran ® 1100).Before the field measurement, the accuracy of two HMP 155A sensors was evaluated after periods of side-by-side measurements.
The two DFCs were cleaned by 10 % HNO 3 and Milli-Q water prior to field deployment.Chamber blanks performed at the field site were consistently low for both DFCs (TDFC: and not subtracted upon calculation of fluxes.The REA-system enabled a mode during which air is sampled synchronously with both conditional inlets.This reference mode provides an automated QC-measure to regularly check for gas sampling path bias, while the gradient-based MM techniques require manual testing by collocating gas sampling inlets and sensors.Such side-by-side tests were performed before or after a campaign.Post-processing of collected 10 Hz EC raw data was performed for each of the 20 min flux averaging periods using Eddypro TM 5.0 flux analysis software package (LI-COR Biosciences Inc.).A series of standard data corrections were implemented following Sommar et al. (2013b) including the Webb-Pearman-Leuning (WPL) correction.Moreover, tests were applied on 20 min fast time (10 Hz) series raw data to qualitatively assess turbulence for the assumptions required of applying MM methods (steady-state conditions and the fulfillment of similarity conditions).The basic flag system of Mauder and Foken ( 2004) was utilized to indicate limitation in turbulence mixing, quality indices of 0, 1 and 2 denoting high, moderate and low quality.
For neutral stratification, this ratio is approximately constant at 1.13-1.35(Nemitz et al., 2009).The median σ w /u * was 1.28 and 1.24 during IC#1 and IC#2.However, the variability introduced by diabatic condition is comparatively more pronounced during IC#1.Hg 0 observations at the sampling site showed a wide range of 1.20 to 8.17 ng m −3 (medians 3.12 ng m −3 and 3.50 ng m −3 during IC#1 and IC#2, respectively).The medians were elevated compared to the hemispheric background (1.5-1.7 ng m −3 ), but nevertheless appeared representative of a semi-rural area of North China plain (∼ 3.2 ng m −3 , Zhang et al., 2013).The angular distribution of Hg 0 observations (Fig. 3b, d) indicated a weak Hg 0 concentration dependence on wind direction during IC#1 but a more manifest dependence appeared during IC#2, with elevated concentrations associated with southerly and southwesterly winds (4.04-4.88ng m −3 , 45-130 % higher than those associated with easterlies, 2.12-2.79ng m −3 ).

Characteristics of DFCs Hg 0 fluxes
Descriptive statistics of the DFC Hg 0 flux observations are presented in Table 1.In a comparison, NDFC-derived Hg 0 fluxes spanned over a broader range and exhibited a higher mean.Figure 4a displays the time series of Hg 0 fluxes gauged by the two DFC methods.Both series showed similar diurnal features with daytime evasion (maximum occurred at The median ± MAD (median absolute deviation) of Hg 0 flux was −0.9 ± 3.2 and −1.7 ± 4.3 ng m −2 h −1 for TDFC and NDFC, respectively.Probability plots of both DFC data sets showed positive kurtosis (3.0 and 4.1) and skewness (1.6 and 2.1) (Fig. 5) as a consequence of stronger emission and increased friction velocity at daytime.The substantial fraction of NDFC data points elevated in magnitude outlying the 1.5 • IQR (interquartile range) bound is associated with periods of high wind speed (i.e., showing the dependence of friction velocity in Eq. 2).Moreover, as indicated in Fig. 5, the shortest half (50 %) of the chamber flux data is positioned more towards dry deposition for the novel compared to the traditional chamber technique.Nevertheless, the intrinsic divergence of the microenvironment inside enclosures in relation to that of near-surface air layer tends to promote efflux.

Comparison of Hg 0 fluxes obtained from DFCs measurement
In the Nevada STORMS campaign, seven flow-through enclosures (DFCs) with different operational parameters and designs were located in an arid area with naturally Hgenriched substrate.The observed DFC Hg 0 fluxes showed similar diurnal profiles but diverged in magnitude by an order of magnitude (Gustin et al., 1999;Wallschläger et al., 1999; Gustin and Lindberg, 2000).The observed difference was partially attributed to the substrate heterogeneity with respect to Hg content.In this study, the surface soil Hg content within the methodological footprint range is at large homogeneous and therefore does not pose an interfering factor.Eckley et al. (2010) examined experimentally a series of operational and instrumental factors that may influence DFCderived flux.The DFC flushing flow rate was identified to have substantial positive influence.In the present study, the TOT of TDFC is 50 % smaller than that of the NDFC.Moreover, the footprint of the traditional type is about two-thirds of the NDFC footprint and therefore a higher flux is expected using the NDFC method (Eckley et al., 2010;Lin et al., 2012).Figure 6 shows a scatterplot of the fluxes measured by the NDFC and TDFC approach before and after turbulence correction.The data were significantly positive correlated (R = 0.93, R = 0.95 between TDFC and NDFC fluxes calculated with Eq. ( 2) and Eq. ( 1) respectively; p<0.01).Quantitatively, direct measured flux was consistent for the two chambers (slope 1.01).After accounting for the atmospheric boundary shear condition by Eq. ( 2), the welldeveloped turbulence (higher friction velocity, Fig. 2) during daytime caused the NDFC-inferred Hg 0 flux to be approximately 2.5 times higher than the TDFC flux.Given that fluxes derived from a DFC of conventional type do not allow for rescaling to represent natural surface shear stress conditions, TDFCs are prone to underestimate the soil Hg emission, particularly when operated at low air exchange rates.The ability to incorporate an atmospheric turbulence property such as friction velocity makes the NDFC method a more favorable approach for estimating Hg 0 gas exchange over soils compared to the TDFC method.) and Hg 0 flux (ng m −2 h −1 ) derived from the turbulent diffusion methods (MBR and AGM).Hg 0 concentration gradients were observed in the similar ranges of −0.49 to 0.33 and −0.48 to 0.25 ng m −4 in both campaigns (Table 1 and Fig. 4), though the more occasionally shifting conditions of weak and developed turbulence in IC#1 tend towards promoting a higher scale of diurnal gradient variability (IC#1 vs. IC#2 standard deviation: 0.09 vs. 0.06).Our gradient observations are in alignment with measurement over temperate grasslands (−0.40 to 0.27 ng m −4 ) (Fritsche et al., 2008b).Basic statistics of the MM Hg 0 flux observations is presented in Table 1.The variability in our observations is similar with those reported from previous studies using the MM flux measurement technique over uncontaminated croplands (corn, soybean and rice paddy fields) (Baya and Van Heyst, 2010;Cobos and Baker, 2002;Kim et al., 2003;Cobbett and Van Heyst, 2007).The MM fluxes exhibited strong temporal variability during daytime and much weaker variability under low-quality turbulence during nighttime.In a typical campaign day, the turbulent flux data sets included both periods of emission and dry deposition.The median of nighttime flux was much smaller than the daytime flux for all MM methods (Mann-Whitney U test, MBR and AGM p<0.001, p<0.10 for REA).
The distribution of the turbulent fluxes and Hg 0 concentration gradient in Fig. 4 deviated significantly from a Gaussian distribution in the Hg 0 concentration gradient and in the derived MBR and AGM fluxes (Shapiro-Wilk's test rejected the hypothesis of normality of the distributions, p<0.01).The statistical MM fluxes (median ± MAD) in IC#1 (Fig. 7a) were −0.5 ± 8.9 and 0.1 ± 3.2 ng m −2 h −1 for AGM and MBR measurement, and 2.8 ± 29.0, 1.4 ± 15.2, and 8.8 ± 45.3 ng m −2 h −1 for AGM, MBR and REA in IC#2 (Fig. 7b), respectively.All the distributions of MM turbulent flux were associated with a positive kurtosis (3.8-16.2) and a slightly positive skewness (0.8-1.5).The observed flux frequency distributions for AGM and MBR peaked more strongly than that of REA (Fig. 7), with the MBR method giving the most confined distribution.Broader flux distribution measured by the REA sampling method has been reported in the measurements of turbulent fluxes for other gases (Fowler et al., 1995;Beverland et al., 1996;Nemitz et al., 2001).Previous studies suggest that vegetation canopy in the growing stage acts as an Hg 0 sink by net uptake of Hg 0 into foliage and therefore contributes to dry depositional flux (Bash and Miller, 2009;Stamenkovic and Gustin, 2007).However, the three MM techniques in this study derived significant higher average Hg 0 emission fluxes in IC #2 compared to IC#1, indicating that the vegetation sink strength was not sufficient to offset the efflux from underlying soil surface for croplands.Even though not measured, it is credible to assume that the soil Hg 0 efflux was higher during the warmer IC#2 due to higher temperature (Table 1) (Baya and Van Heyst, 2010;Gustin, 2011).

Comparison of Hg 0 fluxes derived from micrometeorological methods
The larger variability in REA-compared to the gradientderived fluxes is associated with a combination of methodological, instrumental and site-specific constraints influencing primarily the resolution of C REA (Eq. 3) as identified and discussed in Part II of this paper series (Zhu et al., 2015b).Nevertheless, a Friedman two-way analysis of variance by ranks (a non-parametric method) showed that the median fluxes by the three MM methods were not significantly different (χ 2 = 1.29 < χ 2 p=0.05 = 5.99).This indicated that AGM, MBR and REA methods produced comparable results with respect to the median location of Hg 0 turbulent flux during the inter-comparison.
The MBR method relies on scalar similarity (similarity in the scalar time series throughout the scalar spectra, Kaimal et al., 1972) between Hg 0 and temperature used as the proxy in this study.Since we have no means of explicitly characterizing Hg 0 scalar spectra, it is important to address the distribution of sources and sinks within the footprint area (Foken, 2008).By choosing a large flat and uniform fetch with confined Hg content in the soil substrate, significant divergence from scalar similarity between Hg 0 and temperature is less likely to occur.Nevertheless, non-stationary effects (e.g., advection of Hg polluted air-masses and related changes in concentration with time) bias the measured tur- cess (see Sect. 3.4).The MBR method becomes uncertain and may significantly overestimate flux when the numerator and denominator in the formula of eddy diffusivity approach small numbers, which typically occurs in periods at dawn, dusk and during nighttime (Eq.6, see Converse et al., 2010).As shown in Fig. 8, the 20-min-averaged AGM-and MBR-derived fluxes were well correlated during both campaigns (slopes of 0.76 and 0.86).However, when the sensible heat flux becomes small (small temperature gradient) at |H | < 20 w m −2 , the correlation coefficient diminishes drastically with a fall-off in slope (F AGM /F MBR = 0.35-0.36)implying that the MBR method can significant overestimate turbulent Hg 0 fluxes.MBR flux data collected in the presence of small scalar gradients (often during dawn and dusk transition periods) are therefore of questionable quality and should be considered for omission.AGM fluxes were on an average 26.1 % lower than MBR fluxes during IC #1, but 13.8 % higher during IC#2.The disparate results may largely stem from methodological issues (Fritsche et al., 2008b).In some previous studies using the AGM method to gauge various trace gas fluxes in-cluding Hg 0 (Edwards et al., 2001;Edwards et al., 2005;Simpson et al., 1997), normalization of Eq. ( 5) was introduced to mitigate for systematical failure of obtaining energy budget closures (Twine et al., 2000) by a factor of 1.3-1.35.The AGM method involves momentum flux, and an atmospheric stability parameterization in the flux calculation.For the conditions of weak developed turbulence to a greater extent prevailing under nocturnal stable stratification, where u * is very low, the AGM and MBR methods are prone to large uncertainties and corresponding fluxes are suggested to be flagged by applying wind or friction velocity thresholds (namely u * < 0.07-0.1 m s −1 ) (Fritsche et al., 2008b;Foken, 2008).During IC#2, when the REA system was included, the agreement between REA and the gradient-based methods was worst for small fluxes, which is inherently connected with the lower precision of the former system.As to be discussed in Zhu et al. (2015b), the non-constant (i.e., concentration and time dependent) sampling channel bias, which is difficult to entirely account for, is relatively more aggravating for the REA approach.For other gases (e.g., NH 3 , CH 4 ) that have been studied with this triad of MM techniques, higher  variability in REA flux is generically observed (Nemitz et al., 2001;Fowler et al., 1995;Moncrieff et al., 1998).In addition, systematic flux differences between a suite of NH 3 -REA systems as well as collocated AGM system inter-compared have been reported (Hensen et al., 2009).

Footprint of flux measurement
While the footprint (enclosed soil surface) of the chamber methods is fixed and very small (0.06 m 2 for TDFC and 0.09 m 2 for NDFC), MM methods derive fluxes from a footprint of comparatively large spatial extension upwind of the sampling tower.The MM footprint is not constant over time but a complex function of the sensor height, surface roughness length and canopy structure together with changing meteorological conditions.The predicted source area (using the models of Kljun et al., 2004 andKormann andMeixner, 2001) tends for upper sampling level (z 2 = z REA ) to be extensive for flux periods associated with weakly developed turbulence (Flag 2).In contrast, ∼ 70 % and ∼ 86 % of the data cleared for good turbulence quality, x70% (along-wind distance providing 70 % cumulative contribution to turbulent flux) fall within the unbroken field (150 m) for IC#1 and IC#2 respectively.For the lower sampling height (z 1 ), the footprint falls almost entirely within the primary fetch.Nevertheless, heterogeneous structures (roads, streams, tree stands and low buildings) existing outside the primary fetch (> 150 m) are of minor spatial extent, and within a radius of ∼ 2 km the sampling tower can be regarded to be surrounded by unbroken farmlands.

Diel variations
Figure 9 shows box and whisker plots of the diurnal variation of Hg 0 flux obtained by the five examined methods.Consistent in both campaigns, the MM methods exhibited highly variable fluxes, especially during daytime, where the magnitude in a single 20 min turbulent flux can exceed the flux derived by the chamber methods by many times.DFCs fluxes followed a well-defined diurnal pattern with consistent daytime emission and slight nighttime deposition.The pattern is similar to those for solar irradiance and temperature and reflects that the air-soil Hg 0 flux derived from the DFC technique is primarily governed by thermal and light-induced controls (e.g., Bahlmann et al., 2006).In contrast, flux from MM measurements is subject to the constant changes of atmospheric turbulence within the planetary boundary layer.
To facilitate a comparison between the DFC and MM data set on a diurnal basis, a Savitzky-Golay filter was applied on hourly averaged turbulent Hg 0 flux data to smooth out the short-term variability.In Fig. 10, where the diurnal courses of flux are given by smoothing spline fits, there is a 2 h lag in the time of the day when turbulent and chamber-derived flux peaked (IC#1).For the DFCs, the observed Hg 0 flux peaked within the period P2 (Fig. 10, IC#1) in concert with soil temperature, which is consistent with diurnal cycles reported for chamber measurements in the literature (Fu et al., 2008(Fu et al., , 2012;;Gustin, 2011;Zhu et al., 2013a).The smoothed mean diurnal cycle derived by the gradientbased methods over the same period exhibits peaking Hg 0 fluxes shortly before midday (P1 in Fig. 10, IC#1) but also includes a subsequent shoulder in the flux profile in the early afternoon (within P2 in Fig. 10, IC#1).The pattern resembles to an extent that of latent heat flux (evapotranspiration) (Liu and Foken, 2001) and may be interpreted as an effect of photo-reduction of previously deposited Hg II to Hg 0 into soil in conjunction with the presence of a water film (frost and dewfall) and emerging incoming solar radiation and temperature-driven air-surface exchange of soil Hg 0 pool (Fritsche et al., 2008b).Nevertheless, measurement of airsurface Hg 0 fluxes under the marked varying Hg 0 concentrations in air is challenging.Under such conditions, the measured turbulent fluxes are altered by non-stationary bias, and thus they do not represent actual fluxes to surface.The rates of change in Hg 0 concentration (up to ∼ ± 1.1 ng m −3 h −1 ) at the storage height of nearly 3 m relevant to this study imply vertical Hg 0 flux divergence in the range ± 3 ng m −2 h −1 .At low turbulence, advection in addition may as well gain some importance.However, to fully quantify the advection term for Hg 0 requires an array of instrumentation and such an investigation was not feasible in this study.
The mean diurnal cycles calculated for the three coevally examined MM methods (Fig. 10, IC#2) are based on a significantly smaller set of input data (∼ 30 % of IC#1) and therefore plausibly less robust to provide adequate representativeness after smoothing.Moreover, the campaign is a composite of periods where near-neutral conditions prevailed on daytime as well as adjacent nights and periods with weakly developed turbulence during nighttime respectively.Accordingly, the MM methods unanimously gauged maximum fluxes slightly after noon-time (P2, IC#2).However, there are features (P1 and P3) in the constructed cycles that are difficult to fully couple to environmental responses.
The cumulative flux derived by the examined methods is presented in Fig. 11a  rived a ∼ 25 % lower Hg 0 net evasion.This indicates that the flux correction with synchronized surface shear properties in NDFC partially bridges frequently observed disparities in magnitude between the MM-and conventional chamberderived fluxes (e.g., Gustin et al., 1999).Figure 12a shows the scatterplot of hourly flux specifically for MBR versus NDFC/TDFC -the correlation between individual hourly data points is weak.While in Fig. 12b, the deviation between MBR cumulative fluxes and NDFC/TDFC cumulative fluxes during the sampling campaign suggests that NDFC measurement shows a great advantage in bridging the flux gap between DFCs and MBR measurement.The significant scattering in Fig. 12a stems substantially from the inherent high variability in MBR flux prevalent during daytime.The difference between chamber and MBR flux depends to a certain degree on the diurnal variation of the atmospheric conditions.During daytime, the chamber produces a delay in the daytime flux evolution and fluxes become sustained in the late afternoon due to an artificial reduction in surface cooling within the chamber (Fig. 10).
During IC#2, the gradient-based MM techniques were evaluated together with the REA technique.The temporal features of the convoluted MBR and AGM cumulative fluxes are by and large concordant albeit the latter technique gauged ∼ 20 % higher Hg 0 net flux (1.78 vs. 1.43 µg m −2 ).The relative magnitude of MBR and AGM flux showed an inverse order during the two campaigns, possibly caused by method-  ological limitations given by the diverging micrometeorological conditions (Zhu et al., 2015b).For an extended period, the cumulative flux of REA given in Fig. 11b evolved in a similar way to those of the gradient-based methods (18-21 April).However, considerably different fluxes, occasionally in reverse directions, occurred after 21 April.In particular, during 16-17 April (Fig. 11b), a large net emission event was observed by all three techniques but at different magnitude.

Correlation between Hg 0 flux observation and environmental factors
It has been shown that the air-surface exchange of Hg can be influenced by solar irradiation, temperature, humidity, moisture, wind shear condition, and biotic processes (Choi and Holsen, 2009;Eckley et al., 2010;Fu et al., 2008;Gustin, 2011;Zhu et al., 2013a;Lin et al., 2010), as also observed in our field (Figs. 9 and 10).Table 2 shows the Pearson correlation coefficients between Hg 0 fluxes measured by the different methods and meteorological variables.DFC Hg 0 fluxes were positively correlated with solar radiation, soil temperature, soil moisture, friction velocity (R ∼ 0.4-0.9,p<0.05), and negatively correlated with air Hg 0 concentration and air humidity (p<0.05).The correlations between the MM fluxes and environmental variables were generally weaker (|R| < 0.5) in both campaigns.It is evident that DFC is less sensitive to surrounding atmospheric conditions that control the MM flux.In contrast, the Hg 0 flux controls in the ecosystem enclosed by the chamber are subject to microenvironment conditions that are significantly perturbed foremost by solar heating.The five techniques gauged unanimously positive net Hg 0 fluxes cumulated over the campaign periods.For the investigated triad of MM techniques, the Hg 0 -REA system has a general tendency to derive fluxes largest in magnitude.Over most of the campaign time, REA reported 20-60 % higher cumulative flux compared to the AGM method next to REA.Intriguingly, the Hg 0 flux budget magnitude examined by AGM and MBR methods was reversed during the two campaigns with a difference of ∼ 20 %, which may result from the atmospheric conditions and proxy scalar behavior.The traditional DFC method systematically measured the lowest Hg 0 net emission (42 % and 31 % of AGM-and MBR-derived net emission, respectively).The NDFC technique measured averaged fluxes similar to turbulent Hg 0 fluxes obtained by the MBR method (5.3 % difference).Although not entirely coupled to the atmospheric conditions that control the flux, the NDFC technique nevertheless represents a significant progress and improvement in contemporary enclosure-based Hg 0 flux measurement.
It was feasible to obtain a gradient measurement height ratio at the recommended bound (Foken, 2008).Given the lower precision of REA, gradient-based methods are consequently recommended for atmosphere-ecosystem Hg 0 flux measurements over low vegetation.REA has its niche over tall canopy, where gradient methods have frequently been found impracticable.In future applications, concerning foremost MM flux measurement technique, where the capacity to resolve small concentration differences is critical, it is recommended to implement analysis of synchronously collected samples for various heights (AGM, MBR) and conditionally segregated air parcels (REA) to avoid uncertainties induced by non-uniform ambient air Hg 0 concentration during the flux-averaging period.It has recently been argued that direct measurement of Hg 0 ecosystem air-canopy gas exchange is difficult and potentially subject to larger uncertainties (Zhang et al., 2012).Nevertheless, it is practicable for Hg 0 as it is for other trace gases and aerosols for which continuous MM flux measurement systems are key tools in ecosystem sciences.Our results show that improvement in resolving small Hg 0 concentration differences for the MM systems is required to further reduce uncertainties in the flux estimation.

Figure 1 .
Figure 1.Schematic illustration of the collocated MM and DFCs instrumentation set-ups.P, MFC and FM indicate a pressure transmitter, mass flow controller and flow meter of rotameter type respectively.

Figure 2 .
Figure 2. General meteorological parameters and ambient GEM concentration in the two campaigns.Upper panel: relative humidity (blue open circles), canopy leaf wetness (light blue line filled down), air temperature (red filled diamonds) and rainfall (black bar).Middle panel: wind speed (green line) and wind direction (dark green open circles filled down).Lower panel: ambient GEM concentration (dark purple open circles), global radiation (orange squares filled down) and σ w / u * (magenta line).

Figure 4 .
Figure 4. Time series of GEM gradients, GEM fluxes measured in: (a) IC#1 using MM and DFCs techniques; (b) IC#2 using MM techniques.The color code (green-yellow-red) denotes the quality (high-moderate-low) of turbulent flux data derived from general tests and black bars given in corresponding plots represent absolute flux uncertainties.

Figure 5 .
Figure 5. Distributions of Hg 0 flux derived from DFC measurements (upper panel: TDFC, lower panel: NDFC).The tripartite panels consists from left to right of a shadowgram (a suite of overlaid histograms with different bin widths), a box and whisker plot (the ends of the box represent Q1 and Q3 and the whiskers denote ± 1.5 times the interquartile range, IQR = Q3−Q1; sample points further away are given as individual markers) and the corresponding normal quantile plot (the unbroken solid line signifies the expected normal cumulative distribution and the dashed intervals the Lilliefors confidence bounds.The scale of the upper and lower abscissa indicates normal quantile and probability).Furthermore, in the box and whiskers plot, mean is indicated by a filled diamond while the median is the line within the box.The bracket outside of the box identifies the shortest half, which is the most dense 50 % of the observations.
Figure4aand b show the time series of normalized vertical Hg 0 concentration gradient (ng m −4 ) and Hg 0 flux (ng m −2 h −1 ) derived from the turbulent diffusion methods (MBR and AGM).Hg 0 concentration gradients were observed in the similar ranges of −0.49 to 0.33 and −0.48 to 0.25 ng m −4 in both campaigns (Table1and Fig.4), though the more occasionally shifting conditions of weak and developed turbulence in IC#1 tend towards promoting a higher scale of diurnal gradient variability (IC#1 vs. IC#2 standard deviation: 0.09 vs. 0.06).Our gradient observations are in alignment with measurement over temperate grasslands (−0.40 to 0.27 ng m −4 )(Fritsche et al., 2008b).Basic statistics of the MM Hg 0 flux observations is presented in Table1.The variability in our observations is similar with those reported from previous studies using the MM flux measurement technique over uncontaminated croplands (corn, soybean and rice paddy fields) (Baya andVan Heyst,  2010;Cobos and Baker, 2002;Kim et al., 2003;Cobbett and Van Heyst, 2007).The MM fluxes exhibited strong temporal variability during daytime and much weaker variability under low-quality turbulence during nighttime.In a typical campaign day, the turbulent flux data sets included both periods of emission and dry deposition.The median of nighttime flux was much smaller than the daytime flux for all MM methods (Mann-Whitney U test, MBR and AGM p<0.001, p<0.10 for REA).The distribution of the turbulent fluxes and Hg 0 concentration gradient in Fig.4deviated significantly from a Gaussian distribution in the Hg 0 concentration gradient and in the derived MBR and AGM fluxes (Shapiro-Wilk's test

Figure 7 .
Figure 7. Overview of the distributions of turbulent Hg 0 flux measured by the MM techniques: (a) IC#1, (b) IC#2.See Fig. 5 for a detailed description of the composite plots.

Figure 8 .
Figure 8. Scatterplots of 20 min MBR versus AGM flux during IC#1 (upper panel) and IC#2 (lower panel).The plots on the right-hand side depict specific data for which |H | < 20 w m −2 .

Figure 9 .
Figure 9.Diurnal variation of Hg 0 flux measured with various techniques represented as box and whisker plots.The two box horizontal border lines represent 25th and 75th percentiles from bottom to top, and whiskers indicate the 10th and 90th percentiles.Bold line and fine line in the box indicate mean and median flux.

Figure 10 .
Figure 10.Smoothed diurnal cycles of Hg 0 flux and Hg 0 concentration derived from hourly averaged input data.
, b.During IC#1, the cumulative fluxes measured by MBR and AGM fell between the fluxes measured by the two DFC methods.A period of divergence in the magnitude between the derived turbulent exchange parameters (eddy diffusivity of heat and υ tr ) resulted in intersected courses of MBR and AGM cumulative flux (17 November).MBR flux then stayed beyond the AGM flux on a cumulative basis for the rest of the campaign.The cumulative flux gauged by the TDFC method was the lowest (approximately 1/3 of MBR flux).Over the duration of IC#1, the net Hg 0 flux estimated by MBR and NDFC methods was in good agreement (2.90 vs. 3.02 µg m 2 ) while the AGM method de-

Figure 11 .
Figure 11.Time series cumulative Hg 0 flux using various techniques for (a) IC#1 over bare soil and (b) IC#2 over wheat canopy.

4
Conclusions and implicationsIn this study, we performed a comprehensive intercomparison of five contemporary Hg 0 flux quantification techniques through collocated measurements over an agricultural field.The flat terrain and homogeneous soil Hg content at the experimental site are ideal for the intercomparison of the DFC and MM techniques.MM-and DFC-derived Hg 0 fluxes showed distinct temporal characteristics.The former exhibited a highly dynamic variability while the latter had gradual temporal features.Diurnal trends showed that MM and DFC measurements diagnosed a similar daytime emission peak with different peaking times.Such differences were driven by separate sets of environmental factors influencing the DFC (irradiance and temperature) and MM (atmospheric turbulence properties) measurements.The three MM methods (REA, AGM and MBR) observed statistically significant, inseparable median Hg 0 fluxes (p<0.05)albeit REA flux was distributed over a much broader scale.Gradient and DFCs methods inter-compared favorably with respect to the confined location of median fluxes.Instantaneous fluxes measured by NDFC and TDFC and by MBR and AGM methods respectively were highly correlated (R>0.8,p<0.05) as the pairwise techniques are based on the same theoretical concept.However, the comparability between individual DFC and MM fluxes was poor to moderate (R ∼ 0.1-0.5)indicating the risk of utilizing sporadic (non-diurnally resolved) flux measurements as representative of an ecosystem.

Zhu et al.: Mercury flux data comparability and method characteristics 687 of
0.06 m 2 has been used extensively in our group and elsewhere were inter-compared.The hemicylindrical TDFC made of quartz with an open bottom area Atmos.Chem.Phys., 15, 685-702, 2015 www.atmos-chem-phys.net/15/685/2015/W.

Table 1 .
Summary of observed meteorological variables, Hg 0 concentrations, vertical Hg 0 concentration gradients and Hg 0 fluxes for two campaigns.

Results and discussion 3.1 Meteorological conditions
Meteorological observations and ambient Hg 0 concentration during the two campaigns are presented in Fig.2and summarized in Table1.The weather was predominantly sunny and temperate (−3.5 to 15.1 • during IC #1 and 0.8 to 17.4 • during IC#2).A rain shower yielding 3.4 mm precipitation occurred during IC#1.No precipitation was recorded during IC#2 (Fig.2upper panel).Leaf wetness and RH displayed clear diurnal variation (RH dropped to 40 % and leaf wetness to 0 % during daytime) except during the precipitation event when both were near saturation.Due to the high RH and sometimes sub-zero temperature at night, the ground and wheat possessed intermittently a light frost cover in early morning time.The wind speed was relatively high during daytime and turned moderate/calm at night.The wind direction was more variable from south to northeast with an average wind speed at 1.52 m s −1 (daytime mean: 1.98 m s −1 , nighttime mean: 1.05 m s −1 ) in IC#1, and changed to southwest and northeast with a mean of 2.69 m s −1 (daytime mean: 3.34 m s −1 , nighttime mean: 1.97 m s −1 ) in IC#2.The wind directions in IC#2 were more consistent than in IC#1:

Table 2 .
Pearson correlation analysis of hourly Hg 0 flux from various field measurement techniques and environmental parameters for two campaigns.Top-right segment of data are from IC#2.Bold font denotes a statistically significant correlation coefficient (p < 0.05).