Interactive comment on “ Quantitative assessment of upstream source influences on total gaseous mercury observations in Ontario , Canada ”

General comments: 1. In this paper, a modified STILT model is used to predict ambient TGM concentrations at three sites in North America. The result shows that the modified STILT model could better illustrate the near-field influences compared to CMAQ model. I noticed that the three sites, at which comparisons were made, are all remote sites. Are there some nearby sources (in the sub-grid scale) at one or two sites? I wonder why the authors chose these sites. If the authors want to show the model’s ability to account for influence of nearby sources, I think at least one site which is affected by nearby sources (e.g. industrial or urban site) in the sub-grid scale should be used here.


Introduction
Mercury was one of the first priority PBT (Persistent, Bioaccumulative and Toxic) pollutants identified by the US EPA, due to its significant health influences, especially damage to the nervous and reproductive systems.Human activities, such as mining and burning of fossil fuels, are important sources of atmospheric mercury, but there are also many "natural" sources, including soil, minerals and water.(Here, "natural" sources include mercury that is released into the soil from primary mineral sources and re-emission of mercury which was deposited from historical anthropogenic activities.) In the atmosphere, mercury can exist in three major chemical forms: elemental mercury (Hg 0 ), particulate mercury and "reactive" mercury (Hg ++ , also referred to as Reactive Gaseous Mercury, or RGM).Because of its high vapour pressure and low reactivity, Hg 0 can be transported through the atmosphere over very long distances before being returned to the surface by wet and dry deposition.Deposition from the atmosphere is the primary source of mercury contamination to most threatened aquatic ecosystems (Bullock, 2000).
Measurements of atmospheric mercury levels are important for monitoring health effects, but measurements alone are insufficient for a complete understanding of mercury sources, sinks and transport.To assess which source regions and types are responsible for observed mercury pollution, it is necessary to use numerical models (Cheng and Schroeder, 2000;Cohen et al., 2004;Lim et al., 2001;Lynam and Keeler, 2006).In this context, STILT, a receptororiented model, has proven to be especially useful for identifying transport pathways and estimating surface emission fluxes.STILT has been used to study terrestrial carbon fluxes at the regional scale using observations of CO 2 and CO over North America (Gerbig et al., 2003) and to estimate fluxes of halocarbons using gridded CO emissions and measured CO/halocarbon emission ratios (Hurst et al., 2006).STILT D Wen et al.: Quantitative assessment of upstream source influences has also been used to study transport in deep convective clouds (Xueref et al., 2004).
STILT simulates the transport of particles (representing air parcels of equal mass) with both deterministic and stochastic velocities, enabling a more detailed and accurate representation of air parcel trajectories, particularly in the lower atmosphere where turbulence is strong and where the traditional approach of simulating atmospheric transport with single mean-wind trajectories can be subject to large errors.The particles simulated by STILT comprise the observation's transport history, providing invaluable information for the interpretation of atmospheric observations.Since the particles are transported by wind vectors interpolated down to their sub-grid scale point locations they have the potential to resolve sub-grid scale influences, which are particularly important in cases where strong and variable sources/sinks are found in the near-field of a measurement site.Advection by Lagrangian particles is also known to minimize numerical diffusion that is often found in Eulerian advection schemes (Odman, 1997).
The objective of this study is to present a new method which can be used to simulate Hg concentrations at the location of a monitoring site and assess quantitatively its upstream source influences.For this purpose, we adapted STILT to simulate hourly total gaseous mercury (TGM) concentrations at three monitoring sites in Ontario, Canada.The simulated concentrations were evaluated with the measurements and also compared with the concentrations previously simulated by the Eulerian CTM CMAQ (Gbor et al., 2007).The differences in simulation ability between the two models, which were driven with the same meteorological fields and emission grids, were also investigated.

Total Gaseous Mercury measurements
Total Gaseous Mercury (TGM) is mainly composed of elemental mercury vapour with a minor fraction of RGM.TGM concentrations have been monitored continuously at the sites of the Canadian Atmospheric Mercury Measurement Network (CAMNet) since 1994, using the Tekran Model 2537A Ambient Mercury Vapour Analyser, which has a monitoring error of less than 2% (Poissant, 2000).Hourly TGM measurements are available from the Canadian National Atmospheric Chemistry (NAtChem) database (Environment Canada, 2002).We used these TGM measurements for the year 2002 at three monitoring sites in Ontario: Point Petre, Egbert and Burnt Island.Results from a previous CMAQ modeling study (Gbor et al., 2007) were also included as a comparison to the STILT results.The locations of the three sites are shown in Fig. 1.Burnt Island (45.8083   (44.2317 • N, 79.783 • W) are designated as "rural-affected" (Kellerhals et al., 2003).They are all situated between a large area of very low Hg emissions to the north and a region to the south having high Hg emissions stemming from high population densities and industrial activity.Hence episodic plumes of elevated Hg concentrations contrasted with low background concentrations are expected at these sites.

CMAQ model and simulation
A brief summary of the previous CMAQ model study (Gbor et al., 2007) will be provided here for reference.The model was based on CMAQ V4.3 (Bullock and Brehme, 2002), which was modified by including the dry deposition of Hg 0 and RGM (Gbor et al., 2006) and adding a model that calculates natural mercury emissions from soil, water and vegetation canopies.The simulation was conducted for the year 2002.The domain was the same as that used here.It covers most of North America using 132×90 grid cells with a horizontal spacing of 36 km and 15 vertical layers.Meteorology was provided by version 3.6 of the PSU/NCAR MM5 model with the PX land surface model (LSM) and indirect soil nudging (Xiu and Pleim, 2001).
Inventories of anthropogenic criteria pollutants (O 3 , NO 2 , particulate matter (PM), SO 2 , CO and Pb) for the United States were obtained from the 1999 National Emissions Inventory (NEI) version 3 IDA Files (US EPA, 2004a).The 1995 inventory of anthropogenic criteria pollutants for Canada was obtained from the Ontario Ministry of the Environment (OMOE) (A. Chtcherbakov, personal communication, 2003).Emissions of criteria pollutants from biogenic sources were processed using the BEIS3 program of the Sparse Matrix Operator Kernel Emission (SMOKE) modeling system with a gridded land use file for North America obtained from the OMOE (A. Chtcherbakov, personal communication, 2004).Anthropogenic Hg emission data for the US and Canada were obtained from the US EPA (2004b).
The natural mercury emissions were calculated using the emission model of Gbor et al., (2004Gbor et al., ( , 2006) ) and merged with the anthropogenic mercury and criteria emissions using SMOKE.Mercury boundary conditions were taken from the global mercury simulation of Seigneur et al. (2004), which included both natural and anthropogenic Hg emissions.The CMAQ model used these gridded emissions from SMOKE with the photolysis rate and initial and boundary conditions to simulate the atmospheric Hg concentrations on an hourly basis.For the purpose of comparison with measured TGM concentrations, modeled Hg 0 and RGM concentrations were summed and converted from ppmv to ng/m 3 using the temperature and pressure for each grid cell.

STILT model
The STILT model is built on source code from the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) model system (Draxler and Hess, 1998;Lin et al., 2003).In order to satisfy the well-mixed criterion in the strongly inhomogeneous environment of the PBL where the simple drift correction does not work (Lin et al., 2003;Thomson et al., 1997), the STILT model employed a reflection/transmission scheme for Gaussian turbulence instead.The parameterization for the PBL height was a modified Richardson number method that generalizes to unstable, neutral, and stable conditions (Lin et al., 2003;Vogelezang and Holtslag, 1996).The model simulates the transport of air parcels using ensembles of fictitious particles advected with mean wind velocities as well as stochastic velocities parameterized to capture the effects of turbulent transport.Further details of the STILT model can be found in Lin et al. (2003).
STILT is a receptor-oriented transport model that simulates tracer concentrations at a receptor and identifies upstream source regions based on a "footprint" concept.A footprint, f ( x r ,t r |x i ,y j ,t m ), in units of ppm/(µmole/m 2 /s), represents the sensitivity of the mixing ratio C( x r ,t r ) at receptor location x r at time t r to the surface flux F (x i ,y j ,t m ) from location x i ,y j at time t m .Thus it is a simulation of the mixing ratio at the receptor from a source of unit strength in each grid cell of the domain.The footprint is derived from the local density of particles by counting the number of particles (out of a total number N tot ) in surface-influenced boxes and determining the amount of time t p,i,j,k each particle p spends in each surface volume element (i,j,k) during each time step.The mathematical definition of a footprint (Lin et al., 2003) is given by.f ( concentrations were summed and converted from ppmv to ng/m 3 using the temperature and pressure for each grid cell.

STILT model
The STILT model is built on source code from the Hybrid Single Particle Lagrangian Integrated Trajectory (HYSPLIT) model system (Draxler and Hess, 1998;Lin et al., 2003).In order to satisfy the well-mixed criterion in the strongly inhomogeneous environment of the PBL where the simple drift correction does not work (Lin et al., 2003;Thomson et al., 1997), the STILT model employed a reflection/transmission scheme for Gaussian turbulence instead.
The parameterization for the PBL height was a modified Richardson number method that generalizes to unstable, neutral, and stable conditions (Lin et al., 2003;Vogelezang and Holtslag, 1996).The model simulates the transport of air parcels using ensembles of fictitious particles advected with mean wind velocities as well as stochastic velocities parameterized to capture the effects of turbulent transport.Further details of the STILT model can be found in Lin et al. (2003).
STILT is a receptor-oriented transport model that simulates tracer concentrations at a receptor and identifies upstream source regions based on a -footprint‖ concept.A footprint, ) , , | , ( Thus it is a simulation of the mixing ratio at the receptor from a source of unit strength in each grid cell of the domain.The footprint is derived from the local density of where m air is the molar mass of air.h is the height below which turbulent mixing is strong enough to mix the surface flux thoroughly, and ρ(x i ,y j ,t m ) is the average air density below h.Information about a footprint comes from computing the transport of an ensemble of particles backward in time using winds and turbulence statistics from meteorological fields.
The footprints can be integrated for different time periods and different areas depending on specific applications.The footprints can also be multiplied by surface fluxes to yield simulated concentrations.Through footprint elements f ( x r ,t r |x i ,y j ,t m ), STILT links the surface fluxes F (x i ,y j ,t m ) to concentration changes C m,i,j ( x r ,t r ) at a receptor with the following equation (Lin et al., 2003): Integrated Trajectory (HYSPLIT) model system (Draxler and 7 order to satisfy the well-mixed criterion in the strongly inh 8 PBL where the simple drift correction does not work (Lin et 9 the STILT model employed a reflection/transmission scheme 10 The parameterization for the PBL height was a modified 11 generalizes to unstable, neutral, and stable conditions (Li 12 Holtslag, 1996).The model simulates the transport of air par 13 particles advected with mean wind velocities as well as stoc 14 capture the effects of turbulent transport.Further details of t 15 Lin et al. (2003).16 STILT is a receptor-oriented transport model that simulates t 17 and identifies upstream source regions based on a -f 18 The parameterization for the PBL height was a 11 generalizes to unstable, neutral, and stable cond 12 Holtslag, 1996).The model simulates the transpor 13 particles advected with mean wind velocities as w 14 capture the effects of turbulent transport.Further 15 Lin et al. (2003) STILT calculates hourly concentrations by averaging the concentrations of all particles arriving at a receptor after a specific period.The relative importance of influences from pollution sources can be revealed by mapping footprints onto the emission inventory.

Treatment of Hg deposition
STILT was originally developed for atmospheric transport simulations of an inert tracer.Unlike inert tracers, Hg can be deposited to the surface, so we added a module in STILT to account for the effect of deposition on the Hg concentrations.For this purpose, the changes in atmospheric Hg concentration due to dry and wet depositions are expressed in terms of time constants: where C is the atmospheric concentration, t is time and β dry and β wet are time constants for dry and wet deposition respectively.
The time constant for dry deposition can be expressed as: where V dry (cm/s) is dry deposition velocity.Dry deposition velocities of Hg 0 and RGM were explicitly calculated by the Meteorology-Chemistry Interface Processor (MCIP) in the CMAQ simulation, and directly specified in the input.Z s (m) is the depth of the pollutant layer.Dry deposition is assumed to occur from the lowest model layer in CMAQ, which is 75 m.In order to be consistent with the CMAQ modeling, we set the surface layer depth Z s in STILT to 75 m for dry deposition calculations.The wet deposition of gases depends upon their solubility.For non-reactive gases it is a function of the Henry's Law constant.The gaseous wet deposition velocity can be defined as (Draxler and Hess, 1997): where R is the universal gas constant (0.082 atm/M/K) and T and P are, respectively, air temperature and precipitation rate in the grid box containing the particle.We used the Henry's Law constants 1.11×10 −1 M/atm for Hg 0 (Sanemasa, 1975) and 1.4×10 6 M/atm for RGM (Lindqvist and Rodhe, 1985).Gaseous wet removal only occurs for the fraction of the pollutant below the cloud top.The gaseous wet removal time constant is given by: where Z p is the depth of the meteorological layer in which the particle is found.F t is the fraction of the layer that is below the cloud top; it is explicitly determined by the MCIP in the CMAQ simulation.

Treatment of high stack Hg emissions
The concentration change C m,i,j ( x r ,t r ) of a tracer in the air parcel due to a surface emission F (x i ,y j ,t m ) (µmole/m 2 /s) is incremented whenever a parcel dips below a specific height h which is determined in STILT as a fraction of the planetary boundary layer height.STILT assumes that below the height h, surface flux is mixed thoroughly during one model time step and calculates the change in tracer concentration by vertically diluting the surface flux over h using a combination of Eqs. ( 1) and ( 2) C m,i,j ( x r ,t r ) = F (x i ,y j ,t m ) m air h ρ(x i ,y j ,t m ) where F (x i ,y j ,t m )m air /(h ρ(x i ,y j ,t m )) represents the dilution of the surface flux.This formulation does not apply for emission sources above the PBL.The Hg emissions from elevated sources were processed by SMOKE, which calculates the plume rise (Briggs, 1971(Briggs, , 1972) ) and thus provides the vertical distribution of the Hg elevated point-source emissions.In this process, we assume that the Hg emissions are mixed thoroughly in each grid cell during one model time step and estimate the dilution D(x i ,y j ,z k ,t m ) in each grid cell (i, j, k) that is higher than the height h, using its flux F (x i ,y j ,z k ,t m ): where ρ and L, the air density and the height of the grid box respectively, are obtained from the meteorological input data.Thus a concentration change C m,i,j,k ( x r ,t r ) of a tracer at a receptor due to emissions above the PBL is given by: = F (x i ,y j ,z k ,t m ) The chemical transformation of TGM is assumed to be unimportant and is not considered in the STILT modeling.This assumption is based on: (1) gas phase reactions of elemental Hg, which constitutes more than 90% of the Hg in the atmosphere (Slemr et al., 1985;Schroeder et al., 1991), are very slow; (2) although its chemical conversion is fast, RGM generally only accounts for less than 3% of TGM (Lindberg and Stratton, 1998;Poissant et al., 2005); and (3) the sixday particle simulation period is quite short compared to the atmospheric residence time of TGM.

Estimating the effect of transport errors
Errors in atmospheric transport lead to errors in simulated tracer concentrations.To investigate this in the present context, we employed an error analysis method reported by Lin and Gerbig (2005) and Gerbig et al. (2008) to determine uncertainties in modeled TGM concentrations caused by transport errors.According to this method, a transport error is introduced by incorporating wind field uncertainties ε into stochastic motions of air parcels.The uncertainties in winds are then propagated through stochastic motions of the air parcels in STILT.Transport error δ ε (C) in the modeled concentrations, C, can then be obtained simply from the square root of the difference between the variance of C in simulations with and without adding ε.The statistics of uncertainties ε in wind fields are determined by direct comparison of the Eta Data Assimilation System (EDAS) winds to radiosonde observations.These statistics include standard deviation in horizontal wind errors σ x ; their correlation time scales (l t ) and length scales in the horizontal (l x ) and vertical (l z ); the standard deviations in their mixed layer height errors (σ z i ); their correlation time scales (z i t ) and horizontal length scales (z i x ).The details of the method used to determine these statistics can be found in Lin and Gerbig (2005).We adopted the same values for the aforementioned error statistics as in Lin and Gerbig (2005) and Gerbig et al. (2008) in this paper, with the simplifying assumption that MM5 and EDAS meteorological fields yield, on average, similar errors as a first attempt to simulate the effect of transport error on the TGM concentrations.

STILT simulation
We used STILT to simulate the transport and deposition of Hg observed at the three monitoring sites for four periods (in UTC) during 2002: winter (18-27 February), spring (26 May-4 June), summer (20-29 July) and autumn (13-22 November).These simulation periods were chosen for seasonal coverage and because they contain episodes during which the Hg fluctuations were not simulated accurately by the original CMAQ Eulerian calculations.In the STILT simulations, ensembles of 3000 air parcels (hereafter called particles) were released from the three locations every hour.The choice of 3000 particles will be explained in Sect.3.1.These particles were run backward in time for a period of 6 days, which usually allowed them to reach lateral boundaries of the simulation domain (Fig. 1) -far from sources near the receptors -where they are assigned background concentrations from the global mercury simulation of Seigneur et al. (2004).For the small subset (approximately 9% of total number) of particles that did not reach the lateral boundary, we assigned Hg background concentrations at locations determined by extrapolating a line connecting the receptor to the particles based on the fact that these particles' endpoints were close to the boundary (430 km away, on average).Thus, even if extrapolated, they are still mostly found within the same coarse-grained gridcell of the global model (8 • latitude ×10 • longitude).Assigning the background concentrations for all particles was also necessary to establish the background contribution in this study to be compared against the natural and anthropogenic contributions (see below).STILT was driven by the same hourly MCIP output generated for the CMAQ modeling, after conversion from netCDF into NOAA ARL format to meet the STILT input format requirements.Dry deposition velocities for Hg 0 and RGM generated by the CMAQ modeling were also used.
To investigate the contributions from different sources, the STILT model was also run for each simulation period and receptor for three different scenarios as follows: 1. Background only: all Hg emissions were set to zero.
Only Hg background concentrations were included.
2. Natural only: Hg background concentrations and anthropogenic emissions were set to zero.Only natural Hg emissions were included.
3. Anthropogenic only: Hg background concentrations and natural emissions were set to zero.Only anthropogenic emissions were included.

Sensitivity to particle number
Due to the stochastic nature of STILT particle trajectories, the accuracy of the STILT simulation depends on the number of particles used.Theoretically, if an infinite number of particles were used, STILT would represent completely the ensemble properties of the transport to a given measurement location, assuming perfect parameterizations and input meteorology.A limited particle number leads to incomplete sampling of surface emissions and the magnitude of this sampling error varies with particle number.To quantify this error, we examined the standard deviation (i.e.sampling error) as a function of particle number for the simulated TGM concentrations at Egbert during the period from 20 to 29 July.The following particle numbers were examined: 50, 100, 200, 300, 500, 700, 1000, 2000, 3000, 4000 and 5000.The result (shown in Fig. 2) confirms that the error due to the stochastic nature of the model decreases with the number of particles used.The error varies most significantly for particle numbers less than about 1000, beyond which increases in particle numbers have only small effects on the sampling error.Since the model run time is proportional to the number of particles, we chose 3000 particles for the present simulations, which yielded a sampling error less than 6%.

Modeled and observed concentrations
The observed and modeled hourly Hg concentrations at all three sites are shown in Fig. 3 (winter and spring) and Fig. 4 (summer and autumn).Inspection of Figs. 3 and 4 shows that while STILT and CMAQ often showed correspondence to one another, CMAQ frequently lacked the elevated values and higher variability found in the measurements, which are better represented by STILT.
In order to quantify the models' abilities to capture concentration changes at different scales of variability, spectral densities of the Hg concentrations simulated by the two models were calculated using a Fourier transform method (Venables and Ripley, 2002).In Fig. 5, these are compared with observations between 13-22 November.The general characteristics for other periods are very similar to Fig. 5.Both STILT and CMAQ captured the low frequency variations, but there is a significant difference for the high frequencies (subdaily variations), where STILT is much closer to the observations.This shows that STILT has better skill than CMAQ in reproducing short term variations, such as those occurring in plumes.Because the same wind and emission fields were used by STILT and CMAQ, we believe that the difference stemmed from numerical diffusion is minimized in the Lagrangian model, preventing dilution of the footprint that dampens fluctuations in tracer concentrations.The artificial dilution is also reduced by the Lagrangian model by advecting particles backward in time from the receptor point rather than an entire gridcell.
In order to evaluate the predictions of TGM concentrations quantitatively from both modeling systems, three traditional statistical measures recommended by the US EPA (Doll, 1991) have been used: the mean normalized bias error (MNBE); mean normalized gross error (MNGE) and unpaired peak accuracy (UPA).These are defined in Table 1 and their calculated values for the present simulations are shown in Table 2.The simulation errors (MNGE and MNBE) from both modeling systems are smaller than 20%.CMAQ underpredicted the peak values of TGM plumes (indicated by a large negative value of UPA) in spring and summer periods, while STILT captured these plumes (UPA closer to 0), demonstrating the better performance of STILT for a period when local emissions (especially natural sources) are strong.There is no significant difference in simulation results between the two models for winter and autumn periods.
Performance of both models was especially poor at Point Petre, particularly in spring.Throughout all 4 periods, however, STILT showed larger errors at Point Petre in UPA, MNGE, and MNBE, although it mostly captured the timing and magnitude of the enhancements and the associated higher variability (Fig. 5).We speculate that the main source of error at Point Petre is the overestimation of Hg background concentrations, as reflected by the predominantly positive MNBE.Similarly, an underestimation of background concentrations may explain the underestimated concentrations at Burnt Island in the autumn period.Inaccuracies in local background concentrations are probably due to the coarse resolution (10 • ×8 • ) of the global mercury simulation by Seigneur et al. (2004).
The transport errors in STILT-modeled TGM concentrations are presented as error bars in Figs. 3 and 4. We are not aware of any currently available method to quantify transport errors in CMAQ.The relative importance of transport errors as a percentage of the modeled concentrations for all simulation episodes are listed in Table 2.These results show a clear seasonal variation: small in the winter and autumn episodes and approximately 3 times larger in spring and summer episodes.Larger transport errors are expected for spring and summer, when stronger natural sources (Sect.3.3) lead to more pronounced emission gradients that are improperly sampled by erroneous wind vectors (Lin and Gerbig, 2005).Comparing the transport errors to MNGE and MNBE, we can see that transport errors account for large portions of the discrepancies between the STILT simulations and the observed values.

Source contributions to the selected episodes
Figures 3 and 4 also show the contributions from natural (teal dashed curves) and anthropogenic (brown dashed curves) mercury sources to each site and period.The natural contributions are significant in spring and summer, with an average 0.24 ng/m 3 and small in winter and autumn, with an average 0.04 ng/m 3 .This reflects the strong seasonal variation of the natural emissions, which are dominated by vegetation and water sources.Natural emissions in winter and autumn are much lower due to smaller leaf area, lower temperature, weaker solar radiation and larger percentage of snow cover.For all simulation periods, the average observed Hg concentration is 1.71 ng/m 3 and the STILT-simulated Hg concentration is 1.78 ng/m 3 .The average contribution from natural sources is 0.14 ng/m 3 , almost twice as large as the anthropogenic sources.Roughly speaking, about 8% of the Hg comes from natural sources and 4% from anthropogenic sources.The contributions from natural sources are approximately equivalent to the anthropogenic emissions in winter and autumn.
Some episodic high Hg concentrations (plumes) at the three study sites were not captured well by the CMAQ simulation.To explore the possible causes for the failure of the Eulerian description in these cases, we carried out STILT simulations of the Hg plumes and also of temporally proximal periods in which the Hg concentration is low for each site and season.The time periods are designated in Figs. 3  and 4 by black dashed lines for the Hg plumes and black solid lines for the low episodes.STILT was used to calculate the footprints and identify the source regions for these episodes.
The footprints for the episodes observed at Burnt Island in February are displayed as an example in Fig. 6.In this and later Figures, significant Hg  cles.The major sources identified for the three measurements sites are all located in the region shown in Fig. 6, so we will discuss the results for this region in detail.Footprints, which are deduced solely from air parcel trajectories, are indicated by the colour contours and the scale at the bottom.These show the sources of the air parcels detected during the low and high episodes.
Footprints alone are insufficient to assess upstream source contributions at a receptor; the corresponding source fluxes are required as well.Figs.7-10 show the source contribution maps (derived by multiplying the footprint by the emissions) for all three sites as calculated by STILT for all episodes.These show that in most cases, the low concentration episodes occurred when the transport was from regions of low population density and where there are few anthropogenic mercury sources.The contributions from both anthropogenic and natural sources are negligible, with an average around 0.005 ng/m 3 in winter and autumn episodes.The background is the most important contributor at these times.In spring and summer, even for periods of low concentration, natural sources must be taken into account due to their large contributions, which average 0.13 ng/m 3 (maximum 0.2 ng/m 3 ).
Figures 7-10 show that episodes of high Hg concentration coincided with the arrival of air parcels from regions of higher population density in southern Ontario and north eastern United States, where there are many more sources of anthropogenic mercury.For the episodes of high concentrations, natural sources contribute much more than anthropogenic sources even for cases where there are many large anthropogenic sources upwind.This is illustrated by the dashed curves in the lower panels of Figs. 3 and 4. Natural and anthropogenic sources together contribute on average about 0.47 ng/m 3 to the Hg plumes, accounting for about 25% of the average observed Hg.

Summary and conclusions
Hg concentrations at three monitoring sites were simulated using the STILT Lagrangian particle transport model, which was modified for this study to deal with Hg deposition and high stack Hg emissions.The modeled Hg concentrations were compared with the observations, as well as with previously modeled results using the US EPA CMAQ model.While a comparison over longer time periods than the limited episodes in this study is necessary to establish unequivocal results, STILT-modeled air concentrations of Hg generally agreed well with observations and, on average, exhibited better performance than the Eulerian CTM CMAQ for the examined episodes.In particular, STILT reproduced the high frequency variability present in the data, which was not present in the CMAQ results (Fig. 5).Since the same meteorological fields and Hg emission inventory data were used as inputs for STILT and CMAQ, the better performance of STILT can likely be ascribed to its ability to capture nearfield influences by treating the receptor as a point rather than a gridcell and by minimizing numerical diffusion.A unique strength of STILT is its ability to estimate quantitatively and account for "transport errors" -i.e., the impact of uncertainties in the driving wind fields on simulated tracer concentrations.To our knowledge, such a capability is missing from CMAQ, as well as other Lagrangian models.The transport errors estimated in this study can reach approximately 10% of the simulated value (Table 2) and account for large fractions of the total errors.These results underscore the importance of taking transport errors into consideration in simulation studies.Future studies will improve the transport error estimates by carrying out an independent estimate of transport errors that will compare simulated wind fields against radiosonde observations rather than adopting error statistics based on previous comparisons, which was a simplification adopted for the purposes of the current work.
Upwind contributions to Hg concentrations at the three sites were also simulated using STILT for several selected episodes.Simulation results show that the major contributions to observed low Hg concentrations at the three sites were from some sparsely populated regions toward the north of the three sites and that the major sources of observed high concentrations were generally from more heavily populated areas with large point emission sources.Natural mercury sources contribute more than anthropogenic sources to the observed concentrations.
Finally, we note that a unified framework that combines a Lagrangian model (STILT) with an Eulerian model (CMAQ) is a particularly powerful approach for this kind of simulation.While this paper has mostly focussed on contrasting the two approaches, their strengths are actually complementary.The Eulerian approach yields a full three-dimensional picture of pollutant concentrations, including chemical transformations, over the entire model domain and simulates baseline conditions well, but misses much of the short-term variation.The Lagrangian approach (currently with no chemistry) identifies the more localized source regions of the designated receptors, which enables the researcher to "zoom into" specific measured pollution events to identify their origins.

Fig. 1 .
Fig. 1.Simulation domain and the locations of the three monitoring sites: Burnt Island (B), Egbert (E) and Point Petre (P).Annual average Hg emission rate including anthropogenic and natural Hg for 2002.
in units of ppm/(µmole/m 2 /s), represents the sensitivity of the mixing ratio ) it is a simulation of the mixing ratio 21 unit strength in each grid cell of the domain.The footprint is 22 particles by counting the number of particles (out of a total n 23 boxes and determining the amount of time Δt p,i,j,k each pa 24 volume element (i, j, k) during each time step.The mathema 25 et al., 2003) is given by.26 r ) = f ( Integrated Trajectory (HYSPLIT) model system (D 7 order to satisfy the well-mixed criterion in the st 8 PBL where the simple drift correction does not wo 9 the STILT model employed a reflection/transmissi 10

Fig. 2 .
Fig. 2. Dependence of STILT simulation on the number of particles.

Fig. 3 .
Fig. 3. Total gaseous mercury (TGM) concentration comparisons among observed (red), STILT modeled (blue), and CMAQ modeled (green) for two simulation periods: 18-27 February (left) and 26 May-4 June (right), 2002 for Burnt Island (top), Egbert (middle) and Point Petre (bottom).Background concentrations derived from the lateral boundary condition is in violet.The anthropogenic (brown) and natural (teal) contributions to the simulated concentration are shown as dashed curves in the lower panels.Black dashed lines in the vertical designate a selected high Hg episode; black solid lines designate a selected low Hg episode.

Fig. 4 .
Fig. 4. Total gaseous mercury (TGM) concentration comparisons among observed (red), STILT modeled (blue), and CMAQ modeled (green) for two simulation periods: 20-29 July (left) and 13-22 November (right), 2002 for Burnt Island (top), Egbert (middle) and Point Petre (bottom).Background concentrations derived from the lateral boundary condition is in violet.The anthropogenic (brown) and natural (teal) contributions to the simulated concentration are shown as dashed curves in the lower panels.Black dashed lines in the vertical designate a selected high Hg episode; black solid lines designate a selected low Hg episode.

Fig. 5 .
Fig. 5. Spectral densities of observed and simulated TGM concentrations for the simulation period of 13-22 November 2002.

Fig. 6 .
Fig. 6.Modeled footprints [log 10 (ppm/(µmole/m 2 /s))] for the selected Hg plume (right) and low episode (left) in 18-27 February for Burnt Island site (green diamond).Important Hg point sources are designated with circles and the greyscale shadings of the circles represent their emission strengths.
Figures 3 and 4 also show the contributions from natural (teal dashed curves) and anthropogenic (brown dashed curves) mercury sources to each site and period.The natural contributions are significant in spring and summer, with an average 0.24 ng/m 3 and small in winter and autumn, with an average 0.04 ng/m 3 .This reflects the strong seasonal variation of the natural emissions, which are dominated by vegetation and water sources.Natural emissions in winter and autumn are much lower due to smaller leaf area, lower temperature, weaker solar radiation and larger percentage of snow cover.For all simulation periods, the average observed Hg concentration is 1.71 ng/m 3 and the STILT-simulated Hg concentration is 1.78 ng/m 3 .The average contribution from natural sources is 0.14 ng/m 3 , almost twice as large as the anthropogenic sources.Roughly speaking, about 8% of the Hg comes from natural sources and 4% from anthropogenic sources.The contributions from natural sources are approximately equivalent to the anthropogenic emissions in winter and autumn.Some episodic high Hg concentrations (plumes) at the three study sites were not captured well by the CMAQ simulation.To explore the possible causes for the failure of the Eulerian description in these cases, we carried out STILT simulations of the Hg plumes and also of temporally proximal periods in which the Hg concentration is low for each site and season.The time periods are designated in Figs.3 and 4by black dashed lines for the Hg plumes and black solid lines for the low episodes.STILT was used to calculate the footprints and identify the source regions for these episodes.The footprints for the episodes observed at Burnt Island in February are displayed as an example in Fig.6.In this and later Figures, significant Hg point sources (≥0.04 tonnes/year) are designated by circles.Their emission strengths are indicated by the greyscale shading of the cir-

Fig. 7 .
Fig. 7. Source contributions [log 10 (ppm)] derived by multiplying the footprint with gridded emissions for the selected Hg plumes (right) and low episodes (left) in 18-27 February for Burnt Island (top), Egbert (middle) and Point Petre (bottom).Important Hg point sources are designated with circles and the greyscale shadings of the circles represent their emission strengths.Note the different scales between the plumes (right) and low (left) episodes.

Table 1 .
Definition of US EPA recommended statistic parameters: UPA, MNGE and MNBE.: prediction at time i.O i : observation at time i.N : total number of observations.P u peak : maximum predicted concentration.O peak : maximum observed concentration. i

Table 2 .
Statistics for predicted TGM concentrations using CMAQ and STILT.
i : transport error in prediction at time i.P i : prediction at time i.