Modeling the global atmospheric transport and deposition of mercury to the Great Lakes

Modeling the global atmospheric transport and deposition of mercury to the Great Lakes Mark D. Cohen, Roland R. Draxler, Richard S. Artz , Pierrette Blanchard, Mae Sexauer Gustin, Young-Ji Han, Thomas M. Holsen, Daniel A. Jaffe, Paul Kelley, Hang Lei, Christopher P. Loughner, Winston T. Luke, Seth N. Lyman, David Niemi, Jozef M. Pacyna, Martin Pilote, Laurier Poissant, Dominique Ratte, Xinrong Ren, Frits Steenhuisen, Alexandra Steffen, Rob Tordon, Simon J. Wilson


Introduction
Mercury contamination is an ongoing concern in the Laurentian Great Lakes region, with public health and wildlife toxicology ramifications (Evers et al., 2011;Wiener et al., 2012). Fish consumption advisories due to high mercury concentrations are widespread throughout the region (USEPA, 2016). Atmospheric deposition is likely the largest contemporary mercury loading pathway to the Great Lakes (Drevnick et al., 2012;Jeremiason et al., 2009;Lepak et al., 2015;Mason and Sullivan, 1997;Rolfhus et al., 2003;Sullivan and Mason, 1998) and recently deposited mercury may be more bioavailable than legacy contamination (Lepak et al., 2015;Orihel et al., 2008). It is therefore important to quantify the amount of mercury being deposited to each of the lakes to interpret spatial and temporal trends in ecosystem contamination and to compare atmospheric deposition to other loading pathways. Further, in order to prioritize actions to ameliorate the problem, it is important to quantify the relative importance of different source types and source regions to each lake's mercury loading.
Atmospheric measurements, by themselves, cannot adequately provide this needed information. While there are robust methods for wet deposition, there are no widely accepted methods to quantify mercury dry deposition by measurements alone. Perhaps more importantly, as there are practical limits on the number of sites that can be employed, measurements cannot fully characterize the spatial variations in atmospheric deposition -and therefore, the actual deposition to a given lake or watershed -over a large region with numerous significant, irregularly-spaced mercury emissions sources. Finally, while measurement-based chemical mass balance (e.g., Keeler et al., 2006), back-trajectory (e.g., Han et al., 2005Han et al., , 2007, and isotopic (e.g., Lepak et al., 2015;Sherman et al., 2015) methods have been used to estimate source-attribution information, these methods have not been able to fully provide quantitative, detailed source-attribution estimates for atmospheric mercury loading to the Great Lakes.
Comprehensive atmospheric fate and transport models have been used to estimate the amount and source attribution for mercury deposition to the Great Lakes and/or surrounding regions and sites (e.g., Cohen et al., 2004Cohen et al., , 2007Grant et al., 2014;Lin et al., 2012;Seigneur et al., 2004;Zhang et al., 2012b). However, as significant uncertainties persist in atmospheric mercury modeling (e.g., Ariya et al., 2015;Bieser et al., 2014) ambient measurements are essential to evaluate model accuracy. Thus, while models or measurements by themselves are significantly limited in their abilities to provide needed mercury deposition information for the Great Lakes and other receptors, they offer the most powerful analytical potential when used together.
In this study we used a new, global version of the HYSPLIT-Hg atmospheric fate and transport model to estimate the amount and source-attribution of 2005 mercury deposition to the Great Lakes. Processes in the model were characterized using a novel atmospheric lifetime analysis. Model accuracy was evaluated through extensive comparison with mercury measurements in ambient air and precipitation. To examine the sensitivity of the results to uncertainties in key inputs and assumptions, simulations were carried out with 11 different model configurations. This study provides unusually detailed source-attribution results for mercury deposition, giving lake-by-lake estimates of contributions arising from country-specific direct anthropogenic emissions (e.g., from the United States, Canada, China, etc.) and from specific natural and re-emissions processes (e.g., oceanic, land/vegetation, biomass burning, etc.).

HYSPLIT-Hg Model
We used HYPSLIT-Hg, a version of the HYSPLIT model (Draxler and Taylor, 1982;Draxler andHess, 1997, 1998;Stein et al., 2015) with special features added to simulate atmospheric mercury. A Lagrangian version of HYSPLIT-Hg was used to estimate the atmospheric transport and deposition of mercury to the Great Lakes from anthropogenic sources in the U.S. and Canada (Cohen et al., 2004(Cohen et al., , 2007, and throughout a European domain in a model intercomparison experiment (Ryaboshapko et al., 2007a(Ryaboshapko et al., , 2007b. With the incorporation of global Eulerian capability (Draxler, 2007), HYSPLIT and HYSPLIT-Hg can be run in a Lagrangian, Eulerian, or combination mode (Stein et al., 2015). HYSPLIT-Hg simulations presented here were carried out with an Eulerian-only configuration using a grid with horizontal resolution 2.5 o x 2.5 o and 17 vertical levels up to a height of ∼30 km (10 hPa). Meteorological data to drive the model (e.g., wind speed and direction, precipitation, relative humidity, etc.) were based on the NCEP/NCAR Global Reanalysis dataset (NCEP-NCAR, 1948-present), converted to HYSPLIT format (NOAA-ARL, 2003-present), specified every 6 hours on the same grid as used for mercury fate/transport.
A model spin-up period of 5 years with meteorological data for 2000-2004 was used, with 2005 emissions. Increasing the spin-up period from 2 to 3 years, 3 to 4 years, and 4 to 5 years resulted in marginal increases in modeled concentrations and deposition by ∼3%, ∼1%, and ∼0.2%, respectively. Thus, while a longer spin-up period could have been used, the differences in estimated 2005 results would likely change less than 0.2%.
It is recognized that due to changing emissions over

Chemical transformations
Gas-phase Hg(0) is converted to Hg(II) and Hg(p) by reaction with O 3 , OH•, H 2 O 2 , HCl, and Cl 2 in the HYSPLIT-Hg model. The current version of HYSPLIT-Hg does not include the potential oxidation of Hg(0) by bromine (Br/BrO), in part because of uncertainty in estimating the concentrations of reactive bromine in the atmosphere Saiz-Lopez and von Glasow, 2012;Simpson et al., 2015) (discussed further in Text S1). As summarized by Ariya et al. (2015), the partitioning of Hg(0) oxidation products among Hg(II) and Hg(p) forms varies from 0%-100% among atmospheric Hg models. In the absence of quantitative experimental measurement information, it was assumed that 10% of the product of the gas-phase oxidation by O 3 , OH•, H 2 O 2 is Hg(p) and 90% is Hg(II), while 100% of product of the HCl and Cl 2 oxidation reactions is Hg(II). In the aqueous-phase, Hg(0) is oxidized to Hg(II) by reaction with O 3 , OH•, HOCl, and OCl -1 , while Hg(II) is reduced to Hg(0) by photolysis of Hg(OH) 2 and by transformation of HgSO 3 -1 . The rate and equilibrium parameters used in the simulations are summarized in Tables 1 and  S1. Additional details about the estimation of reactant concentrations and phase partitioning are provided in Text S2. The HYSPLIT-Hg model does not currently include a simulation of bromine-mediated arctic atmospheric chemistry (e.g., Steffen et al., 2008).
In preliminary simulations using nominal literature values of the oxidation and reduction reactions shown in Table 1, the net lifetime of Hg(0) (discussed in more detail in conjunction with Figure 1, below) was unrealistically short (∼0.3 years) compared to the expected lifetime of ∼0.7-1.3 years (Ariya et al., 2015;De Simone et al., 2014;Holmes et al., 2010a;Lamborg et al., 2002;Seigneur et al., 2006;Selin et al., 2007). As a result, model-predicted atmospheric concentrations of Hg(0) were unrealistically low compared to typical measured concentrations (e.g., ∼1.3-1.7 ng m -3 in the Northern Hemisphere; Driscoll et al., 2013). The most significant oxidation processes in the chemical mechanism used here are the gas-phase reactions with OH• and O 3 (e.g., see Figure 1 below). When the oxidation rate constants for these two reactions were reduced by a factor of 5, realistic atmospheric lifetimes and concentrations of Hg(0) were obtained. This scaling approach is similar to that used in the GEOS-Chem model, where the Hg(II) reduction rate is adjusted so that reasonable Hg(0) concentrations are obtained (e.g., Holmes et al., 2010a). Here, we adjusted Hg(0) oxidation rates so that the model produces reasonable Hg(0) concentrations, but left Hg(II) reduction rates unchanged. The rates of the OH• and O 3 gas-phase reactions are considered to be highly uncertain. It has been argued that the rates may be significantly less than experimentally determined or may not occur at all (Ariya et al., 2015;Calvert and Lindberg, 2005;Subir et al., 2011). Therefore, we believe that the 1/5 scaling of these reaction rates used in the base case is within the range of uncertainty in the reaction rates for these two reactions.

Dry and wet deposition
Dry deposition of the various mercury forms from lowest-layer cells to terrestrial surfaces was estimated using a parameterized resistance-based approach (Wesely, 1989;Wesely and Hicks, 1977). For water surfaces, the approach of Slinn and Slinn (1980) was utilized. Dry deposition of both gas-phase and particle/droplet phase mercury was estimated. As a first approximation, a constant atmospheric particulate surface area of Elementa: Science of the Anthropocene • 4: 000118 • doi: 10.12952/journal.elementa.000118 Global atmospheric transport and deposition of mercury to the Great Lakes 4 3.5E-06 cm 2 cm -3 was utilized equal to the "background + local sources" value estimated by Whitby (1978). A typical particle size distribution based on Whitby (1975) -as cited by Prospero et al. (1983) -was used in this modeling. The assumed distribution was divided into 14 particle size bins, whose mid-point particlesize diameter ranges from 0.001-20 microns. The details of the distribution are summarized in Figure S1. To estimate the mass and mass fraction in each bin, based on the assumed surface area distribution, it was assumed that the particles were spherical with a density of 2 g cm -3 . With these assumptions, the mass loading of the entire distribution corresponds to 39 µg m -3 . Approximately 90% of the mass in the assumed distribution has a diameter of less than 10 µm; thus, the PM-10 (total particle mass < 10 µm) concentration associated with the assumed distribution is on the order of 35 µg m -3 .
When no liquid water was present -i.e., the particles were dry -Hg(0) and Hg(II) were assumed to be entirely in the gas phase, while Hg(p) and Hg2s were assumed to be entirely in the particle phase. In this case, Hg(p) and Hg2s were apportioned to the different particle size bins based on the fraction of the total surface area in each size bin. When liquid water was present, the condensed-phase concentrations of Hg(0), Hg(II), and Hg2s were estimated via thermodynamic calculations as described above. For Hg(0) and Hg(II), the total droplet phase mass was apportioned among the different size ranges based on the estimated volume fraction in each size range. With or without the presence of liquid water, Hg(p) and Hg2s were apportioned among the different size ranges based on the fraction of the total aerosol surface area in each size range.
Wet deposition was estimated based on the vertical location of a given cell relative to the cloud layer during precipitation events. If the cell was above the cloud layer, no wet deposition occurred. If the cell was within the cloud layer, the particle-phase pollutants Hg(p) and Hg2s were wet deposited at a rate governed by an estimated volume-based scavenging ratio (grams Hg per m 3 of precipitation / grams Hg per m 3 of air). As summarized by Gatz (1976) and Slinn et al. (1978), scavenging ratios for particle-phase pollutants associated with relatively small particle sizes -like Hg(p) -are relatively small, with typical values (in these units) less than 100,000. A scavenging ratio of 60,000 was used in these simulations, similar to the 40,000 ratio used for particle-phase pollutants in earlier, related HYSPLIT modeling (Cohen et al., 2002(Cohen et al., , 2004 (Pal and Ariya, 2004) b (Hall, 1995) c As discussed in the text, the Hg(0) oxidation rate is scaled to 20% of the nominal, literature value shown here. In other configurations, the rate is scaled to 33% and 10% of this value. d 10% of the product of this reaction assumed to be Hg(p) and 90% Hg(II). e (Tokos et al., 1998) (Seigneur et al., 1994) g (Calhoun and Prestbo, 2001), as cited by Bullock and Brehme (2002) h (Lin and Pehkonen, 1997) i (Munthe, 1992) j (Lin and Pehkonen, 1998) k Hg2s is Hg(II) adsorbed to soot, as described in the text. Equilibrium ratio shown coupled with 1 st -order time constant (60 minutes) for rate of approach to equilibrium. Follows the approach of Bullock and Brehme (2002), based on experimental results from Seigneur et al. (1998) l (Van Loon et al., 2000) (temperature T in degrees K) m (Bullock and Brehme, 2002;Xiao et al., 1994) Rate shown is maximum at peak insolation. In the simulation, the actual rate at any given location is scaled according to the ratio of local insolation to peak insolation. doi: 10.12952/journal.elementa.000118.t001 wet deposition of Hg(0) and Hg(II) was estimated using the precipitation rate and the thermodynamically estimated aqueous-phase concentrations.
For mercury in cells below a precipitating cloud layer, different approaches were used depending on the mercury form. Gas-phase Hg(II) and Hg(0) were scavenged assuming thermodynamic partitioning between the gas-phase and falling precipitation. Hg(p) and Hg2s associated with dry aerosol particles were scavenged using a size-dependent scavenging coefficient estimated for falling drops in the range of 0.04-0.4 mm (Seinfeld and Pandis, 1986). For each size range, the geometric mean value of the scavenging coefficient estimated for collectors of 0.04 mm and 0.4 mm was used. Hg(p), Hg2s, Hg(II), and Hg(0) associated with deliquesced aerosol particles below the cloud layer were scavenged using the same size-dependent approach used for dry particles, but the wetted (rather than dry) particle size was used to estimate the scavenging coefficient.

Model-estimated lifetimes for specific processes
A special set of simulations was carried out to investigate specific processes and process combinations. In these simulations, the entire three-dimensional domain was initially filled with a constant mixing ratio (1 ng / kg of air) of a particular form of mercury, i.e., Hg(0), Hg(II), or Hg(p). During the ensuing simulation, no additional mercury was added to the system. In any given simulation, one or more chemical transformations and/or deposition processes was allowed to occur, while others were disabled. Each specific case was simulated using 12 four-week ("monthly") simulations, with a different simulation started at the beginning of each month of 2005. The change in mass of the initial mercury form was tracked throughout.
Since there are complex, non-linear spatial and temporal variations in mercury's atmospheric fate and transport processes, different methodologies were used to estimate the 1/e lifetimes from the simulation results. The results are summarized in Figure 1 for Hg(0) and Figure 2 for Hg(II) and Hg(p). In one methodform ing the basis of the box-plot components of Figures 1 and 2 -lifetimes were estimated from the hourby-hour decreases in mass over each simulation. Assuming a quasi-first-order removal process, governed by dm/dt = -k m, hourly values of the removal rate constant k were estimated from the change in mass (dm/dt) and the mass in the system (m). The "instantaneous" (i.e., hourly) 1/e atmospheric lifetime τ is equivalent to the hourly value of k -1 . For each of the 12 monthly simulations we estimated the median of the instantaneous hourly-τ values. The box plot elements in Figures 1 and 2 show the statistical distribution (minimum, 25 th percentile, median, 75 th percentile, and maximum) of the 12 monthly median hourly-τ estimates. Also shown are lifetimes estimated from the hourly-τ values based on (a) the numerical average of the 12 monthly-medians, (b) the numerical average of all hourly-τ values in the first week of the 12 monthly simulations, and (c) massweighted average of all of the hourly-calculated instantaneous hourly-τ values over the 12 monthly simulations.  Values shown are estimates of the 1/e lifetime based on the instantaneous hour-by-hour decreases in the Hg(0) mass in the system, distributed initially throughout the entire threedimensional model domain at a constant mixing ratio. For simulations examining wet and/or dry deposition alone, all chemical transformation processes were turned off. For estimates of specific chemical transformation processes, wet and dry deposition processes were disabled, unless explicitly stated (i.e., the 1 st group of lifetimes shown is for net oxidation/reduction and deposition). For all estimates including OH• and O 3 oxidation reactions, values for 100%, 33%, 20%, and 10% reaction rate scaling are shown. For each process or group of processes, 12 four-week simulations were performed, starting at the beginning of each month of 2005. The box-plot elements show the statistical distribution of the medians of the individual hourly-calculated lifetimes for each monthly simulation, i.e., the 25 th , 50 th , and 75 th percentile, are shown as rectangles and the minimum and maximum are shown as whiskers of the 12 monthly medians are shown for each process. Also shown are the following statistical summary values: (a) the numerical average of the above monthly medians; (b) the overall average of the hourly-calculated estimates over just the 1 st week of the 4-week simulations; and (c) a mass-weighted average over all of the hourly-calculated estimates for each process.
As can be seen from the figures, the different methodologies sometimes show a range of approximate 1/e lifetime estimates for a given process. The ranges arise from numerous factors, including the fact that the processes are not spatially or temporally uniform. Further, the mass distribution of mercury is changedaffecting the efficiency of removal processes -as the simulation proceeds. This is particularly important for the estimation of deposition processes. For example, only material in the surface layer is removed during dry deposition, and so, the surface layer mass is depleted at a relatively fast rate at the beginning of the simulation. Once the initial surface layer material is depleted, the rate of dry deposition depends more on the rate of dispersion from layers aloft to the surface layer. A similar situation occurs with wet deposition for material at or below typical cloud heights. Once precipitation has removed this material, further wet deposition will depend more on the rate of dispersion from above layers. Thus, the 1 st -week average values may be the most relevant for wet and dry deposition processes. However, particularly for wet deposition, the variability of precipitation means that any given week will not necessarily be representative of the long-term average. Based on the above discussion, the atmospheric lifetimes presented in Figures 1 and 2 should be regarded as rough estimates.
In Figure 1 it can be seen that if the OH• and O 3 rates are specified at 100% of their potential values, the net modeled lifetime for Hg(0) due to oxidation and reduction would be on the order of 0.3 year (3-4 months). As noted above, because this is much shorter than the expected lifetime for Hg(0), we specified the rates of these reactions to be 20% of the potential value in our base case and considered variations of 10% and 33% in a sensitivity analysis. The net oxidation lifetime of Hg(0) considering all oxidation/reduction reactions considered in the model is ∼1 year in the base case, with the OH• and O 3 reactions scaled to 20%. For the comparable 10% and 33% scaling variations, the net Hg(0) lifetime against all oxidation/reduction in the model is ∼1.3 and ∼0.7 years, respectively. OH• oxidation shows more relative seasonal variation than other processes, i.e., the variation in the monthly median values are larger than other processes. The lifetime of Hg(0) with respect to dry deposition -with all chemical transformation processes turned off -is ∼3 years. We note that this lifetime is for the case where Hg(0) is distributed at a constant mixing ratio throughout the three-dimensional model domain. The average lifetime with respect to dry deposition for Hg(0) in the lower levels of the atmosphere would be somewhat less. The comparable lifetime of Hg(0) with respect to wet deposition is extremely long (∼30 years) owing to the minimal solubility of Hg(0) in water. The estimated lifetimes for this process -and other extremely slow processes explicitly simulated for Hg(0) -are less reliable, as the hourly changes in mass were too small to accurately quantify due to limits in numerical precision. While not shown here, numerical experiments with longer simulations suggested that the ∼30-year lifetimes shown for these very slow processes likely represent lower-bound estimates of their true values. The relatively long deposition-related lifetimes for Hg(0) suggest that deposition is a less important fate process in the model, as compared to chemical transformations, and, that the model will be less sensitive to uncertainties in the simulation of Hg(0) deposition. The relative unimportance of Hg(0) deposition processes can also be seen in Figure 1 by comparing the 1 st group of lifetimes (for oxidation/reduction and deposition) with the 2 nd group of lifetimes (for oxidation/reduction without deposition). When deposition processes are added in, the lifetimes for Hg(0) are only decreased on the order of ∼10-20%.

Figure 2
Lifetimes of Hg(II) and Hg(p) against selected chemical and physical processes in the HYSPLIT-Hg model.
Values shown are estimates of the 1/e lifetime of Hg(II) (left) and Hg(p) (right), based on the instantaneous hour-by-hour decreases in the Hg(II) or Hg(p) mass in the system, distributed initially throughout the entire three-dimensional model domain at a constant mixing ratio. As with the estimates for Hg(0): (a) all chemical transformation processes were turned off for simulations examining wet and dry deposition processes; and (b) all wet and dry deposition processes were disabled during simulations investigating specific chemical transformation processes. Interpretation of the box plot and symbol elements are as described in Figure  The lifetime of Hg(II) with respect to reduction in the model is ∼0.3 year and shows considerable seasonal variation (Figure 2). Wet and dry deposition processes are effective removal processes for Hg(II) with lifetimes of ∼0.1-0.2 year. The wet deposition of Hg(II) is not influenced by the value of the in-cloud particle scavenging ratio (WETR). For Hg(p), wet deposition lifetimes are ∼0.1 year and as expected, are somewhat dependent on the scavenging ratio. Hg(p) dry deposition is slower, characterized by a lifetime of ∼0.7 year.
These process-specific simulations demonstrate the relative importance of different chemical and physical phenomena to the simulation. For any given mercury form, the fastest processes (i.e., those with the shortest lifetimes) will have the greatest impacts on its atmospheric fate and transport. Accordingly, uncertainties in the fastest processes will have the greatest impacts on the accuracy of the simulation. Therefore, these results may be helpful in prioritizing efforts to reduce modeling uncertainties.

Mercury emissions
Mercury emissions used as model input included the following components: anthropogenic, biomass burning, geogenic, soil/vegetation, ocean, and prompt reemissions. The anthropogenic component was subdivided into emissions of Hg(0), Hg(II), and Hg(p), as described below. Emissions from all other components were considered as Hg(0). All inventory components were ultimately assembled on a global 2.5 o x 2.5 o grid, equivalent to the horizontal spacing of the global meteorological data used for the modeling.
Point-source anthropogenic mercury emissions for the U.S. were assembled from the USEPA 2005 National Emissions Inventory (NEI) (USEPA, 2009). For relatively minor, small, and widespread "area" sources (e.g., mobile sources) specified at the county level in the U.S., the 2002 NEI was utilized (USEPA, 2007), as it formed the predominant basis for the 2005 NEI area-source inventory. For point and area sources whose emissions were not separated into Hg(0), Hg(II), and Hg(p) in the NEI, EPA-recommended process-based "speciation" factors were utilized to estimate the emissions partitioning (USEPA, 2006). For point-source mercury emissions in Canada, Environment Canada's 2005 National Pollutant Release Inventory (NPRI) was utilized (Environment Canada, 2010). For Canadian area sources, 2000 data from Environment Canada were utilized, defined on a 100-km grid (Smith, 2008), as this was the latest data available for this analysis. While the use of 2000 (rather than 2005) data is a limitation, it is unlikely that changes in Canadian area source mercury emissions between 2000 and 2005 were large enough to significantly affect the results of this analysis. For point-and area-source mercury emissions in Mexico, the latest detailed inventory that was available was a 1999 inventory prepared for the Commission for Environmental Cooperation (CEC) (Acosta-Ruiz and Powers, 2001). The area-source emissions were geographically apportioned to each of the 32 Mexican states based on year-2000 population. Since mercury emissions in the Canadian and Mexican inventories were not separated into different mercury forms, the EPA-recommended speciation factors noted above were utilized to estimate emissions partitioning. For anthropogenic mercury emissions in the remainder of the world, the 2005 Arctic Monitoring and Assessment Program (AMAP) global inventory of Pacyna and colleagues was used (Pacyna et al., 2010;Wilson et al., 2010). The AMAP inventory includes global emissions of 880 Mg yr -1 from fossil-fuel combustion, 350 Mg yr -1 from artisanal and small-scale gold mining, ∼300 Mg yr -1 from other metallurgical processes, and lesser amounts from other source types. It is specified on a 0.5 x 0.5 degree grid (approximately 50 km x 50 km), with total emissions of Hg(0), Hg(II), and Hg(p) for each grid cell. Summing the North American inventories and non-North-American portions of the AMAP global inventory, global direct anthropogenic emissions totaled 1327, 475, and 125 Mg yr -1 , respectively, for Hg(0), Hg(II) and Hg(p). Temporal variations were not available in the above data sources and so anthropogenic emissions were assumed constant throughout the year.
Global mercury emissions from biomass burning were assumed to be 600 Mg yr -1 as utilized by Lei et al. (2013Lei et al. ( , 2014, consistent with the estimate by Friedli and colleagues (2009). This is slightly higher than the 200-400 Mg yr -1 values used in some other recent modeling analyses (e.g., Chen et al., 2014;Holmes et al., 2010b;Kikuchi et al., 2013;Song et al., 2015), but the difference is small compared to the total global emissions of mercury (on the order of 6000-7000 Mg yr -1 , as discussed below). Global mercury emissions from geogenic processes were assumed to be 500 Mg yr -1 , as used by Lei et al. (2013Lei et al. ( , 2014, Kikuchi et al. (2013), Holmes et al. (2010a), and Selin et al. (2008).
For soil/vegetation emissions and similar processes, the surface exchange of Hg(0) is bidirectional. In the HYSPLIT-Hg model simulations, emissions are specified as the gross, or one-way, upward flux, as opposed to the net, or bidirectional, upward flux. The downward component of the surface exchange is estimated as the simulation proceeds via run-time deposition modeling. Global, annual one-way mercury emissions from soil/vegetation were taken to be 1100 Mg yr -1 , similar to many other studies: e.g., 1100 Mg yr -1 was used by Selin et al. (2008), 890 Mg yr -1 was used in the base simulation of Kikuchi et al. (2013), and an optimized flux of 860 Mg yr -1 was recently estimated by Song et al. (2015). Prompt re-emissions of deposited Hg(II) after reduction to Hg(0) were assumed to be 30% of the total Hg(II) deposition to terrestrial surfaces, similar to the fraction used in other modeling studies (e.g., Jung et al., 2009;Selin et al., 2008). In this analysis, prompt re-emissions amounted to ∼400 Mg yr -1 , consistent with the 260-600 Mg yr -1 range used in other modeling studies (e.g., Holmes et al., 2010b;Kikuchi et al., 2013;Selin et al., 2008;Song et al., 2015). Taken together, the one-way Hg(0) emissions from soil/vegetation and prompt re-emissions totaled ∼1500 Mg yr -1 . As described below in the results section, gross, one-way dry deposition flux of Hg(0) to land surfaces was modeled to be ∼700 Mg yr -1 ; thus, the net Hg(0) emissions from land surfaces in the model was ∼800 Mg yr -1 . This is very similar to the recent estimate of Agnan et al. (2016) of the global net flux of Hg(0) from terrestrial surfaces of ∼600 Mg yr -1 , with a range of -500 to 1650 Mg yr -1 . This estimate includes roughly ∼100 Mg yr -1 of emissions from naturally enriched substrate, and thus, cannot be directly compared to the net terrestrial emissions used here, as we combined those emissions with volcanoes and other geogenic processes in a separate sub-inventory.
Global, annual, gross (one-way) mercury emissions from the ocean were taken to be 4350 Mg yr -1 . As described below, the gross, one-way deposition flux of Hg(0) to the global ocean surface was estimated to be 1650 Mg yr -1 . Thus, the net Hg(0) emissions flux from the ocean surface in this modeling analysis was 2700 Mg yr -1 , the same as the bottom-up estimate of 2700 Mg yr -1 developed from flux measurements (Pirrone et al., 2010), and consistent with the range of 2000-3600 Mg yr -1 used in numerous other modeling analyses (Amos et al., 2012;Chen et al., 2014;Corbitt et al., 2011;Holmes et al., 2010a;Kikuchi et al., 2013;Selin et al., 2008;Song et al., 2015).
Spatial and temporal (monthly) variations for the biomass-burning, geogenic processes, soil/vegetation, ocean, and prompt-reemission inventory components were adapted from the results of the Lei et al. (2014) analysis. Total emissions used in this analysis, using the net exchange of Hg(0) from surfaces, was 6500 Mg yr -1 . In Figure 3, the base-case emissions utilized in this study are compared with those used in other analyses (Bergan and Rodhe, 2001;Chen et al., 2014;Corbitt et al., 2011;Holmes et al., 2010b;Kikuchi et al., 2013;Lei et al., 2013;Mason and Sheu, 2002;Pirrone et al., 2010;Selin et al., 2007Selin et al., , 2008Shia et al., 1999;Song et al., 2015). Variations in addition to the base-case emissions, shown in this figure, will be described below. An overall map of total, global mercury emissions used in this modeling, in the base case, is shown in Figure 4. In carrying out the simulations used in this analysis, the HYSPLIT-Hg model was run separately for each of the sub-inventories (i.e., anthropogenic, biomass burning, ocean, etc.) in order to estimate the contribution of each emissions component to the overall deposition.  For each study, the net ocean and net terrestrial emissions are combined with the anthropogenic emissions to show the total emissions used in the analysis. The "net terrestrial" values shown represent the sum of net Hg(0) emissions from soil/vegetation, biomass burning, geogenic processes, and prompt reemission. In the Kikuchi et al. study (2013), several variations were presented in addition to the base case: M1 (with a new soil-emissions parameterization); M2-1 (with O 3 as an atmospheric oxidant); M2-2 (same as M2-1 but with a different treatment of polar emissions). The base case (1A) and other variations shown for this work are described in more detail in the text and involve variations in the rates of certain Hg(0) oxidation reactions, the extent of prompt Hg(II) reduction in plumes, and emissions. doi: 10.12952/journal.elementa.000118.f003

Model evaluation
Atmospheric mercury observations for 2005 were utilized to evaluate the model results, with particular emphasis on measurement sites in the Great Lakes region. The comparisons themselves are shown in the Results section below, but methodological aspects are briefly described here.
There are limited ambient mercury concentration monitoring data available for 2005, but we were able to obtain 2005 data for Hg(0), Hg(II), and/or Hg(p) at the sites summarized in Table 2 and Figure 5. For several of the sites, relatively complete data were available for the entire year for Hg(0) or Total Gaseous Mercury (TGM). For seven sites, only Hg(0) or TGM measurements were available during 2005. For the seven sites with Hg(II) and/or Hg(p) measurements, the fraction of the year covered by measurements was relatively low -from as low as 4% to as high as 37%. Each of the datasets had a particular sample averaging period, i.e., 1, 2, 3, or 24 hours, as shown in Table 2. In order to compare the model predictions against any particular dataset, the model results were assembled with a matching averaging period, in order to make the "fairest" comparison. For example, hourly-average measurement data were compared against hourly-average model results, while 24-hour average measurements were compared against 24-hour average modeling data. Comparisons were assembled only for the times that the measurements were made during the year at any given site. Further, as measurements are routinely reported at a standard temperature and pressure (STP) of 0 o C and 1 atm, the modeling results were converted to the same STP so that the comparisons could be made. In assembling the comparisons for the TGM measurements, model output of Hg(0) was used.
As can be seen from Table 2, seven of the sites had continuous, hourly TGM measurement data for essentially the entire year. For these sites, model vs. (versus) measurement comparisons were made for daily averages as well as for the hourly data. The median and mean concentrations of the hourly and daily data for any given site are essentially the same, but the hourly results show much greater variability. Due to the inherent inability of the simulation to capture sub-grid phenomena, these coarse-grid Eulerian model results are expected to have difficulty matching the variability of the hourly concentration values, but might have some skill in representing that of the daily-average values.
We compared Hg(II) model output concentrations with "Reactive Gaseous Mercury" (RGM) or "Gaseous Oxidized Mercury" (GOM) measurements. We did not include the model-output Hg2s concentrations in these Hg(II) comparisons, as it is unclear where in the measurement system Hg2s would register. For Hg(p) comparisons, we used the model output Hg(p) -including material on all particle sizes -to compare against the measured Hg(p) values. The measurements of Hg(p) are usually for relatively small particles only, i.e., particles ∼2.5 microns and smaller. We did not include the modeled Hg2s concentrations in the Hg(p) comparisons. For a few sites, both Hg(II) and Hg(p) were measured simultaneously (e.g., Mt. Bachelor and Reno-DRI). In these cases, we compared the total non-elemental mercury measured [i.e., the sum of Hg(II) and Hg(p)] against the total non-elemental mercury modeled. In this "non-Hg(0)" comparison, we included modeled concentrations of Hg2s, along with Hg(II) and Hg(p). A few of the sites are located in relatively complex terrain, and as discussed in Text S3, there is uncertainty as to which output model level should be compared to the measurements. In these cases, we have shown individual comparisons for all potentially relevant vertical model levels.
Measured mercury wet deposition fluxes and ambient concentrations were compared with HYSPLIT-Hg model output at the Mercury Deposition Network (MDN) (National Atmospheric Deposition Program, 2012) sites also shown in Figure 5. Since only total annual wet deposition results were evaluated, only the 86 MDN sites that operated for the entire year 2005 were considered. Of these 86 sites, 32 were in the Great Lakes region, 17 were in the Gulf of Mexico region, and 37 were elsewhere in North America. Precipitation measurement data at the 86 sites are compared in Figure S2 with the gridded precipitation data from the NCEP/NCAR meteorological reanalysis used to drive the HYSPLIT-Hg model. Due primarily to the coarseness of the model output grid (2.5 o x 2.5 o ), but also due to model and measurement uncertainties, an exact match is not expected. In comparing HYSPLIT-Hg model output with measured mercury wet  Base-case atmospheric mercury emissions from all source categories.
Annual total emissions of all forms of mercury (Hg(0), Hg(II), and Hg(p)) on the 2.5 o x 2.5 o global grid used in this modeling. Emissions shown in this map are gross "oneway" emissions used as input to the HYSPLIT-Hg model, as opposed to net emissions, and include contributions from anthropogenic, biomass burning, soil/vegetation, re-emissions, oceanic, and geogenic sources.  The locations of these sites are shown in Figure 5, labeled with the 3-letter abbreviations given in parentheses.  (Cole et al., 2013(Cole et al., , 2014Temme et al., 2007) e (Han et al., 2004(Han et al., , 2005(Han et al., , 2007 f (Zhang et al., 2012a) g (Swartzendruber et al., 2006) h (Peterson et al., 2009) i (Lyman and Gustin, 2008) doi: 10.12952/journal.elementa.000118.t002

Figure 5
Measurement sites with 2005 data used for model evaluation.
Sites with atmospheric mercury concentration data used in this analysis shown as larger white circles with 3-character site abbreviation, described in 2. deposition, the modeled flux at measurement locations was multiplied by the ratio of measured to modeled precipitation, in a post-processing procedure. It is recognized that sub-grid variations in precipitation will introduce complex, non-linear deviations in the simulations and that using the measured/modeled precipitation ratio at any given site is an approximation. In light of the above, and because there are other subgrid phenomena not captured in the simulation -e.g., the near-field impact of local sources -we would not expect an exact agreement between modeled and measured mercury wet deposition. We acknowledge that it can be informative to compare model results against wet deposition measurements with higher temporal resolution (e.g., as carried out by Bieser et al., 2014). However, given the coarseness of the meteorological data and Eulerian grid used in the analysis, comparison of the model results against weekly MDN data was not considered a realistic approach here.

Additional model configurations investigated
In addition to the base-case atmospheric chemistry and emissions configuration described above (configuration "1A"), several additional simulations with different configurations were carried out to investigate the impact of differing assumptions on model results (Table 3). In all of the "1" variations (1A-1C), base-case oxidation rate parameters were used, i.e., the gas-phase Hg(0) oxidation reactions by OH• and O 3 were scaled by 20%. In all of the "2" variations (2A-2D), the gas-phase Hg(0) oxidation reactions by OH• and O 3 were scaled by 33%. In all of the "3" variations (3A-3D), the gas-phase Hg(0) oxidation reactions by OH• and O 3 were scaled by 10%. For all 33% and 10% oxidation configurations (except 2D and 3D as noted below), the net oceanic and land-based emissions of Hg (0) were adjusted as shown in Table 3 so that the model would produce realistic Hg (0) concentrations. The net surface exchange of Hg(0) from oceanic and land surfaces is relatively uncertain and can be used as a "tuning parameter" for models (e.g., Song et al., 2015) as has been done here. The net oceanic emissions used in these alternate configurations (∼1400 and ∼4200 Mg yr -1 in 3A-3C, and 2A-2C, respectively, compared to the base case ∼2700 Mg yr -1 ) are within the overall range of 800-5500 Mg yr -1 used in other modeling analyses . The net terrestrial emissions used in the alternate configurations (∼350 and ∼1300 Mg yr -1 in 3A-3C, and 2A-2C, respectively, compared to the base case ∼800 Mg yr -1 ) are well within the range of -500 to 1650 Mg yr -1 estimated from field measurements (Agnan et al., 2016). As noted above in the Mercury emissions section, this latter estimate includes ∼100 Mg yr -1 of emissions from naturally-enriched substrates, which are considered part of a separate category of emissions in this modeling. Even with this difference in basis, the range in model configuration values is still likely well within the empirically-estimated range.
In configurations 1B and 1C, it was assumed that 33% and 67%, respectively, of anthropogenic emissions of Hg(II) were reduced to Hg(0) immediately after release. The prompt reduction of emitted Hg(II) in plume has been hypothesized by some to be more consistent with observations (e.g., Kos et al., 2013;Zhang et al., 2012b). In configurations 2B and 2C, these same plume reduction assumptions were combined with the 33% oxidation scaling assumption of configuration 2A. Comparable plume reduction assumptions were made in configurations 3B and 3C for the 10% oxidation scaling configurations.
For two alternative oxidation scaling configurations -2D (33% scaling) and 3D (10% scaling) -the base case emissions were used. While these configurations are somewhat implausible -considering the relatively unrealistic Hg(0) concentrations they produce -they are included to illustrate the sensitivity of modeling results in two different ways. First, configuration 2D is equivalent to 2A except that it has lower net surface exchange of Hg(0). This could result from increased one-way downward dry deposition or decreased oneway upward emission/re-emissions from the surface. Configuration 3D and 3A can be compared in a similar (but opposite) manner. Second, configuration 2D is equivalent to 1A except for increased Hg(0) oxidation. Analogously, 3D is equivalent to 1A except for decreased Hg(0) oxidation. Table 3 also shows a few model evaluation metrics for each configuration, e.g., the average bias in the model-estimated Great Lakes region Hg(0) atmospheric concentration and mercury wet deposition. More detailed model evaluation results are provided below.

Model evaluation
Comparison of modeled vs. measured mercury wet deposition Given the relatively coarse computational grid used in this work (2.5 o x 2.5 o ), it is not expected, as noted above, that any grid-average model result will match the measurements at any given measurement location. As is the general case with studies such as this, sub-grid-scale phenomena such as the near-field impacts of "local" sources on a given site will not be captured by the coarse Eulerian computational grid. With these tempered expectations, measurements are compared with grid-averaged model results. Figure 6 shows, and Tables S2 and S3 include, a comparison of modeled (base case) vs. measured 2005 mercury wet deposition at the 86 MDN sites for which essentially complete data were available for 2005. The comparison is differentiated among sites in the Eastern Great Lakes region, the Western/Central Great Lakes region, the Gulf of Mexico region, and all other MDN sites. As seen in Table S3, overall modeled deposition for the base case is within 11% of measured at sites in the Great Lakes region. Over all 86 MDN sites, the modeled deposition is within 28%. The poorest performance is for sites in the Gulf of Mexico region, where the model underestimates wet deposition by an average of 58%. As this work focused on the Great Lakes region, the results for the Gulf of Mexico region were not analyzed in depth. Model underestimates of mercury wet deposition in the Gulf of Mexico region has been found in other modeling studies (e.g., Amos et al., 2012;Bullock et al., 2009;Holmes et al., 2010a;Lin et al., 2012;Selin et al., 2007;Zhang et al., 2012b), and may be due to an inadequate characterization of the complex phenomena associated with deep convective thunderstorm scavenging of mercury (Nair et al., 2013). Without the 17 Gulf of Mexico sites, the model underestimated wet deposition by an average of 14%. Figure S3 and Tables S2-S3 show comparable summaries of modeled vs. measured wet deposition values for each configuration. Figure 7 shows that configurations with the slowest oxidation (3A-3C, the 3 leftmost symbols in each comparison) tend to produce the lowest average wet deposition in any given region. The configurations with the fastest oxidation (2A-2C, the 3 right-most symbols) show the highest modeled wet deposition. Within each series, as plume reduction decreases (from 67% to 33% to 0%), wet deposition increases slightly (e.g., the first 3 symbols: 3C, 3B, and 3A). Cases 2D and 3D where only oxidation rates (but not emissions) were changed show results similar to the base-case simulations. The model tends to underpredict mercury wet deposition in the Western/Central Great Lakes region (e.g., bias of -27% to -35% in cases 1A-1C), but has a tendency to over-predict in the Eastern Great Lakes region (e.g., bias of +3% to +16% in cases 1A-1C). Over all configurations, the model shows a moderate under-prediction (range in bias -38% to +9%) for the Great Lakes region as a whole.   Tables S4, S5, S6, S7, and S8. Results for different vertical model output levels are provided for sites located in more complex terrain. The extreme cases 2D and 3D show the expected tendencies. Case 2D, using a faster oxidation rate for Hg(0), but base-case emissions, under-predicts Hg(0) for most sites (median bias -17%). Conversely, case 3D, with slower oxidation kinetics, overestimates Hg(0) for most sites (median bias +26%). Results for all other configurations are relatively similar, generally falling within 100 pg m -3 of one another. For the three complex-terrain sites (Underhill, VT; Mt. Bachelor, OR; and Reno, NV), the model results for the highest potential level appear closer to the measured values, suggesting that the site is impacted by relatively elevated air, i.e., at the higher end of the range of possible model output levels to use in the comparison.
For sites with relatively complete 2005 data coverage in the Great Lakes (St. Anicet, QU; Burnt Island, ON; Egbert, ON; Point Petre, ON) and nearby regions (Bratts Lake, SK; Kejimkujik, NS), model results (except for implausible cases 2D and 3D) are encouragingly close to the measured concentrations (bias generally within +/-10%). For sites outside of the Great Lakes region, and sites with more limited data coverage, the results are more mixed (bias from -34% to +22%). For the Stockton (NY), Potsdam (NY), Paradise (NV), and Gibbs Ranch (NV) sites, the coverage of data during 2005 was relatively limited, with less than 10-20% of the year represented by measurements. The relatively high variability in measurements at these sites -especially for Paradise and Gibbs Ranch -can be seen in Figure S5. In these situations, the influence of sub-grid plume impacts may have been relatively large. Over the course of a year, if measurements were made continuously, the influence of sub-grid plume impacts might potentially balance out to some extent. But when measurements are made on a more sporadic basis, the influence may be more significant.

Comparison of modeled vs. measured Hg(II) and Hg(p)
Similar comparisons of modeled Hg(II), Hg(p), and total non-Hg(0) concentrations with measurements are shown in Figures 9, S6, S7, and S8 and Tables S8, S9, S10, and S11. As with Hg(0), the different model configurations give relatively consistent results at any given site, and the small differences among configurations seem reasonable. At any given site, Hg(II) and total non-Hg(0) concentrations increase with as oxidation rates increase (3A-3C < 1A-1C < 2A-2C), similar to the wet deposition results. The concentrations of Hg(p) estimated by the model are not strongly affected by the oxidation rates since only a small amount of Hg(p) is assumed to be formed during atmospheric Hg(0) oxidation in these simulations.
Overall, it is seen that the model tends to estimate higher than reported concentrations of Hg(II), Hg(p), and total non-Hg(0), with median modeled concentrations on the order of 2-3 times the measured concentrations. This discrepancy has been found in other modeling studies. Zhang et al. (2012b) found that the GEOS-Chem model overestimated Hg(II) and Hg(p) by a factor of 4 and 2, respectively, unless it was assumed that a significant fraction (∼75%) of Hg(II) emitted by coal-fired power plants was immediately reduced to Hg(0) in the downwind plume. With the assumed plume reduction of emitted Hg(II), the model overestimated Hg(II) and Hg(p) by ∼40%. In another example, "reactive mercury", defined as the sum of Hg(II) and Hg(p), was overestimated by a factor of 2.5 relative to measurements (Weiss-Penzias et al., 2015). A summary of model vs. measured concentrations of Hg(II) and Hg(p) in the Great Lakes region showed that models generally overestimated measurements by a factor of 2 to 10 (Zhang et al., 2012a). In a related study, large model overestimates of Hg(II) and Hg(p) were also found (Kos et al., 2013). In a detailed model evaluation study, the CMAQ model overestimated measured concentrations of Hg(II) and Hg(p) by an average factor of ∼3.6 (Holloway et al., 2012).
In the simulations carried out here, the extent of assumed Hg(II) plume reduction did not dramatically influence model-estimated Hg(II) concentrations at model evaluation sites. Figure 9 and Table S8 show that the highest model-predicted Hg(II) concentrations at any site were estimated from configurations 2A, 2A, and 2C -with the fastest oxidation rate (33% scaling) -but there was relatively little difference among these "2-series" configurations where the assumed plume reduction varied from 0% (2A) to 33% (2B) to 67% (2C) (median model/measured ratios fell from 3.7 to 3.5). Similarly, there is little difference seen among the plume reduction variations within the three 20% oxidation scaling configurations (1A, 1B, 1C, with median ratios 2.8-2.6), or the three basic 10% scaling configurations (3A, 3B, 3C, with median ratios from 2.1-1.9).
T h e a ve r a g e m e a s u re d concentrations (M), shown as black circles, are shown along with the model estimates for the base case (1A) and other model configurations for each site with 2005 data included in the analysis. Data for Hg(II) (left), Hg(p) (center) and total non-Hg(0) (right) mercury forms are shown. The total non-Hg(0) model estimates include Hg2s, as well as Hg(II) and Hg(p). Hg(II) model results are compared against measurements generally reported as Gaseous Oxidized Mercury (GOM). Hg(p) model results -representing the total particle-associated mercury over all size ranges -are compared against fine-particle mercury measurements, which typically include only particles smaller than 2.5 um. Unless otherwise specified, model concentrations are reported for the lowest model layer (L2), i.e., the layer closest to the earth's surface (L1 is defined as the surface itself ). For sites in more complex terrain, model results for additional model elevations or "levels" (L2, L3, etc.) are shown. All concentrations are reported at 0 o C and 1 atm, and the comparisons are made only for the specific times during 2005 that measurements were made. doi: 10.12952/journal.elementa.000118.f009 It is interesting to consider these model vs. measurement concentration discrepancies for Hg(II), Hg(p) and total non-Hg(0) together with the wet deposition discrepancies discussed earlier. If the model underestimated wet deposition of these mercury forms as a result of inaccurate parameterizations or unrealistic assumptions, then this could have contributed to the overestimated atmospheric concentrations found for these same mercury forms. Given that wet deposition was found to be very fast for these mercury forms in the pathway-specific simulations discussed earlier (Figure 2), uncertainties in precipitation and the model's wet deposition estimation methodologies could have a significant effect. Figure 2 shows that the scavenging ratio (one of the deposition parameters in the model) influences the deposition rate of Hg(p), and this is just one of many potential uncertainties. Further, since the wet deposition of Hg(p) may be faster than that of Hg(II) (Figure 2), errors in the assumed distribution of products from the atmospheric oxidation of Hg(0) could be important. In this modeling, we assumed that 10% of the oxidation products were Hg(p) and 90% Hg(II) ( Table 1). If the true proportion of Hg(p) in the oxidation products was higher, this would likely lead to higher overall mercury wet deposition.
As the discussion above illustrates, some or even most of the discrepancies between modeled and measured concentrations may be due to uncertainties in emissions and the characterization of atmospheric fate and transport processes in the model. However, there are a few additional possible reasons for the tendency of the model-estimated concentrations to be greater than the reported measurements.
First, the measurement of Hg(p) is somewhat uncertain. In a comparison of Hg(p) measured with manual (filter) and automated (Tekran) methods, Talbot and colleagues (2011) found that manually-measured concentrations were ∼21% higher than automated concentrations, on average, and that as much as 85% or more of the data showed a difference of greater than 25%. Further, the measurements of Hg(p) are typically limited to Hg(p) on particles less than about 2.5 microns in diameter. The reported measurements are more accurately called "Fine Particulate Mercury", rather than "Total Particulate Mercury". In the modeling results, we are showing the estimates for total particulate mercury. There are very limited data on the size distribution of Hg(p). In a related study investigating aerosol mercury at a marine and a coastal site, significant mercury was often found in coarse aerosol size fractions (Feddersen et al., 2012). In some cases, the majority of Hg(p) appeared to exist on particles larger than 2.5 microns in size, although this was not always found. Keeler et al. (1995) found that, on average, about 88% of Hg(p) was associated with particles less than 2.5 microns in measurements in Detroit, but this varied from 60%-100%, depending on conditions. In the light of the above, one would generally expect the modeled total Hg(p) to be somewhat higher than the reported measurements of fine-particle Hg(p), even if both were "perfect".
For Hg(II), experimental evidence is emerging that suggests the commonly used, operationally-defined measurements of "Gaseous Oxidized Mercury" (GOM) [also sometimes called "Reactive Gaseous Mercury" (RGM)] using coated denuders may be underestimating the true concentration in the atmosphere. For example, Lyman et al. (2010) found that atmospheric ozone reduced the efficiency and retention of GOM capture on KCl-coated denuders. Denuders exposed to ozone lost from 29-55% of loaded HgCl 2 and HgBr 2 . Collection efficiency of denuders decreased by 12-30% for denuders exposed to 50 ppb of O 3 . Model vs. measurement comparisons are also complicated by an incomplete understanding of the actual chemical speciation of "GOM". Additional examples and discussion of the above issues are provided by Ariya et al. (2015), Gustin et al. (2015), Jaffe et al. (2014), Kos et al. (2013), andWeiss-Penzias et al. (2015).

Global mercury budget
The modeled 2005 emissions and deposition fluxes in the base case (1A) are summarized in Figure 10. The base case simulated a total global one-way Hg(0) deposition rate of ∼2400 Mg yr -1 . Excluding anthropogenic emissions, but including biomass burning and geogenic processes, the base case included total one-way Hg(0) surface emissions of ∼7000 Mg yr -1 . Thus, the "net", non-anthropogenic Hg(0) evasion flux at the earth's surface in the base case was estimated to be 4600 Mg yr -1 . This value is similar to that found in other analyses, e.g., 4800 Mg yr -1 (Selin et al., 2007), 4000 Mg yr -1 (Selin et al., 2008), 5200 Mg yr -1 (Pirrone et al., 2010), 3000 Mg yr -1 (Holmes et al., 2010a), and 3800 Mg yr -1 (Corbitt et al., 2011). These and other comparisons can be seen graphically in Figure 3.
In the base-case simulation, the total deposition was ∼3% less than the total emissions, and this was found with all model configurations (as shown in Table 3). The apparent tendency of the model to underestimate mercury wet deposition may have contributed to this imbalance. At the same time, one would not expect the actual atmospheric emissions and deposition to exactly balance out in any given year for any of the plausible model configurations (i.e., all except 2D and 3D), given the overwhelming array of stochastic processes contributing to the totals and the associated fact that emissions are continually changing in space and time.   Figure 11 shows the total model-estimated 2005 deposition to each of the Great Lakes and the fraction of this deposition arising from different source types. The values shown are for deposition directly to the lakes. The figure also shows estimates for an area-weighted average over all five of the Great Lakes. In each panel of the figure, estimates are shown for the base case as well as the 10 alternative model configurations considered. For all model configurations, the contributions to each of the Great Lakes due to the direct, anthropogenic emissions from the United States, Canada, Mexico, China, Russia, India, and all other countries were estimated and are also shown in Figure 11. In most cases, overall deposition fluxes decreased as the oxidation rate decreased (e.g., 2A > 1A > 3A) and deposition decreased as plume reduction increased (e.g., 1A > 1B > 1C).

Source-attribution for Great Lakes mercury deposition
As discussed at numerous points above, there are significant uncertainties in atmospheric fate/transport processes as well as in the spatial, temporal, and chemical distribution of emissions. While configurations 2D and 3D are considered implausible due to the somewhat unrealistic Hg(0) concentrations they produced, they provide some information about the model's sensitivity to these uncertainties. These configurations used base-case emissions but different Hg(0) oxidation rates. Their deposition amounts and source attribution estimates were similar to the base case (1A), suggesting that uncertainties in oxidation rates alone may not strongly influence the results. Also, configurations 2D and 3D are the same as configurations 2A and 3A, respectively, except for significant differences in the net surface exchange of Hg(0) at oceanic and terrestrial surfaces. The pair-wise comparison of these two sets of variations (2D vs 2A, and 3D vs. 3A) shows that the qualitative aspects of the results for any given lake are not dramatically affected by uncertainties in net surface exchange of Hg(0).
Direct anthropogenic emissions tended to have the highest impact on Lake Erie, with 58% of the total estimated deposition in the base case and a range of 38-64% among the different configurations [58% (36-64%)]. These emissions had the lowest impact on Lake Superior, contributing an estimated 27% (22-33%) of deposition, with intermediate impacts on the other lakes.
United States anthropogenic sources directly contributed on the order of 25% (12-29%) of 2005 deposition to the Great Lakes basin, on average, using different simulation methodologies. The importance of US contributions to individual lakes varied from 46% (24-51%) for Lake Erie to 12% (6-13%) for Lake Superior.
China was generally the country with the 2nd highest anthropogenic contribution, contributing on the order of 6.0% (5.1-9.7%) of 2005 deposition to the Great Lakes basin. The importance of Chinese anthropogenic emissions contributions to individual lakes was much more uniform than that for the U.S. (ranging from 4-11% over all lakes and configurations), as might be expected given their distance from the lakes.
Canada's direct anthropogenic sources were estimated to contribute on the order of 2.0% (1.0-2.3%) of the total model estimated deposition to the Great Lakes during 2005. India, Mexico, and Russia were all estimated to contribute on the order of ∼0.5-1% of the total 2005 deposition. The total direct anthropogenic contribution from all "other" countries in the world during 2005 was on the order of 4.3% (3.6-7.2%) of the total model-estimated deposition. The remainder of the deposition for the Great Lakes, on average, came from oceanic natural and re-emissions of previously deposited mercury [32% (27-42%)], terrestrial natural emissions and re-emissions [17% (14-21%)], biomass burning [5.1% (4.3-6.9%)], and geogenic emissions [(6.4% (5.6-8.4%)].
As expected, the extent of plume reduction of emitted Hg(II) affects the modeled source-attribution results. An illustration of this is provided in Figure 12 which shows the impact of U.S. and Chinese anthropogenic emissions on Lakes Erie and Superior as a function of the assumed plume reduction fraction. The relative contribution of U.S. emissions decreases with increasing plume reduction. Emitted Hg(II) from sources in the U.S. upwind of the Great Lakes can have a relatively large depositional impact, but if reduced in plumes to Hg(0), the impact will be lessened. The effect of plume reduction is larger for Lake Erie than for Lake Superior as there are much greater local/regional Hg(II) emissions upwind of Lake Erie. Conversely, the relative contribution of Chinese emissions increases with increasing plume reduction. In this case, if less emitted Hg(II) is deposited locally and regionally due to plume reduction, more mercury is available to be transported and ultimately deposited to the Great Lakes. The differences among the different configuration series shown in Figure 12 are due to increased oxidation of the global atmospheric Hg(0) pool and the resulting increased relative deposition to the Great Lakes. There is more such Hg(0) oxidation in the 33% oxidation scaling series (2A, 2B, 2C), and so, the relative contribution of U.S. anthropogenic sources is decreased. Conversely, there is less oxidation in the 10% oxidation scaling series (3A, 3B, 3C), and so the relative impact of U.S. emissions is increased. A very similar relative pattern among the different configuration series was found for Chinese anthropogenic emissions, but the overall impacts and absolute magnitudes of differences between series are much smaller and cannot be easily seen.

Figure 12
Fraction of total modeled deposition attributed to direct anthropogenic emissions from the USA and China.
Results are shown for Lake Erie and Lake Superior, to illustrate differences in the modeled source-attribution patterns for 2005. Each panel shows the results for a particular series of model configurations. The assumed oxidation scaling for each series (20%, 33%, or 10%) is shown, and within each series, the data are plotted as a function of the assumed amount of plume reduction of emitted Hg(II). doi: 10.12952/journal.elementa.000118.f012 These overall source-attribution results are similar to the findings of other source-attribution studies. As one example, Zhang et al. (2012b) found that 12-21% of 2008-2009 mercury deposition to the contiguous U.S. was attributed to North American anthropogenic sources, but 22-39% of deposition in the Ohio River valley region was attributed to North American anthropogenic sources. In another GEOS-Chem-based analysis (Selin et al., 2008), 20% of contiguous U.S. mercury deposition for 2004-2005 was attributed to North American anthropogenic sources, but as much as 50-60% of the deposition in the industrial Midwest and Northeast was attributed to these sources.
In another analysis, the GEOS-Chem model was used to develop source-attribution estimates for North America as a whole -i.e., all of Mexico, the United States, and Canada -for 2005-2011(Chen et al., 2014. Nine percent of North American deposition was attributed North American anthropogenic emissions, while 27% and 64% of North American deposition was attributed to anthropogenic emissions from other regions and "natural" emissions, respectively. Natural emissions in this analysis included emissions and re-emissions from all oceanic and terrestrial emissions pathways other than direct anthropogenic emissions. As discussed above, these "natural" pathways include re-emissions of mercury from previously deposited anthropogenic emissions. Simulations with the CAM-Chem/Hg model (Lei et al., 2013) found that 22% of total 1999-2001 mercury deposition in the United States could be attributed to direct U.S. anthropogenic emissions, while in industrial regions, ∼50% of the deposition was due to these domestic emissions.
The CMAQ model (with GEOS-Chem boundary conditions) was used to estimate source-attribution for 2005 mercury deposition to 6 regions within the U.S. (Lin et al., 2012). Results showed that 25% of deposition in the East Central and 29% of deposition in the Northeast U.S. could be attributed to U.S. direct anthropogenic emissions, while the overall average over the continental U.S. was 11%. CMAQ model simulations with boundary conditions generated by the Mozart global model (Grant et al., 2014) found that U.S. fossil-fueled power plants contributed significantly more to Lake Erie than to other Great Lakes, with contributions ranging up to ∼50% in portions of the lake during some seasons. Contributions to mercury deposition from Chinese anthropogenic emissions varied from ∼1-10% over the Great Lakes basin, with the largest relative impacts on Lake Superior.
In a measurement-based study, approximately 70% of the wet deposition of mercury in Eastern Ohio (in the vicinity of Lake Erie) could be attributed to local and regional sources (Keeler et al., 2006).

Discussion
The analysis has considered "direct" anthropogenic emissions, as well as emissions from other terrestrial and oceanic processes. It is important to note that past anthropogenic emissions (and subsequent deposition) have contributed to most of these other processes. For example, emissions from the ocean reflect "natural" levels of mercury that would have been present without anthropogenic influences, of course, but also reflect mercury in the ocean system that is present because of previous "direct" anthropogenic emissions. Similarly, emissions from soil/vegetation, prompt re-emissions, and biomass burning all include "indirect" anthropogenic contributions. So, the actual anthropogenic contributions to Great Lakes mercury deposition are larger than those reported above.
Model results presented here are estimates for 2005, and the amounts and speciation of mercury emissions have changed since that time. For example, overall emissions and emissions of Hg(II) from coal-fired power plants in the U.S. have been significantly decreased over the past decade (Castro and Sherwell, 2015;Zhang et al., 2016). It is expected that the fraction of Great Lakes deposition attributed to direct emissions from U.S. anthropogenic sources will be decreased due to these changes. The results presented here suggest that U.S. emissions played a significant role in contributing to atmospheric mercury deposition to the Great Lakes in 2005. If this analysis is valid, then reductions in U.S. emissions since that time could contribute to reductions in atmospheric mercury in the region. Recent analyses have suggested that reductions in U.S. emissions have indeed led to observable, decreasing trends in environmental mercury concentrations at many sites in the Great Lakes and surrounding regions (Castro and Sherwell, 2015;Giang and Selin, 2016;Hutcheson et al., 2014;Sunderland et al., 2016;Weiss-Penzias et al., 2016;Zhang et al., 2016), lending credibility to the results presented here.

Conclusions
An Eulerian version of the HYSPLIT-Hg model was used to simulate the transport and deposition of globally emitted mercury to the Great Lakes for 2005. A novel and unusually detailed process-specific atmospheric lifetime analysis was carried out, providing important insights into the fate/transport behavior of atmospheric mercury. Individual estimates have been provided for each of the lakes, allowing differences among the lakes to be examined. Simulations were carried out for a base-case and 10 alternative model configurations, allowing for a detailed exploration of the sensitivity of model results to uncertainties. Model results and sensitivities were evaluated by comparisons with 2005 ambient measurements of mercury wet deposition and atmospheric concentrations in the Great Lakes and other regions. In making these comparisons, it is recognized that the coarse Eulerian grid used in the analysis will not have captured potentially significant near-field impacts from sources. Aside from two configurations that were included as a sensitivity analysis and not regarded as being plausible, all of the configurations gave relatively similar results: (a) a tendency to underestimate mercury wet deposition, e.g., by ∼10% in the Great Lakes region, on average; (b) relatively good agreement with Hg(0) measurements, especially for sites with essentially continuous monitoring throughout the year; and (c) a tendency to produce non-Hg(0) concentrations greater than those measured, a discrepancy found in other modeling studies. The disagreement is surely due partly to model inaccuracies but may also in part be due to different atmospheric Hg fractions being compared (e.g., measured Fine Particulate Mercury (< 2.5 µm) compared with modeled total particulate mercury, and in part due to the measurements potentially being biased low.
Total mercury deposition to each of the Great Lakes was estimated for the base case and alternative configurations. Additionally, detailed source-attribution estimates have been provided, including configurationby-configuration, and lake-by-lake results for 6 specific countries and for 6 overall source types. Differences were found among the Great Lakes, and among the different configurations, in both total deposition and source-attribution for that deposition. Lake Erie, downwind of significant local and regional mercury emissions sources, was estimated by the model to be the most impacted by direct anthropogenic emissions (∼ half of the deposition), and Lake Superior, with fewer nearby emissions, was estimated to be the least anthropogenically impacted (∼ one quarter of the deposition). Among direct anthropogenic emissions contributions for the base case simulation, the U.S. was the largest contributor, followed by China, contributing 25% and 6%, respectively, on average, for the Great Lakes. The contribution of U.S. direct anthropogenic emissions to total mercury deposition varied between 46% for the base case (with a range of 24-51% over all model configurations) for Lake Erie and 11% (with a range of 6-13%) for Lake Superior. These results illustrate the importance of atmospheric chemistry, as well as emissions strength, speciation, and proximity, to the amount and source-attribution of mercury deposition. the provision of Canadian gaseous mercury data, accessed from the Government of Canada Open Data Portal (Environment and Climate Change Canada). We thank the HYSPLIT modeling group at the NOAA Air Resources Laboratory (ARL), including Glenn Rolph, Barbara Stunder, Ariel Stein, Fantine Ngan, Alice Crawford, and Tianfeng Chai for helpful discussions, the ongoing co-development of the HYSPLIT model, and conversion of the NCEP/ NCAR meteorological data to HYSPLIT format. We thank Rick Jiang and colleagues in the ARL IT department for essential assistance in establishing and maintaining the computational resources required in this project. Finally, we thank Johannes Bieser, Ian Hedgecock, and an anonymous reviewer for thoughtful and constructive review comments.

Funding information
MC received partial funding for this work from the U.S.EPA through the Great Lake Restoration Initiative.

Competing interests
The authors have declared that no competing interests exist.