Eddy covariance flux measurements of gaseous elemental mercury over a grassland

Direct measurements of the net ecosystem exchange (NEE) of gaseous elemental mercury (Hg⁰) are important to improve our understanding of global Hg cycling and, ultimately, human and wildlife Hg exposure. The lack of long-term, ecosystem-scale measurements causes large uncertainties in Hg⁰ flux estimates. It currently remains unclear whether terrestrial ecosystems are net sinks or sources of atmospheric Hg⁰. Here, we show a detailed validation of direct Hg⁰ flux measurements based on the eddy covariance technique (Eddy Mercury) using a Lumex RA-915 AM mercury monitor. The flux detection limit derived from a zero-flux experiment in the laboratory was 0.22 ng m⁻² h⁻¹ (maximum) with a 50 % cutoff at 0.074 ng m⁻² h⁻¹. We present eddy covariance NEE measurements of Hg⁰ over a low-Hg soil (41–75 ng Hg g⁻¹ in the topsoil, referring to a depth of 0–10 cm), conducted in summer 2018 at a managed grassland at the Swiss FluxNet site in Chamau, Switzerland (CH-Cha). The statistical estimate of the Hg⁰ flux detection limit under outdoor conditions at the site was 5.9 ng m⁻² h⁻¹ (50 % cutoff). We measured a net summertime emission over a period of 34 d with a median Hg⁰ flux of 2.5 ng m⁻² h⁻¹ (with a −0.6 to 7.4 ng m⁻² h⁻¹ range between the 25th and 75th percentiles). We observed a distinct diel cycle with higher median daytime fluxes (8.4 ng m⁻² h⁻¹) than nighttime fluxes (1.0 ng m⁻² h⁻¹). Drought stress during the measurement campaign in summer 2018 induced partial stomata closure of vegetation. Partial stomata closure led to a midday depression in CO₂ uptake, which did not recover during the afternoon. The median CO₂ flux was only 24 % of the median CO₂ flux measured during the same period in the previous year (2017). We suggest that partial stomata closure also dampened Hg⁰ uptake by vegetation, resulting in a NEE of Hg⁰ that was dominated by soil emission. Finally, we provide suggestions to further improve the precision and handling of the “Eddy Mercury” system in order to assure its suitability for long-term NEE measurements of Hg⁰ over natural background surfaces with low soil Hg concentrations (< 100 ng g⁻¹). With these improvements, Eddy Mercury has the potential to be integrated into global networks of micrometeorological tower sites (FluxNet) and to provide the long-term observations on terrestrial atmosphere Hg⁰ exchange necessary to validate regional and global mercury models.

1 Introduction

Mercury (Hg) is a top priority environmental pollutant that is transported through the atmosphere in its gaseous elemental form (Hg⁰; > 95 % of total atmospheric Hg). Anthropogenic Hg emissions into the atmosphere exceed natural emissions by approximately a factor of 5 (Outridge et al., 2018). Atmospheric Hg has a lifetime of 8–13 months (Saiz-Lopez et al., 2018); this allows for long-range transport before deposition back onto the Earth's surface, which also occurs at remote locations far from pollution sources. Once deposited, Hg can be transformed into methylmercury that can bioaccumulate in the freshwater and marine food webs, thereby posing a threat to human and ecosystem health (Watras et al., 1998; Fitzgerald et al., 2007; Mason et al., 2012; Braune et al., 2015).

Atmospheric Hg deposition to terrestrial surfaces occurs predominantly as Hg⁰ dry deposition through stomatal uptake by vegetation or as wet or dry deposition after oxidation in the atmosphere to more soluble reactive mercury, Hg(II) (Lindberg et al., 2007; Driscoll et al., 2013; Jiskra et al., 2018). Wet deposition of Hg(II) via rain and snowfall is relatively well quantified by Hg deposition networks such as the National Atmospheric Deposition Program (NADP), the European Monitoring and Evaluation Programme (EMEP) and the Asia-Pacific Mercury Monitoring Network (APMMN). Dry deposition of Hg(II) is difficult to measure, and its contribution to total Hg deposition remains uncertain (Gustin et al., 2013; Jaffe et al., 2014; Miller et al., 2018; Lyman et al., 2020). Mercury stable isotope fingerprints have identified Hg⁰ as the dominant deposition pathway to terrestrial surfaces. Dry deposition of Hg⁰ through vegetation uptake contributes 65 %–90 % of the total Hg deposited to soils (Demers et al., 2007; Jiskra et al., 2015; Enrico et al., 2016; Zhang et al., 2016; Zheng et al., 2016; Obrist et al., 2017). However, Hg⁰ dry deposition remains poorly constrained due to the lack of long-term monitoring networks (Obrist et al., 2018). Reduction of Hg(II) in terrestrial surface pools and the subsequent emission of Hg⁰ back to the atmosphere prolongs the cycling of anthropogenic Hg emissions in the environment and can thereby delay the effects of curbing primary anthropogenic emissions on human Hg exposure (Zhu et al., 2016; Wang et al., 2016; Obrist et al., 2018). Net ecosystem exchange (NEE) of Hg⁰, which refers to the balance between Hg⁰ dry deposition and emission from foliage and soils, represents a major factor in how fast the environment will recover from anthropogenic Hg pollution. On the global scale, estimates of the terrestrial NEE of Hg⁰ remain uncertain. In the most recent global mercury assessment, soil emission estimates were lowered to 1000 Mg a⁻¹ (UNEP, 2019) relative to 2200 Mg a⁻¹ in the 2013 assessment (UNEP, 2013); however, the associated uncertainties remain large. A recent review of 132 direct flux measurement studies revealed a NEE Hg⁰ flux between −513 and 1653 Mg a⁻¹ (the range of 37.5th and 62.5th percentiles, which refers to the central 25 % of the distribution) (Agnan et al., 2016). The database predominantly contains Hg⁰ flux measurements performed with dynamic flux chambers (85 % of all studies) that are ideal for short-term, mechanistic studies but are less suitable for quantitative flux estimations, especially over vegetated surfaces (Gustin et al., 1999; Eckley et al., 2016; Osterwalder et al., 2018). Year-round NEE measurements of Hg⁰ at the landscape scale are compelling to reduce measurement uncertainties. However, only four year-round whole-ecosystem Hg⁰ flux studies have been published, and all of these works used micrometeorological techniques, including the modified Bowen ratio and aerodynamic gradient methods (Fritsche et al., 2008a; Castro and Moore, 2016; Obrist et al., 2017), and the relaxed eddy accumulation (REA) technique (Osterwalder et al., 2017). These approaches use instruments that do not fulfill the fast Hg sensor response criterion that is required for eddy covariance (EC) flux measurements. Hence, these are not direct flux measurements and are consequently dependent on a number of assumptions. The main difficulty in using the modified Bowen ratio and aerodynamic gradient method is resolving a significant concentration gradient during turbulent conditions. During calm conditions, in contrast, it is challenging to determine a significant eddy diffusivity. Further drawbacks are (1) the potentially different sink/source characteristics of the footprint due to the two measurement heights, (2) temporally intermittent sampling between the two sampling inlets, and (3) the fact that transport characteristics are based on reference scalars like heat, water or CO₂ (Businger et al., 1986; Stannard et al., 1997; Edwards et al., 2005; Sommar et al., 2013a). The REA technique (Businger and Oncley, 1990) circumvents most of these difficulties. However, uncertainties in the Hg⁰ flux calculations are introduced by the determination of the proportionality coefficient (β-value) and system-dependent shortcomings, such as a biased offset between the updraft and downdraft sampling lines or difficulties in controlling the air flow from the air inlets to the analyzer. Thus, it remains challenging to accurately measure very small concentration differences with REA (typically < 0.1 ng m⁻³) between updrafts and downdrafts over natural surfaces with low substrate Hg concentrations (Cobos et al., 2002; Bash and Miller, 2008; Sommar et al., 2013b; Osterwalder et al., 2016, Kamp et al., 2018).

The EC technique has been under development since the late 1940s to measure the surface–atmosphere exchange of heat, mass and momentum in the surface boundary layer, which refers to the lowest 20–50 m of the atmosphere (Montgomery, 1948; Obukhov, 1951; Swinbank, 1951). In order to estimate a vertical turbulent flux, the covariance of two concurrently measured variables is calculated, (1) the scalar quantity of interest (in our case Hg⁰) and (2) the turbulent fluctuations of the vertical wind velocity, which are both measured at high temporal resolution. Since the 1990s, a new generation of digital three-axis ultrasonic anemometers, infrared gas analyzers and comprehensive software packages have facilitated land–atmosphere exchange measurements of CO₂ and H₂O (McMillen, 1988). Today, the EC technique is considered the standard method of determining evapotranspiration and the NEE of energy and trace gases such as CO₂, CH₄, N₂O, O₂, O₃ and volatile organic compounds using high-resolution (10–20 Hz), sometimes portable and generally very reliable equipment (Aubinet et al., 2012).

The first application of the EC technique to measure NEE of Hg⁰ reported an emission flux of 849 ng m⁻² h⁻¹ over contaminated soils (85 mg Hg kg⁻¹ dry soil) during a pilot campaign in Nevada, USA (Pierce et al., 2015). The EC system was based on a fast response (25 Hz), field-deployable pulsed cavity ring-down spectrometer (CRDS; Faïn et al., 2010; Pierce et al., 2013). However, the minimum detection limit of 32 ng m⁻² h⁻¹ did not allow for Hg⁰ flux measurements over soils exhibiting background Hg concentrations (typically < 100 ng Hg g⁻¹; Ericksen et al., 2006) (Pierce et al., 2015).

Here, we present EC measurements of the NEE of Hg⁰ over a grassland with typical background soil Hg concentrations. Our novel EC system makes use of a Lumex RA-915 AM mercury monitor (Lumex Ltd., St. Petersburg, Russia) atomic absorption spectrometer with Zeeman background correction, which allows for the measurement of Hg⁰ in ambient air at a relatively high sampling frequency of 1 Hz (Sholupov and Ganeyev, 1995; Sholupov et al., 1995). Ambient air Hg⁰ measurement comparison studies between the more frequently used Tekran^® 2537 analyzer (Tekran Inc., Toronto, Canada) and the RA-915 AM were performed by the European Committee for Standardization's (CEN) Technical Committee 264 “Air Quality” EN 15852 and showed good agreement between the two instruments (Brown et al., 2010). Among other applications, the mercury monitor's precursor, the Lumex RA-915+ mercury analyzer was successfully deployed in the Global Mercury Observation System (GMOS) project at two sites in Russia and Suriname (Sprovieri et al., 2016).

The objective of this study was to test the performance of the RA-915 AM as a fast response analyzer as well as its suitability for EC flux measurements with the goal of reliably measuring the NEE of Hg⁰ over terrestrial ecosystems. Hereinafter, the new EC system is referred to as “Eddy Mercury”. We provide a description of the Eddy Mercury system and present the data analysis procedure to calculate the NEE of Hg⁰ in detail. We discuss the patterns in the NEE of Hg⁰ measured over a grassland during a 34 d pilot campaign and provide suggestions to improve the reliability and precision of the Eddy Mercury system for future long-term applications.

Figure 1 (a) Schematic of the experimental setup at the Chamau (CH-Cha) research site showing the exact locations of topsoil samples for total Hg analysis (n=22) and eddy covariance (EC) flux measurements of Hg⁰, CO₂ and H₂O conducted between 20 July and 6 September 2018. Panel (b) shows the footprint contour lines of 10 % to 90 % in 10 % steps representing the flux source area during our measurement period. Numbers indicate the distance in meters from the EC station (black cross). The footprint was calculated by applying the footprint model presented in Kljun et al. (2015). Figure (b) is a direct output from the following online tool: http://footprint.kljun.net/ (last access: 6 February 2020).

2 Material and methods

2.1 Site description and instrumentation

The Eddy Mercury system was tested between 20 July and 6 September 2018 at the Swiss FluxNet site Chamau (CH-Cha), which is located in central Switzerland, about 30 km southwest of Zurich (47^∘12^′36.8^′′ N, 8^∘24^′37.6^′′ E; 393 m a.s.l.). In this study, NEE of Hg⁰ and CO₂ was measured concurrently with two independent EC systems over the intensively managed grassland used for forage production. Details on grassland species composition, harvest, and fertilization practices are described in Zeeman et al. (2010), Merbold et al. (2014) and Fuchs et al. (2018). The tower for long-term EC greenhouse gas measurements was located between two adjacent grassland parcels (Fig. 1a). The northern parcel, measured when up-valley winds prevailed, was oversown with clover in March 2015 and April 2016 to investigate the N₂O emission reduction potential compared with the conventionally fertilized grassland of the southern parcel, measured primarily when down-valley winds prevailed (Fig. 1b). The soil type is a Gleysol–Cambisol, with a bulk density of about 1 g cm⁻³, 30.6 % sand, 47.7 % silt and 21.7 % clay in the top 10 cm (Roth, 2006). A topsoil pH of 5.3 was determined by adding 25 mL of 0.01 M CaCl₂ solution to 10 g of dry soil (Labor Ins AG, Kerzers, Switzerland) in 2014. The 24-year (1994–2017) average annual temperature measured at the nearby SwissMetNet surface weather station in Cham (CHZ, 444.5 m a.s.l.) was 10.1 ^∘C, and the average annual precipitation was 997 mm.

The Eddy Mercury system was mounted approximately 3 m west of a fully equipped long-term EC tower measuring greenhouse gas exchange (CO₂ and H₂O) and meteorological variables at a height of 2 m (Fig. 1). The CO₂ flux system consisted of a three-dimensional (3-D) ultrasonic anemometer (Solent R3-50, Gill Instruments, Lymington, UK) and an open-path infrared gas analyzer for CO₂ and H₂O concentrations running at a 20 Hz resolution (IRGA, LI-7500, LI-COR Biosciences, Lincoln, NE, USA). From the 20 Hz IRGA measurements, 30 min flux averages were calculated using the LI-COR EddyPro^® software. The 30 min CO₂ flux has been recorded continuously since 2005 (Eugster and Zeeman, 2006; Zeeman et al., 2010). The measured meteorological variables included temperature and relative humidity (Hygroclip S3 sensor, Rotronic AG, Switzerland), net all-wave radiation (CNR 1, Kipp & Zonen B.V., Delft, Netherlands), incoming and reflected photosynthetic active radiation (PAR lite, Kipp & Zonen B.V., Delft, Netherlands), and precipitation (height of 0.5 m; tipping bucket rain gauge from LAMBRECHT meteo GmbH, Göttingen, Germany). In addition, soil temperature was recorded at depths of 0.05, 0.1, 0.15, 0.25 and 0.4 m (T107, Campbell Scientific Inc., Logan, UT, USA).

2.2 Soil sampling and total mercury analysis

Topsoil samples (0–10 cm) were taken in a circular arrangement around the EC tower (Fig. 1a) using a core drill. The soil samples were transported to the laboratory in sealed plastic bags and stored in a fridge at 4 ^∘C. The samples were filled into aluminum shells, weighed and dried at 40 ^∘C until their weight remained constant. The samples were pestled and sieved through a 2 mm mesh to separate the fine earth and the skeleton. The fine earth was ground to powder using a laboratory scale ball mill. To get rid of all potential humidity, the ground samples were stored in small paper bags in a desiccator and dried again at 40 ^∘C. The 22 topsoil samples were analyzed for total Hg using a DMA-80 direct mercury analyzer (MLS Mikrowellen GmbH, Leutkirch im Allgäu, Germany). Certified Hg standard solution (NIST-3133) was gravimetrically diluted to concentrations of 10 to 1000 ng g⁻¹ and used for the calibration of the instrument. Repeated measurements of standard reference material (ERM-CC141 loam soil) 90.3±7.8 ng g⁻¹ (mean ± standard deviation, n=3) agreed with the certified value (83±17 ng g⁻¹).

2.3 Description of the Eddy Mercury system

The core of the Eddy Mercury system with respect to the measurement of the NEE of Hg⁰ is the RA-915 AM mercury monitor. The RA-915 AM uses atomic absorption spectrometry (AAS) with Zeeman background correction to continuously measure Hg⁰ in ambient air (Sholupov et al., 2004). The multipath sample cell of the RA-915 AM has an optical path length of 9.6 m and a cell volume of 0.7 L. Baseline corrections (zero drift) were performed automatically by the instrument using Hg-free air at user-defined intervals. Span corrections are carried out using an inbuilt calibration cell that contains Hg⁰ vapor. The measurement range lies between 0 and 2000 ng m⁻³, and the instrument detection limit is 0.5 ng m⁻³ according to the analytical specifications of the manufacturer. The air flow rate was increased to 14.3 L min⁻¹ by bypassing the instrument pump in order to reduce the residence time in the measurement cell (normal flow of 7 L min⁻¹). For this, a stronger external pump was connected (model MAA-V109-MD, GAST Manufacturing, MI, USA). The instrument was placed in a weatherproof, air-conditioned box (Elcase, Marthalen, Switzerland) to protect the sensitive RA-915 AM from rain and to reduce temperature fluctuations. A USB to RS232 serial data interface was used to establish a one-way communication link from the RA-915 AM to the data acquisition computer. The air inlet was mounted 24 cm below the center of the head of the 3-D ultrasonic anemometer (Gill R2A, Solent, UK) used for wind vector measurements that was installed 2 m above the ground. A micro-quartz fiber filter (Grade MK 360, 47 mm diameter, Ahlstrom–Munksjö, Sweden) was installed in a 47 mm perfluoroalkyl polymer (PFA) single-stage filter assembly (Savillex, Eden Prairie, USA) at the air inlet. The air inlet was connected to the RA-915 AM by a 2.8 m intake hose with a 11 mm inner diameter (i.d.) attached to a 0.35 m, 4 mm i.d. sample intake hose. Both hose segments were unheated, insulated PFA tubing. The median lag time of the turbulent airflow (Reynolds number of > 5000) from the tube inlet to the analyzer was in the order of 1.15 s.

2.4 Eddy covariance flux measurements

The RA-915 AM analyzer was configured to measure Hg⁰ concentrations at 1 Hz. The Hg⁰ concentrations and the 3-D wind vectors were measured from 20 July to 6 September 2018 using four different settings of the RA-915 AM analyzer with respect to the length of the measurement interval between two autocalibration cycles (zero and span): (1) 24 h intervals from 20 to 26 July 2018; (2) 4 h intervals from 1 to 26 August 2018; (3) 1 h intervals from 27 to 31 August 2018; (4) 4 min intervals from 31 August to 6 September 2018. The ultrasonic anemometer had an internal sampling frequency of 1000 Hz that was averaged (eight records of each acoustic sensor pair for each direction) to 20.83 Hz. The 1 Hz RA-915 AM data were merged with the ultrasonic anemometer's data stream by oversampling, as described in Eugster and Plüss (2010). Data were collected on a Linux-based Raspberry Pi computer equipped with a real-time clock chip and internet access. Because data transfer via the USB port from the embedded Windows 7 system of the RA-915AM was highly unreliable, only the system time stamps were synchronized with the Linux data acquisition system every second via a Windows PowerShell script. In cases when this communication also failed, an approximate time synchronization was done by polling the RA-915 AM time stamp via the Samba file-sharing protocol. Thus, in addition to the synchronization method described by Eugster and Plüss (2010), the merging of Hg⁰ measurements with wind vector data had to be carried out off-line in a separate data workup step. Fluxes were calculated over 60 min intervals to account for the low sampling frequency of Hg⁰ signals. Thus, using the previously mentioned interval modes (1) and (2), 3600 Hg⁰ measurements were used for each 1 h flux average.

2.5 Eddy covariance Hg⁰ flux calculations

Calculation of the NEE of Hg⁰ required some modifications of the standard procedure that has been established for CO₂ fluxes (e.g., Aubinet et al., 2012). The modifications were carried out according to the five steps described in detail below (Sect. 2.5.1–2.5.5).

Figure 2 Examples of raw data time series over a 1 h data segment: (a, c, e) during period 1 with 24 h calibrations only (21 July 2018, 10:00–11:00), and (b, d, f) during period 4 with 4 min instrument calibration intervals (6 September 2018, 01:00–02:00). Vertical wind speed (a, b) was not detrended. Hg⁰ concentrations are shown before detrending (c, d) and after detrending (e, f). While the 4 min calibration intervals clearly reduce the longer-term drift (d) compared with daily calibrations (c), the gaps during calibrations had to be filled using linear interpolation before calculating fluxes.

2.5.1 Preparation of raw Hg⁰ measurements

The RA-915 AM raw data files provide the following information at a 1 Hz resolution: date and time of measurement, photomultiplier current (arb. unit), air flow rate (L min⁻¹), temperature of analyzed air (^∘C), temperature of RA-915 AM (^∘C), sample cell pressure (kPa), Hg⁰ raw concentration (ng m⁻³, including all online corrections), status code and status description. The status code (a numerical value) and status description (a text variable) are redundant and provide the necessary information to distinguish ambient air concentration measurements from zero and span calibration measurements. The Hg⁰ flux was calculated based on the Hg⁰ raw concentration. To account for drift and baseline drift, which both are unavoidable when longer measurement periods are used between calibration events, we proceeded as follows. After a calibration event, the Hg⁰ raw concentration was considered to be the best empirical estimate of the true Hg⁰ concentration. By the end of a measurement period (the beginning of the next calibration cycle), a linear drift correction was applied to bring the Hg⁰ raw concentration before the next calibration event to the level of the next calibration result (offset correction). As visual inspection of the data clearly indicated that there was more drift than a simple linear trend in the data (see examples in Fig. 2), a high-pass filter approach was used to minimize drift and optimize the determination of Hg⁰ fluctuations for EC flux measurements (Sect. 2.5.4).

2.5.2 Preparation of the ultrasonic anemometer data

The ultrasonic anemometer data contained the three wind speed components of the wind vector (all in m s⁻¹), the speed of sound (m s⁻¹) and the information sent from the RA-915 AM to the data acquisition system via the serial data link. The speed of sound (c) was converted to a virtual sonic temperature ($T_{\mathrm{v}}\approx c^{\mathrm{2}}/\mathrm{403}$) in Kelvin (Kaimal and Gaynor, 1991). The vertical wind speed (w) was despiked using an iterative 7σ filter that discards w outside the range of the 6 h mean ±7 standard deviations.

2.5.3 Merging of ultrasonic anemometer data with the Hg⁰ time series

After the preparation of the two datasets, they were merged by accounting for the time difference between the RA-915 AM and the Linux data acquisition using the information that could be transferred via the serial link from the RA-915 AM to the Linux system (accurate to within 1 s). If no such information was received from the RA-915 AM, the time difference between the two systems was determined using a network time drift fallback option specifically added to the Linux system to overcome the problems with serial output from the RA-915 AM: during the field experiment we polled the most recent data record acquired by the RA-915 AM every 5 min using the Samba file-sharing protocol, and we associated that time stamp with the one of the EC system. This (somewhat less accurate) information was then adjusted during periods where both approaches overlapped to determine the time difference required to shift the Hg⁰ raw data relative to the ultrasonic anemometer data before merging the two datasets. To ascertain that Hg⁰ data were lagging the sonic data, we added a ≈1.5 s safety margin in the interpretation of the available time synchronization information received either via serial link or Samba file-sharing.

2.5.4 Determination of the time lag between vertical wind speed and Hg⁰ fluctuations

The merged dataset was then divided into 1 h segments for Hg⁰ flux calculations. Within each 1 h segment, the time lag between the two time series was fine-tuned using a cross-correlation procedure to find the best positive or negative correlation within a reasonable time window (0–4 s) around the physically expected time difference (the 1.15 s physical delay plus the 1.5 s safety margin used in Sect. 2.5.3). Because considerable nonturbulent drift of the Hg⁰ signal was still present after correcting for online calibration (Sect. 2.5.1), we detrended each 1 h segment using a third-order polynomial fit (Eq. 5) before computing the cross-covariance between the detrended Hg⁰ signal and w (Sect. 3.2.1). To account for the different sampling rates of w (20.83 Hz) and Hg⁰ (1 Hz), we used simple linear interpolation between individual Hg⁰ measurements and to bridge across calibration gaps. After a first automatic run, each best estimate for time lag was visually inspected and updated by a narrower search window for each 1 h segment that narrowed the search procedure down to the most realistic cross-correlation peak (positive or negative). Note that calibration gaps are relevant data gaps with setting 4 (from Sect. 2.4) but less problematic with settings 1–3. In all cases, the lack of variance in Hg⁰ data during the gaps reduces the computed Hg⁰ flux. Thus, our flux estimates are conservative estimates with respect to flux magnitudes.

2.5.5 Computation of Hg⁰ EC fluxes

After all data preparations according to Sect. 2.5.1 to 2.5.4, the Hg⁰ flux ($F_{{\mathrm{Hg}}^{\mathrm{0}}}$) was calculated as the covariance

$$\begin{array}{}\text{(1)}& F_{{\mathrm{Hg}}^{\mathrm{0}}}=\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\chi }}^{\prime }},\end{array}$$

where χ was the calibrated, detrended and linearly gap-filled Hg⁰ concentration in nanograms per cubic meter (ng m⁻³), and w was the vertical wind speed. For improved readability $F_{{\mathrm{Hg}}^{\mathrm{0}}}$ was converted from nanograms per square meter per second to nanograms per square meter per hour (ng m⁻² s⁻¹ to ng m⁻² h⁻¹) before reporting. In the notation used here, primes denote short-term deviations from the mean (after detrending according to Sect. 2.5.4) over an averaging period (1 h), and overbars denote the mean of a variable. Hg⁰ flux computations were carried out using R version 3.5.2 (R Core Team, 2018).

2.6 Determination of the Hg⁰ flux detection limit

To determine whether a calculated Hg⁰ flux is significantly different from a zero flux, we used two approaches: (1) an indoor zero-flux experiment, and (2) a statistical estimate of the flux detection limit following the concept by Eugster and Merbold (2015) that is an improvement of the concept presented by Eugster et al. (2007). The indoor zero-flux experiment was set up in the laboratory for 2 d before all of the equipment was installed in the field. The low-turbulence conditions in combination with the absence of local Hg⁰ sources in the laboratory allowed us to ascertain the fluxes that result from the procedure described above when there is no real Hg⁰ flux. Such zero-flux experiments tend to underestimate the flux detection limit under real-world outdoor conditions, while the second approach quantifies the statistical uncertainty of a calculated flux. The flux (covariance) is the product of the correlation coefficient (r_w,χ) between w and χ and the square root of the variances of the two variables (e.g., Eugster and Merbold, 2015):

$$\begin{array}{}\text{(2)}& \stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\chi }}^{\prime }}=r_{w,\mathit{\chi }}\cdot \sqrt{{w^{\prime }}^{\mathrm{2}}}\cdot \sqrt{{{\mathit{\chi }}^{\prime }}^{\mathrm{2}}}=r_{w,\mathit{\chi }}\cdot {\mathit{\sigma }}_w\cdot {\mathit{\sigma }}_{\mathit{\chi }}.\end{array}$$

The significance of r_w,χ can be estimated using a Student's t test (see Eugster and Merbold, 2015, for details). Thus, for each 1 h period, we computed the value of r_w,χ that is significant at p=0.05, and we multiplied this value by measured σ_w and σ_χ to obtain a more realistic estimate of the flux detection limit. It should be noted that this concept was brought forward by Wienhold et al. (1996) using a visual empirical approach, and Eugster and Merbold (2015) then further developed the visual approach into a more objective time series statistical approach to perform the quantification of the flux detection limit. The threshold of significance of r_w,χ can be estimated as follows:

$$\begin{array}{}\text{(3)}& r_{w,{\mathit{\chi }}_p}={\displaystyle \frac{t_p}{\sqrt{n-\mathrm{2}+t_p^{\mathrm{2}}}}},\end{array}$$

where t_p is the Student's t value for the significance level p (e.g., 0.05), and n is the autocorrelation-corrected number of independent samples in the time series,

$$\begin{array}{}\text{(4)}& n\simeq N{\displaystyle \frac{\mathrm{1}-{\mathit{\rho }}_{\mathrm{1}}}{\mathrm{1}+{\mathit{\rho }}_{\mathrm{1}}}},\end{array}$$

where N is the number of samples in a time series, and ρ₁ is the lag 1 autocorrelation coefficient of the scalar product time series w⋅χ.

2.7 Eddy covariance CO₂ flux calculations and quality control flags

The 30 min CO₂ flux was quantified in the conventional way established in ecosystem studies (see Aubinet et al., 2012) using the EddyPro^® (LI-COR Inc., Lincoln, NE, USA) software (see Fuchs et al., 2018, for specific information related to the Chamau field site). For each 30 min CO₂ flux interval, a flux quality control (QC) flag was determined following Mauder and Foken (2004): 0 (best data quality for detailed investigations), 1 (good data for longer-term studies) and 2 (poor quality). As there are no established quality control procedures for Hg⁰ fluxes yet, we used the QC information from the CO₂ flux measurement to retain or reject concurrent Hg⁰ flux measurements. Thus, we only present Hg⁰ flux measurements with CO₂ flux quality flags < 2. During CO₂ flux processing using the EddyPro^® software, coordinate rotation for tilt correction, angle of attack correction for wind components, Webb–Pearman–Leuning terms for compensation of density fluctuations (Webb et al., 1980), and analytical corrections for high-pass (Eugster and Senn, 1995; Moncrieff et al., 2004) and low-pass filtering effects (Horst, 1997) were applied. Furthermore, a self-heating correction for the open-path gas analyzer was conducted (Burba et al., 2008), and CO₂ fluxes greater than 50 and less than −50 µmol m⁻² s⁻¹ were discarded.

3 Results and discussion

3.1 Environmental conditions

In 2018, the annual mean air temperatures in Switzerland reached 6.9 ^∘C, which was the highest value recorded since the onset of meteorological measurements in 1864 (MeteoSchweiz, 2019). This nationwide average temperature was 1.5 ^∘C warmer than the average of the normal period from 1981 to 2010. Total precipitation measured from April to November 2018 was only 69 % of the long-term average (1981–2010). Thus, the period from April to November 2018 was the third driest period ever recorded in Switzerland (MeteoSchweiz, 2019). From the beginning of the growing season until the end of our measurement campaign (April to September 2018), air temperatures at the Cham (CHZ) SwissMetNet surface weather station were elevated by 2.2 ^∘C compared with the long-term average from 1994 to 2017 (15.8 ^∘C) during the same period. Total precipitation from April to September 2018 was only 72 % (467 mm) of the long-term average (648 mm) calculated for the period between 1994 and 2017. These specific conditions reduced CO₂ uptake compared with the same period in 2017 (Sect. 3.3, Fig. 7) and led to lower grassland productivity and yields of only 6.8 t DM ha⁻¹ a⁻¹ (DM refers to dry matter) in 2018 compared with an average yield of 12.7 t DM ha⁻¹ a⁻¹ quantified from 2015 to 2017 (start of the clover experiment). Over the course of the 34 d campaign (from 20 July 2018 at 02:00 to 24 July 2018 at 08:00 and from 9 August 2018 at 12:00 to 6 September 2018 at 17:00; all times in central European time, CET, or UTC+1), sunny conditions prevailed with a mean solar irradiation (Rg) of 352 W m⁻² during daytime (Rg ≥5 W m⁻²) and a mean irradiation of 606 W m⁻² at 13:00. The hourly mean air and soil surface temperature ranged from 13.6 ^∘C (06:00) to 24.1 ^∘C (15:00) and from 18.1 ^∘C (08:00) to 21.5 ^∘C (18:00), respectively. The median daytime (Rg ≥5 W m⁻²) and nighttime (Rg < 5 W m⁻²) wind speed was 0.97 m s⁻¹ (range of 0.05–5.77 m s⁻¹) and 0.37 m s⁻¹ (range of 0.06–2.49 m s⁻¹), respectively. The prevailing wind direction during the day was north–northwest (47 %), whereas it was east–southeast (55 %) at night.

3.2 Performance of the Eddy Mercury system

3.2.1 High-frequency signal analysis

Two examples of the raw data used to compute fluxes (Eq. 1) are shown in Fig. 2, one from period 1 with 24 h calibration intervals (Fig. 2a, c, e) and one with frequent calibrations every 4 min (Fig. 2b, d, f). Frequent calibrations strongly reduced the instrument drift (Fig. 2d) compared with the long calibration intervals (Fig. 2c), although at the expense of some loss of variance and flux (as will be discussed below). In principle, block-averaging raw data within a sampling interval is the best approach to compute EC fluxes (Aubinet et al., 2012). In the case of substantial instrument drift, as is seen with the RA-915 AM (Fig. 2c), it is necessary to remove the drift using an adequate procedure. Because of the curvature of the drift of the analyzer, a simple linear detrending did not lead to satisfactory results in this case; hence, we used a third-order polynomial regression fit:

$$\begin{array}{}\text{(5)}& {\mathit{\chi }}^{\prime }=\mathit{\chi }+{\mathit{\alpha }}_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{1}}\cdot t+{\mathit{\alpha }}_{\mathrm{2}}\cdot t^{\mathrm{2}}+{\mathit{\alpha }}_{\mathrm{3}}\cdot t^{\mathrm{3}},\end{array}$$

with t elapsed time within the averaging interval of 1 h. The turbulent Hg⁰ fluctuations after this additional detrending led to the time series shown in Fig. 2e and f. Lengthy discussions on the possible shortcomings of such a detrending can be found in Lee et al. (2005) and Aubinet et al. (2012) and are therefore not repeated here. Using the example data shown in Fig. 2e, we produced an artificial dataset with gaps that correspond to the 4 min recalibration scheme used during the period shown in Fig. 2f. This led to a loss in Hg⁰ flux in the order of 12 %. Although nonzero, this should be considered a robust finding given the general understanding that EC flux measurements are accurate to within 10 %–20 %, even with higher quality instrumentation (Aubinet et al., 2012). To obtain higher quality EC fluxes than we present here, the long-term stability of the instrument needs to be improved (Sect. 3.4); however, improving the gap-filling strategy is not expected to contribute significant new insights into Hg⁰ flux calculations. Drift of the current version of the Eddy Mercury system is substantial (Fig. 3a), which is an effect that is common in experimental sensor setups, although it is no longer prevalent in present-day CO₂ sensors. The removal of any drift also reduces the variance of a signal and, subsequently, the flux covariance of interest. Thus, knowledge about the stability of an instrument over which no drift correction is required becomes important. The Allan variance plot (Fig. 3b; Allan, 1966; Werle et al., 1993) indicates that the optimum averaging time is approximately 54 s. For comparison, a CH₄ analyzer tested by one of the authors (Eugster and Plüss, 2010) shows an optimum average time that is roughly 3 times as long (approximately 180 s) before the instrument drift starts to dominate the Allan variance. Figure 3b shows that the Allan variance caused by drift at integration times beyond 550 s exceeds the variance associated with turbulence at the 1 s integration time (see blue arrow in Fig. 3b). In a more ideal instrument, the long-term drift is smaller than the short-term variance of interest for EC measurements (see, e.g., Eugster and Plüss, 2010). Despite these findings, Fig. 3 clearly shows the potential and quality of the instrument for Hg⁰ flux measurements.

Figure 3 Allan variance plot using 14 h of continuous measurements in the laboratory (zero-flux experiment), starting on 18 July 2018 at 18:00. Panel (a) shows the raw time series, and panel (b) shows the Allan variance as a function of integration time.

Figure 4 Panel (a) shows an example spectrum of Hg⁰ fluctuation measurements, and panel (b) displays a cospectrum of the 1 h averaged Hg⁰ flux on 31 August 2018 from 14:00 to 15:00.  (a) The power spectrum before (red line) and after (black line) detrending is shown, and the theoretical slopes in the inertial subrange are shown for ideal conditions (i, solid line, $f^{-\mathrm{2}/\mathrm{3}}$ slope), for a rectangular oversampling at frequencies greater than 1 Hz (r, broken line, f⁻¹ slope) and for a first-order damped spectrum (d, dashed line, $f^{-\mathrm{8}/\mathrm{3}}$ slope). The approximated white noise level is shown using a color band (w, f⁺¹ slope). (b) The flux cospectrum shows the absolute values of the cospectral densities, with black symbols denoting positive contributions to $\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\chi }}^{\prime }}$ and red symbols denoting negative contributions. The light green bold line is a local polynomial regression fit to the data points, and the blue line denotes an idealized cospectrum. The two dashed blue lines show damped cospectra with respective damping constants of 0.1 and 0.3 s.

This interpretation is also supported by spectral and cospectral analyses (Fig. 4). Figure 4a shows an example spectrum of Hg⁰ measurements obtained over a 1 h interval. The difference between the red and black lines in Fig. 4a visualizes the effect of polynomial detrending on the power spectrum of Hg⁰, which is relatively small and of no real concern. As the RA-915 AM only delivers 1 Hz raw data, we had to oversample this digital Hg⁰ signal to match the 20.83 Hz resolution of the ultrasonic anemometer. Spectral densities at high frequencies > 0.5 Hz (the Nyquist frequency of the RA-915 AM is half of the sampling frequency) reflect the effect of oversampling. In the case of the RA-915 AM, oversampling leads to local minima in spectral densities at 1 Hz and all its harmonic multiples (2, 3, 4, … Hz), which is the result of linear interpolation between measurements. Between these local minima, the spectral density obeys the f⁻¹ power law (line “r” in Fig. 4a), which is very close to the inertial subrange slope $f^{-\mathrm{2}/\mathrm{3}}$ (line “i” in Fig. 4a). A damped signal (first-order damping; see Eugster and Senn, 1995) would follow a $f^{-\mathrm{8}/\mathrm{3}}$ power law (line “d” in Fig. 4a); thus, it is obvious that our setup had an adequate flow rate through the RA-915 AM that did not lead to substantial damping of the turbulent Hg⁰ fluctuations. With the oversampling used here, the white noise level (blue band “w” in Fig. 4a) is artificially reduced below the level that we would obtain without oversampling.

After an adequate time lag correction to synchronize the detrended Hg⁰ signal with vertical wind speed fluctuations w^′ was applied, the cospectra of fluxes that were significantly different from a random pattern closely agreed with the theoretical idealized cospectrum for neutral atmospheric stability derived from Kaimal et al. (1972) (see Eugster and Senn, 1995), shown using the solid blue line in Fig. 4b. Some minor signs of damping are seen at higher frequencies where the green spline deviates from the solid blue line (Fig. 4b). The comparison of cospectral densities with theoretical damped cospectra (dashed blue lines in Fig. 4b) clearly confirm the finding from the spectral analysis that the flow rate was high enough in the RA-915 AM sample cell to prevent significant damping effects that tend to be a problem with closed-path EC flux measurements.

Figure 5 Cross-correlation analysis to determine the time lag between the vertical wind speed (w) and Hg⁰ time series (χ). (a, b) An example with a clearly positive Hg⁰ flux (21 July 2018, 10:00–11:00), and (c, d) an example with a marginally positive flux (6 September 2018, 01:00–02:00). Panels (a) and (c) show the cross-correlation within a time lag window of ±10 s, and panels (b) and (d) zoom in to the search window used in this study (vertical red band). The blue horizontal band shows the range of zero fluxes (cross-correlation $r_{w,\mathit{\chi }}\ne \mathrm{0}$ with p≥0.05).

When a clear Hg⁰ flux was statistically different from a zero flux, the cross-correlation peak was well defined (Fig. 5a, b). On some occasions with low fluxes relative to the flux detection limit (Sect. 3.2.2), the automatic detection of the cross-correlation peak was not successful. The peak often does not extend very strongly beyond the (expected) noise level, as shown in Fig. 5c. However, when zooming in (Fig. 5d), the peak becomes rather clear, although it is only marginally above the range of insignificant correlations shown by the blue background in Fig. 5. To minimize erroneous peak detections and, thus, incorrect flux estimates, we fine-tuned the search window (red band in Fig. 5) for each 1 h data segment by visually inspecting and selecting the search window within which the local maximum of the absolute correlation coefficient between w and χ was found.

3.2.2 Flux detection limit

The flux detection limit was calculated for each 1 h flux period (Sect. 2.6). The significance threshold for r_w,χ was calculated for an error probability p=0.05, and the product of this threshold r_w,χ and measured σ_w and σ_χ was determined as the flux detection limit for that specific 1 h period. Figure 6 shows the probability density function of the flux detection limits from all 1 h data segments. For comparison, the results from the 14 h zero-flux experiments in the laboratory are added as a blue boxplot in Fig. 6. This comparison clearly shows that a zero-flux experiment in the laboratory highly overestimated the quality of Hg⁰ flux measurements with a median (maximum) flux detection limit of 0.074 (0.22) ng m⁻² h⁻¹. The more realistic flux detection limits based on statistically significant (p < 0.05) correlations are rather in the order of 5.9 (50 % cutoff) to 24 ng m⁻² h⁻¹ (99 % cutoff) with a 95 % cutoff at 13.7 ng m⁻² h⁻¹. During the 34 d measurement campaign, 49.7 % of the Hg⁰ fluxes (363 out of 731 h) were significantly different from zero. Using the same approach but in a qualitative way, Pierce et al. (2015) estimated the flux detection limit of their system to be around 32 ng m⁻² h⁻¹.

3.2.3 Comparison of detection limits for Eddy Mercury, gradient-based and REA systems

The Eddy Mercury system circumvents major sources of uncertainty compared with gradient-based and REA systems, which are related to assumptions on similarity or equivalence of the eddy diffusivities of the scalar transfer coefficients (sensible heat flux, latent heat flux and trace gases). Generally, land–atmosphere Hg⁰ flux measurements using micrometeorological methods are scarce, and information on detection limits is even rarer. For gradient-based systems, a minimum resolvable Hg⁰ concentrations gradient (MRG) is determined by mounting the sampling lines at the same height for several days (same-air test) and computing the concentration differences between the lines that are used for flux calculations. The MRG threshold is usually defined as the average plus 1 standard deviation of the concentration difference obtained by the same-air test. Fluxes are considered significant when the Hg⁰ concentration difference is above the MRG. Exemplarily, Edwards et al. (2005) derived a flux gradient system-specific MRG of 0.01 ng m⁻³ and a flux detection limit of 1.5 ng m⁻² h⁻¹. To calculate the flux detection limit of the gradient sampling system, site characteristics and atmospheric conditions have to be considered (see Eq. 8 in Edwards et al., 2005). Fritsche et al. (2008a) derived a MRG of 0.02 ng m⁻³ for their setup. The minimum determinable gradient-based Hg⁰ flux was between 0.5 and 4.6 ng m⁻² h⁻¹ (Fritsche et al., 2008b). Converse et al. (2010) and Zhu et al. (2015b) reported a similar MRG for their gradient-based micrometeorological systems of 0.07 and 0.06 ng m⁻³, respectively. During Hg⁰ flux studies over agricultural land in China, 57 % and 62 % of the aerodynamic and modified Bowen ratio measurements were significant (Zhu et al., 2015b). For Hg⁰ flux REA systems, Zhu et al. (2015b) reported that the absolute precision in the updraft and downdraft Hg⁰ concentration difference was concentration-dependent at 0.069±0.022 [Hg⁰] (ng m⁻³), whereas Osterwalder et al. (2017) determined a detection limit of 0.05 and 0.04 ng m⁻³. Over wheat canopy, 55 % of the fluxes were significant (Zhu et al., 2015a), whereas 52 % of the fluxes were significant over a boreal peatland (Osterwalder et al., 2017). The share of significant Hg⁰ fluxes for the gradient-based, REA and Eddy Mercury methods is in a similar range of approximately 50 %. However, the same-air tests applied to determine the detection limit of gradient-based and REA fluxes is more appropriate to compare with our approach to determine the zero flux in the laboratory. With a median zero flux of 0.074 ng m⁻² h⁻¹, the share of significant fluxes measured with Eddy Mercury would increase to 99.7 %, which is not realistic for measurements outside the laboratory environment. Generally, the reported mean fluxes derived from the gradient-based, REA and Eddy Mercury methods should include data below the detection limit because otherwise the magnitudes of the average exchange rates would be overestimated (see Fritsche et al., 2008a; Osterwalder et al., 2016).

Figure 6 Flux detection limit empirical probability distribution of the magnitude of flux measurements under outdoor conditions (black line). The boxplot insert shows the range of the magnitude of measured fluxes during the zero-flux experiment in the laboratory without Hg⁰ sources under very low-turbulence conditions. The black line shows the theoretical detection limit based on the statistical significance (p < 0.05) of the correlation coefficient between the vertical wind speed and Hg⁰ fluctuations.

Figure 7 Hourly aggregated diel cycle of (a) Hg⁰ fluxes and (b) simultaneously recorded CO₂ fluxes (bold green line) and fluxes measured in 2017 (red bold dashed line) for the same period. Each hour of the day represents the quantiles obtained from a 3 h window centered on the respective hour of all technically valid observations. The bold lines represent median flux values. The interquartile range (IQR) is the range of the middle 50 % of the data. The 70 % and 80 % confidence intervals (CI) and the number of measurements per hour (n) are given. The median CO₂ flux in 2017 is displayed (red bold dashed line). The IQR (red vertical lines) and the 70 % CI (lighter vertical lines) are indicated.

3.3 Net ecosystem exchange of Hg⁰ over grassland

The median (interquartile range, IQR) Hg⁰ flux measured at the Chamau (CH-Cha) research site using the Eddy Mercury system was 2.5 (−0.6 to 7.4) ng m⁻² h⁻¹. The Hg⁰ flux revealed a distinct diel pattern with median (IQR) daytime and nighttime fluxes of 8.4 (1.9 to 15) and 1.0 (−0.9 to 3.3) ng m⁻² h⁻¹, respectively. The minimum hourly median Hg⁰ flux (0.5 ng m⁻² h⁻¹) was detected at 21:00 (Fig. 7a). Emission of Hg⁰ reached a maximum between 11:00 and 14:00 (hourly median of 10.8 ng m⁻² h⁻¹). The diel Hg⁰ variation corresponded to solar radiation with the highest mean level of irradiance at 13:00 (606 W m⁻²). The flux of CO₂ changed from net emission during the night to net uptake by vegetation with sunrise (Fig. 7b). At noon, CO₂ fluxes were 26 % lower than the most negative flux occurring between 10:00 and 11:00 (−0.1 mg C m⁻² s⁻¹). The absence of a midday maximum CO₂ uptake indicates a midday depression due to plant stress from exceptionally hot and dry conditions. The partial closure of their stomata during the warmest period of the day minimizes water loss through transpiration with the consequence of lower CO₂ uptake. Overall, the median CO₂ flux during our measurement campaign in 2018 was only 24 % compared with the same period in 2017 which exhibited average climatic conditions (red dashed line in Fig. 7b). The median CO₂ uptake in 2018 was 0.031 mg C m⁻² s⁻¹ compared with 0.127 mg C m⁻² s⁻¹ measured in 2017. We suggest that the increased stomatal resistance of vegetation during the campaign in response to high drought stress not only led to the abovementioned minimized uptake of CO₂ but also damped stomatal gas exchange in general, including the uptake of Hg⁰. Subsequently, soil emission was the dominating factor driving the NEE of Hg⁰ during summer in 2018.

The Hg⁰ flux measured at the CH-Cha site is comparable to Hg⁰ fluxes reported for other grassland sites worldwide (Zhu et al., 2016). A median Hg⁰ flux of 0.4 ng m⁻² h⁻¹ and a flux range between −18.7 and 41.5 ng m⁻² h⁻¹ (site-based average fluxes) was reported for nine studies (Poissant and Casimir, 1998; Schroeder et al., 2005; Ericksen et al., 2006; Obrist et al., 2006; Fu et al., 2008a, b; Fritsche et al. 2008a, b; Converse et al., 2010). Several studies reported net Hg⁰ emission during summer. Converse et al. (2010) reported net average Hg⁰ emission of 2.5 ng m⁻² h⁻¹ from a high-elevation wetland meadow in Virginia, USA. Zhang et al. (2001) measured a Hg⁰ flux of 7.6±1.7 ng m⁻² h⁻¹ from an open background site in Michigan, USA. The average Hg⁰ flux from a grassland in Québec, CA, was 2.95±2.15 ng m⁻² h⁻¹, and a correlation of the diel flux pattern with solar radiation was reported (Poissant and Casimir, 1998). Average net Hg⁰ emission of 1.1 ng m⁻² h⁻¹ was recorded from a pasture in Ontario, CA (Schroeder et al, 2005). The mean Hg⁰ flux from four grassland sites in the USA ranged from 0.3 to 2.5 ng m⁻² h⁻¹ between May 2003 and 2004 (Ericksen et al., 2006). Fu et al. (2008a) reported average Hg⁰ fluxes ranging from −1.7 to 13.4 ng m⁻² h⁻¹ from three grasslands in China in August 2006. The mechanism driving Hg⁰ emission from grasslands is not fully understood. Photoreduction has been reported to enhance Hg⁰ emission from soils and the foliage surface and from Hg within foliar tissue (Gustin et al., 2002; Moore and Carpi, 2005; Choi and Holsen, 2009; Yuan et al., 2019). Soil warming has been suggested to promote Hg⁰ emission (Poissant et al., 1999; Zhang et al., 2001; Gustin et al., 2002; Almeida et al., 2009), which is likely due to the increased decomposition of organic material (Fritsche et al., 2008c) and facilitated mass transfer of Hg⁰ through the topsoil to the atmosphere (Lin and Pehkonen, 1999). Zhang et al. (2001) reported a strong positive correlation of Hg⁰ fluxes with solar radiation and soil temperature. A solar shielding experiment resulted in a 65 % decrease in soil Hg⁰ emission, suggesting that photoreduction is a major factor but also that soil temperature cannot be neglected.

Few grassland studies have shown net Hg⁰ dry deposition. Fritsche et al. (2008a) reported an average Hg⁰ flux of −1.7 (modified Bowen ratio) and −4.3 ng m⁻² h⁻¹ (aerodynamic gradient) during the vegetation period over a sub-alpine grassland at Frübül in central Switzerland, 15 km southwest of our study site. More summertime Hg⁰ fluxes from three central European grasslands were measured on a campaign basis, and average grassland–atmosphere Hg⁰ fluxes ranged from −4.3 to 0.3 ng m⁻² h⁻¹. The highest variability of the fluxes was recorded for the Neustift site in Austria with a range from −76 to 37 ng m⁻² h⁻¹ (Fritsche et al., 2008b). A second full year Hg⁰ flux study was performed at an upland meadow in Maryland, USA (Castro and Moore, 2016). The hourly mean summertime Hg⁰ flux was −1.2 ng m⁻² h⁻¹ and ranged between −224 and 354 ng m⁻² h⁻¹.

We found that the southern source area of our grassland site has a 28 % higher Hg substrate concentration (mean $=\mathrm{59.4}\pm \mathrm{8.4}$ ng Hg g⁻¹) than the northern source area (mean $=\mathrm{46.4}\pm \mathrm{5.1}$ ng Hg g⁻¹) (Wilcoxon two sample t test, p < 0.05; Fig. 8a). The Eddy Mercury system was able to resolve a marginally significant greater daytime Hg⁰ flux (+44 %, p=0.0515; Fig. 8c) and insignificantly greater nighttime Hg⁰ flux (+68 %, p=0.296; Fig. 8b) from the southern source area that is enriched in Hg compared with the northern source area. The proportionality of Hg⁰ emission to soil Hg concentration has been shown across Hg-enriched soils (Eckley et al., 2015; Zhu et al., 2018; Osterwalder et al., 2019), but no significant correlation has been observed for low-Hg level background soils (Agnan et al., 2016). There are two possible explanations for the lack of a significant relationship between Hg⁰ flux and soil Hg concentration: (i) the analytical uncertainty of Hg⁰ flux measurements, or (ii) a masking of Hg⁰ emission by stomatal uptake of Hg⁰ at vegetated surfaces that is independent of the soil Hg concentration.

Figure 8 Boxplots displaying (a) the total topsoil Hg concentration (0–10 cm) in the northern and the southern parcels as well as the Hg⁰ flux over the respective parcels (b) during the night (21:00–05:00) and (c) during the day (10:00–17:00).

3.4 Suggestions to improve the Eddy Mercury system

Here we propose a number of adjustments that are expected to improve the Eddy Mercury system's performance in particular by (1) facilitating data transfer and processing, (2) increasing the measurement frequency and sample air flow through the RA-915 AM and (3) achieving more stable temperature conditions in the field. The length of data gaps that are mainly caused by system calibrations should be reduced to the point where discussions about gap-filling methods and detrending procedures can be considered obsolete.

Improve data transfer. The determination of the time lag between the wind speed measurement and the Hg⁰ concentration measurement introduced a considerable source of uncertainty, and cross-correlation peaks had to be visually verified (Sect. 2.5.4.). In the future, we aim for a real-time transfer of raw data to the serial port instead of data transfer via the USB port on the embedded Windows 7 system of the RA-915 AM. This will allow better synchronization between the Hg⁰ measurements and the ultrasonic anemometer (Sect. 2.4) and will significantly facilitate post-acquisition data treatment.

Increase measurement frequency. The pilot campaign was performed with a measurement frequency of 1 Hz. In the future, we wish to increase the measurement frequency up to between 10 and 20 Hz. Such an increase in measurement frequency is possible via software adaptations of the RA-915AM, will make the oversampling of the Hg⁰ signal performed here (Sect. 3.2) redundant and will result in better counting statistics.

Increase sample flow rate. During this pilot study we connected a more powerful pump to the RA-915 AM and managed to increase the flow rate from standard operation of 7–10 to 14.3 L min⁻¹, resulting in a residence time in the measurement cell that was 2 times lower than the original residence time. The lower residence time in the cell reduced the dampening of the signal (Sect. 3.2). However, this high flow led to a reduction in the cell pressure (approximately 700 mbar) that affected the detection limit for the Hg⁰ concentration measurements. In the future, we propose to further reduce the residence time of the air in the measurement cell by increasing the sample air flow by another 30 % to 20 L min⁻¹ using an external pump. To minimize the pressure drop, we propose to reduce the constrictions present in the RA-915 AM by increasing the internal diameter of the valves and the inlet tubing.

Improve the long-term stability of the instrument. The stability of RA-915 AM Hg⁰ concentration measurements is temperature dependent (Sect. 3.2). We encountered strong diurnal temperature fluctuations of the instrument during the pilot campaign. We took several measures during the campaign to increase the temperature stability (e.g., placing the pump outside the temperature-controlled analyzer box, insulation of the analyzer box and shading it from direct sunlight). To improve the temperature stability in the future, we suggest placing the RA-915AM in an instrument box that has a better insulation and more powerful temperature control or, ideally, placing it in a climate-controlled instrumental hut. For long-term deployments of the Eddy Mercury, the sampling hose can be extended to bridge the distance between the air inlet, which is located close to the sonic anemometer and the instrumental hut where the system is placed. In this case, however, it would be important to guarantee a turbulent flow in the tube (Reynolds number of > 3000–3500; Lenschow and Raupach, 1991; Leuning and King, 1992), an adequate refresh rate in the sampling cell and to ensure that the pressure drop in the sampling cell was within the requirements of the instruments (> 600 mbar; Vladimir Ryzhov from Lumex Ltd, personal communication, 2018).

4 Conclusions

This study demonstrates an application of the EC method for Hg⁰ flux measurements over a grassland site with low soil Hg concentrations (< 100 ng g⁻¹). The maximum flux detection limit derived from a zero-flux experiment in the laboratory was 0.22 ng m⁻² h⁻¹. The statistical estimate of the flux detection limit under real-world conditions was 5.9 (50 % cutoff) to 13.7 ng m⁻² h⁻¹ (95 % cutoff). The Eddy Mercury system overcomes the major uncertainties of other micrometeorological methods that have previously been used for Hg⁰ flux measurements associated with intermittent sampling at two different levels (aerodynamic gradient methods) and the stringent sampling and analytical requirements (relaxed eddy accumulation). The Eddy Mercury system will considerably facilitate ecosystem-scale Hg⁰ flux measurement because it features a fully automated operation, cutting down operation costs for technical maintenance by experienced staff, argon supply and consumables. Eddy Mercury has the potential to be established as a standard micrometeorological method for long-term Hg⁰ flux measurements over grasslands and other terrestrial ecosystems. Such standardization of measurements is strongly required to obtain comparable data and properly evaluate controlling factors on the NEE of Hg⁰ on larger spatial and temporal scales (Obrist et al., 2018). Ultimately, the Eddy Mercury system could complement air pollution and greenhouse gas measurements within the global network of micrometeorological tower sites (FluxNet) (Baldocchi et al., 2001). The Eddy Mercury system also comes at an opportune time to include NEE measurements of Hg⁰ in the joint WHO and UN Environment project to “develop a plan for global monitoring of human exposure to and environmental concentration of mercury”.