Comparison between Back-Trajectory Based Modeling and

Feb 2, 2005 - Matthew J. Dellinger , Michael P. Ripley .... Ralf Ebinghaus , Daniel Engstrom , Xinbin Feng , William Fitzgerald , Nicola Pirrone , Eri...
0 downloads 0 Views 601KB Size
Download PDF
Environ. Sci. Technol. 2005, 39, 1715-1723

Comparison between Back-Trajectory Based Modeling and Lagrangian Backward Dispersion Modeling for Locating Sources of Reactive Gaseous Mercury YOUNG-JI HAN AND THOMAS M. HOLSEN* Department of Civil and Environmental Engineering, Clarkson University, Potsdam, New York 13699-5710 PHILIP K. HOPKE Department of Chemical Engineering and Center for Air Resources Engineering and Science, Clarkson University, Potsdam, New York 13699-5708 SEUNG-MUK YI School of Public Health and the Institute of Health and Environment, Seoul National University, 28 Yeongon-dong, Chongro-gu, Seoul 110-799, Korea

Reactive gaseous mercury (RGM) was measured using an annular denuder coated with potassium chloride at three rural sites (Potsdam, Stockton, and Sterling) in New York State from April 2002 to April 2003. Concentrations of RGM ranged from 0.1 to 84.6 pg m-3 with large spatial and temporal variation. Potential source contribution function (PSCF), a common receptor modeling tool, was used with these measurements, and source-receptor relationships were calculated using back-dispersion and deposition as well as back-trajectories. Modeling results were compared with the RGM emissions inventory, and Spearman rank-order correlation coefficients were calculated. PSCF results incorporating backward dispersion and deposition were better correlated with the emissions inventory than PSCF based on back-trajectories alone. This difference was determined to be mainly due to the inclusion of dispersion rather than deposition. The main sources of RGM were suggested to be coal-fired power plants in New York and Pennsylvania, the large copper smelter in Quebec, and the taconite mining areas around the Great Lakes.

Introduction Mercury can exist in three oxidation states: 0, +1, and +2. Elemental mercury (Hg0) constitutes the majority of mercury in the atmosphere and is relatively long-lived, having a residence time of 0.5-2 yrs, while divalent mercury (Hg2+) has a very short residence time (∼days) (1, 2). Some of the divalent mercury such as HgCl2 or HgBr2 is highly water soluble, reactive, easily deposited, and less volatile than Hg0, so is often called reactive gaseous mercury (RGM). For example, ionic forms make up about 40% of the Hg emitted * Corresponding author phone: (315)268-3851; fax: (315)268-7636; e-mail: [email protected]. 10.1021/es0498540 CCC: $30.25 Published on Web 02/02/2005

 2005 American Chemical Society

from utility boilers, but it is quickly deposited to the surface by wet and dry processes or tends to condense onto atmospheric particulate matter. Therefore the contribution of RGM in total deposition to surfaces is probably very high even though its concentration is typically less than 5% of the total gas-phase mercury (3, 4). Trajectory-based receptor models, including potential source contribution function (PSCF), have been proven to be useful tools in identifying source areas (5). However, in the boundary layer, single trajectories are hardly sufficient to describe transport processes due to turbulent mixing and may benefit from the inclusion of dispersion (6). In addition, for reactive species or species that undergo dry and wet deposition, atmospheric mechanisms including reaction and dry and wet deposition should be considered when identifying the source-receptor relationship. In this regard, total gas phase mercury, which consists mostly of elemental mercury, is of little concern because elemental mercury is not reactive and not usually removed readily by deposition. However, the deposition velocity of RGM is extremely high for both wet and dry mechanisms. Therefore, receptor modeling approaches using only backward trajectories may not correctly identify source locations. To date, there have been few attempts to identify source locations using source-receptor relationships that include atmospheric dispersion, reaction, and dry and wet deposition. One exception is quantitative transport bias analysis (QTBA). However, dispersion in QTBA is treated very simply and may not adequately describe atmospheric turbulence conditions (7). There have been a few studies that obtained sourcereceptor relationships using the Lagrangian particle dispersion model (LPDM) (8-10), but this approach has not been used for locating RGM sources. This study aims to (i) identify RGM source areas which impact New York State using backtrajectory and back-dispersion based hybrid receptor modeling and (ii) determine the effects of the atmospheric dispersion and deposition on RGM source area identification using receptor models.

Data Source and Model Description Ambient RGM Data. RGM samples (from midnight to midnight) were taken using an annular denuder (purchased from URG) coated with KCl every third day at the Potsdam (44.75 N, 75.0 W) and Stockton sites (42.27 N, 79.38 W) in NY from June 2002 to April 2003. In Sterling (43.34 N, 76.66 W), RGM samples (from 9 am to 9 am of the next day) were collected every sixth day starting in April 2002 (Table 1). All three sites are located in remote rural areas and are not significantly affected by local RGM sources (Figure 1). After sampling, denuders were heated to convert Hg2+ to Hg0 and analyzed using a cold vapor atomic florescence spectrophotometer (CVAFS; model 2500 from Tekran). Detailed procedures are described in Han (11) and Han et al. (12). The RGM concentrations measured in this study were in a reasonable range for rural sites (Figure 2). According to previous investigations, the RGM concentration can increase up to 200 pg m-3 in a 1-h sampling period near mercury sources (13-16). However, Malcolm et al. (17) found the average RGM concentrations to be 1.6 ( 1.4 pg m-3 in the marine regime and 4.9 ( 4.0 pg m-3 in the mixed regime (wind direction from both marine and urban area) at the background site in Pompano Beach, Florida, which are in the lower range of concentrations measured in this study. Emissions Data. The latest emissions inventory for mercury sources in the United States (1996 and 1999) and Canada (1995) was used in this study. A mercury emissions VOL. 39, NO. 6, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1715

TABLE 1. Summary of Ambient RGM Data (pg m-3) Used in This Study sampling site

sample number

sampling period

Potsdam, NY

77

Stockton, NY

59

Sterling, NY

62

June 2002 April 2003 June 2002 April 2003 April 2002 April 2003

arithmetic geometric mean ( SD mean ( SD 4.2 ( 6.4

2.5 ( 2.8

5.7 ( 9.2

3.4 ( 2.6

6.0 ( 10.8

3.4 ( 3.1

inventory for the United States was obtained from the U.S. EPA (18, 19). For coal-fired electricity generation, municipal waste incinerators, and medical waste incinerators, the estimates in this inventory were for 1999, while the remainder were representative of 1996. For Canada, a 1995 inventory was used (19). Mercury speciation based on very few measurements for most source categories was done for the U.S. inventory, but the speciation data for the Canadian emissions were not available. It was simply assumed that the proportion of the mercury forms emitted from Canadian sources was the same as the corresponding source categories in the United States. Only direct anthropogenic emissions were included in this inventory. Lagrangian Trajectory Model - HYSPLIT 4. The HYSPLIT 4 model developed by the National Oceanic and Atmospheric Administration (NOAA) that was used in this study is a hybrid between Eulerian and Lagrangian approaches, in which advection and diffusion calculations are made in a Lagrangian framework while concentrations are calculated on a fixed grid. There are other Lagrangian models including FLEXPART (20) and STILT (21). In this section, the basic concept of HYSPLIT 4 modeling system is discussed, and the description of particle dispersion will be described later. The initial version of Hysplit used only rawinsonde observations for meteorological data, but in the most recent version, Hysplit 4, the rawinsonde data were replaced by gridded meteorological data from either analyses or shortterm forecasts from routine numerical weather prediction models (22). Precipitation and surface fluxes are taken from forecast files only. The meteorological data fields are linearly interpolated to a terrain-following coordinate system. The Eta Data Assimilation System (EDAS) having a horizontal resolution of 80 km, twenty three vertical levels up to 2250 hPa, and a temporal resolution of 3 h was used as the meteorological data set. At a minimum Hysplit requires U, V (the horizontal wind components), T (temperature), Z (height), or P (pressure), and the pressure at the surface, P0. Also, in most circumstances, the vertical velocity field (W) is relative to the meteorological model’s native terrain-following coordinate system. Once the basic (U, V, W) meteorological data are processed and interpolated to the internal model grid, trajectories can be computed by advection components. The advection of a particle is computed from the average of the three-dimensional velocity vectors for the initial-position P(t) and the first-guess position P′(t + ∆t). The velocity vectors are linearly interpolated in both space and time. The first guess position is

P′(t + ∆t) ) P(t) + V(P,t) ∆t

(1)

and the final position is

P(t + ∆t) ) P(t) + 0.5[V(P,t) + V(P′, t + ∆t)] ∆t (2) where, V, here, represents the average of the threedimensional velocity vectors. 1716

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 6, 2005

To compute the vertical advection component, the model vertical velocity mode which uses the meteorological model’s vertical velocity fields (kinematic) generated by a diagnostic or prognostic meteorological model was used. It should be noted that the well-mixed criterion can be violated due to inconsistencies in both the turbulent and mean wind transport components, so that the forward and backward simulations would not give the same result. If forward particles were started in a low-turbulence region at the source, and the receptor was in a high-turbulent region, then few forward time particles would reach the receptor. We did not simulate forward trajectories to compare the backward trajectories for this specific project; however, in HYSPLIT, the errors between forward and backward trajectories are usually acceptable for typical applications (http:// www.arl.noaa.gov/faq/hg5.html). Potential Source Contribution Function (PSCF). Potential source contribution function (PSCF) was originally developed by Ashbaugh et al. (23) and has been used in many applications to locate source areas (24-26). It counts each trajectory segment endpoint that terminates within a given grid cell. If N is the total number of trajectory segment endpoints over the study period and if n segment trajectory endpoints fall into the ijth cell, then the probability of these events (P[Aij]) is calculated by nij/N. If mij is the number of segment endpoints in the same ijth cell when the concentrations are higher than a criterion value, the probability of this high concentration event, Bij, is given by P[Bij] ) mij/N. The PSCF value for the ijth cell is then defined as a conditional probability, PSCF(i,j) ) P[Bij]/P[Aij]. Therefore, cells having a high PSCF value are likely to contribute to high concentration events at receptor sites so they are reasonably assumed to be possible source areas. Even though the trajectory segment endpoints calculated by trajectory models such as HYSPLIT (hybrid single particle Lagrangian integrated trajectory) have uncertainties derived by both integration and interpolation, a sufficient number of samples (sufficient number of endpoints) can accurately locate source areas if the errors are random and not systematic. Conventional PSCF with Source-Receptor Relationship using Back-Trajectories. The source-receptor relationship can be calculated using either a trajectory model or a dispersion model (eq 3). In this study, the HYSPLIT 4 model was used to model both cases as

y ) Mx

(3)

where y indicates the discrete observation (e.g., measured concentration) at receptor sites, x is the source emission term, which varies over locations and times, and M is the source-receptor relationship and includes the transport process. In conventional PSCF, the number of end points in the ij-th cell, which is expressed as a residence time (hr), is simply considered to be the source-receptor relationship of the ij-th cell. PSCF with Source-Receptor Relationship using Backward Dispersion HYSPLIT. The PSCF values were computed using trajectories including only an advection component and did not consider dispersion, reaction, and/or deposition. However, in a Lagrangian dispersion model such as HYSPLIT 4, a random component to the motion caused by atmospheric turbulence can be added in each step to the advective motion. In this way, a cluster of particles released at the same point

FIGURE 1. Three monitoring sites of RGM in New York State.

FIGURE 2. Daily RGM concentrations measured at three monitoring sites in New York from April 2002 to April 2003. will expand in space and time simulating the dispersive nature of the atmosphere (22, 27) given by

Xfinal(t + ∆t) ) Xmean(t + ∆t) + U′(t + ∆t)∆tG

(4)

Zfinal(t + ∆t) ) Zmean(t + ∆t) + W′(t + ∆t)∆tZtop-1

(5)

where Xmean and Zmean indicate the mean horizontal and vertical position of a particle, respectively. The contributions of the turbulent wind components (U′, horizontal; W′, vertical) are added to the “mean” position to give a final position from which the advection in the next time step is computed (27, 28). Because the horizontal and vertical positions (X and Z) are given in grid and sigma units, respectively, while the turbulent velocity components are in m sec-1, the unit conversion factors (G and Ztop) are required. The horizontal turbulent velocity components at the current

time U′(t + ∆t) are from the turbulent velocity component at the previous time U′(t), an autocorrelation coefficient (R) that depends on the exponential of the time step, the Lagrangian time scale, and a computer-generated random component (′′).

U′(t + ∆t) ) R(∆t)U′(t) + U′′(1 - R(∆t)2)0.5

(6)

The vertical turbulent velocity is defined as

W′/σW (t + ∆t) ) R(∆t)(W′(t)/σW) + (W′′/σW) (1 - R(∆t)2)0.5 + TLW(1 - R(∆t)) ∂σW/∂z (7) where, σW(t + ∆t) ) σW(t) + W′(t)∆t∂σW(t)/∂z; R(∆t) ) exp(-∆t/TLi), where TLi is either TLw or TLu; U′′(or W′′) ) σi λ where λ is a Gaussian random number with mean of 0 and a standard deviation of 1; σi is the standard deviation of the VOL. 39, NO. 6, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1717

turbulent velocities; and TLi is the Lagrangian time scale (i ) u, horizontal; or i ) w, vertical). During this computation step, reaction and dry/wet deposition can be included. For computational simplicity, the total deposition from both dry and wet removal processes is expressed in terms of (reciprocal) time constants (27).

Dwet+dry ) m{1 - exp[∆t(βdry + βgas + βinc + βbel)]}

(8)

where m is the pollutant particle mass, βdry is a dry deposition term, βgas is for wet removal for gases, βinc is for in-cloud wet removal of particles, and βbel is for below-cloud wet removal of particles. For gaseous pollutants, dry deposition can be directly obtained using a dry deposition velocity or computed using a resistance model (29).

Vd ) [Ra + Rb + Rc]-1

(9)

where Ra is the aerodynamic resistance, Rb is the quasilaminar sublayer resistance, and Rc is the canopy layer resistance. Computation of Ra and Rb depend on meteorological conditions (30, 31) and Rc depends primarily upon a number of surface physiological and the chemical/physical properties of the pollutant (32). Wet deposition is divided into two processes, those in which the polluted air is continuously ingested into a cloud from a polluted boundary layer and those in which rain falls through a polluted layer (29). For a gaseous pollutant, the wet deposition velocity depends on its Henry’s Law constant (M atm-1) and precipitation rate.

Vw,gas ) HRTP

(10)

where R is the universal gas constant (0.082 atm M-1 K-1), T is temperature, and P is precipitation rate. In the final modeling step, the total mass of particles after consideration of mechanisms such as reaction and deposition is summed in the ijth grid cell and divided by the grid-cell volume, so that gridded concentration fields (mass m-3) are computed. To obtain the residence time of each grid (the so-called transmission-corrected residence time), a transformation is necessary (10).

t)

∆TsVs jc µtot

(11)

where ∆Ts is the time during which the source is acting, Vs is the volume of the grid cell, jc is the concentration in a grid cell produced by HYSPLIT, and µtot is the total mass associated with particles relased. Once the transmission-corrected residence time is calculated, the potential source contribution function (PSCF) can be calculated from P[Aij] and P[Bij]. The HYSPLIT dispersion model can estimate the atmospheric turbulence component from calculated stability (which changes at every grid and time step) through a meteorological model and relatively precise dry and wet deposition can be calculated using the meteorological model and measured or estimated precipitation rate (27). Model Application. Two-day backward trajectories and back-dispersion were used since RGM is rapidly deposited after emission and therefore does not undergo long-range transport. Trajectories were calculated every 6 h for all 24hr-averaged samples, therefore four trajectories represented one sample. PSCF used a grid-cell size of 1° by 1° and arrival heights of 100 m. A maximum trajectory height was 5000 m, and a minimum height of 0 m (above ground level) was used. In the backward dispersion HYSPLIT model, all RGM was assumed to be HgCl2. The emission rate was set to a continuous 1 kg hr-1 to obtain the source-receptor relation1718

9

ship. A total of 10,000 particles was released from the receptor site at each simulation to well-represent the dispersion effect in the boundary layer. Averaged concentrations over 48-hr intervals were calculated for each grid cell (1° by 1°) and an output grid level of 50 m was used. A starting height of 100 m above ground level was used in all cases. The molecular weight of HgCl2, surface reactivity ratio, and Henry’s constant were set to be 271.5 g, 1, and 1.4 × 106 M atm-1 (3), respectively. Since HgCl2 is very reactive and sticky to many surfaces (19), the surface reactivity ratio was assumed to be 1. To calculate dry and wet deposition of RGM, gas diffusivity (Da) and an effective Henry’s constant were needed. Since Da for HgCl2 is not available in the database (33), it was calculated to be 4.5 × 10-2 cm2 s-1 using the equation cited in U.S.EPA (3).

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 6, 2005

Da )

1.9 (MW)2/3

(12)

There are a number of aqueous phase reactions affecting the effective Henry’s constant for HgCl2 such as reactions with Cl-, SO3-, and OH-. To determine an appropriate Henry’s constant, a special version of the NOAA HYSPLIT 4 model was used to simulate the concentrations of HgCl2 in aqueous and gas phases. (A detailed description of this special version HYSPLIT is described in Cohen et al. (19)). After generating Hg2+ concentrations about 3000 times using different combinations of chloride ion and S(IV) concentrations and pH, the average value of the effective Henry’s constant was estimated to be 1.2 × 108 M atm-1 (34). Dry deposition was calculated by either directly using a dry deposition velocity or using the three-layer resistance model as is discussed below.

Results and Discussions Conventional PSCF using Back-Trajectories. First, PSCF was performed using 2-day back-trajectories for both single-site and multisite measurements to identify the RGM source locations. Multisite PSCF has the advantage that it may remove the trailing effect that can identify areas upwind and downwind of sources as a source area (11, 35). In this paper only multiple-site PSCF results will be presented but singlesite PSCF results are described in Han (12). For multisite measurements, joint-probability PSCF (JP-PSCF) is often used. In this method, samples are separated into “high” and “low” concentrations using criteria values (in this case, the arithmetic mean) for each site. The trajectories are then pooled using eq 13.

P(JPij) )

P[Bij]potsdam + P[Bij]stockton + P[Bij]sterling P[Aij]potsdam + P[Aij]stockton + P[Aij]sterling

(13)

To remove the large uncertainty caused when a grid cell has a small number of total endpoints and large PSCF values, an arbitrary weight function W(nij) was applied when the number of the end points in a particular cell was less than about three times the average number of the end points for all cells (12, 36, 37) (eq 14). The average number of end points per cell (avg. ) 7.54) was calculated by multiplying the total sample number (n ) 198) by the number of trajectories per sample (m ) 4) and the number of end points in each trajectory (i ) 48). This product was then divided by the number of cells (j ) 5040).

[

1.00 0.70 W(nij) ) 0.42 0.17

25 < nij 10 < nij e 25 5 < nij e 10 nij e 5

(14)

FIGURE 3. JP-PSCF result for RGM measurements taken in Potsdam, Stockton, and Sterling, NY. Forty-eight-hours back-trajectories were used for calculating source-receptor relationships.

FIGURE 4. Locations and emissions of coal-fired power plants in the United States (data from (19)). JP-PSCF identified source areas in southwestern New York, eastern Ohio, Pennsylvania, southern Illinois, New Jersey, and Indiana (Figure 3). Pennsylvania and southeastern Ohio (the Ohio River valley) have some of the largest coal-fired power plants in the United States (Figure 4; 19, 38). According to the National Wildlife Federation (38), the C. R. Huntley Steam Station near Buffalo is the largest New York State source and emits an estimated 274 lb of mercury into the air annually.

These sources were well identified by JP-PSCF. There are numerous waste incinerators and oil-fired power plants as well as coal-fired power plants in the area around St. Louis, MO (3, 18, 19), and those sources were also determined to significantly contribute to high RGM concentrations in New York State. Back-Dispersion-Based PSCF. The PSCF was also calculated for each site using the transmission-corrected VOL. 39, NO. 6, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1719

FIGURE 5. JP-PSCF result for RGM measurements taken in Potsdam, Stockton, and Sterling, NY. Forty-eight-hours backward dispersion was used and dry and wet deposition were considered for calculating source-receptor relationships. source-receptor relationship obtained using HYSPLIT 4 in the back-dispersion mode. Wet and dry deposition for RGM was considered during transport, and the dry deposition velocity was calculated using the three-layer resistance model. Once the transmission-corrected residence time was calculated, PSCF was calculated using the multisite measurements. The transmission-corrected residence time (hr) in the grid cells was significantly smaller than the residence time obtained using the back-trajectory model, because deposition was considered. A weight function was applied when a grid cell had a small residence time (20% of average residence time in a grid cell; τij e 0.1 h; eq 15).

[

1.00 0.70 W(nij) ) 0.42 0.17

0.1 (h) < nij 0.06 < nij e 0.1 0.03 < nij e 0.06 nij e 0.03

(15)

The sensitivity of the results to different weight functions was investigated and similar results were obtained (12). Results for PSCF using transmission-corrected residence time incorporating dispersion and deposition were somewhat different from results using conventional PSCF (Figure 5). The main source areas were identified to be coal-fired power plants located in New York State, a large smelter located in Quebec, the urban areas of New Jersey and Philadelphia, PA, and the taconite mining area around Lake Superior (Figure 5). Generally, the source areas in JP-PSCF results using transmission-corrected residence time were located close to the receptor sites, while in JP-PSCF using back-trajectories they were located in the southern United States. Sensitivity Test for the Back-Dispersion-Based Model PSCF. To determine the importance of atmospheric dispersion on receptor modeling, PSCF using only back-dispersion without deposition was performed. The difference in PSCF results considering only dispersion and considering both dispersion and deposition will indicate the importance of considering deposition in identifying source areas for RGM. A sensitivity test between the different approaches used to calculate dry deposition was performed as well. In this 1720

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 6, 2005

analysis, PSCF results using the three-layer resistance model with an effective Henry’s constant of 1.2 × 108 M atm-1 and using a constant dry deposition velocity of 0.016 m s-1 were obtained. Input of a constant dry deposition velocity assumes that the deposition flux equals the velocity times the groundlevel air concentration. In the three-layer resistance model, the dry deposition velocity was calculated based on landuse (region type) and atmospheric conditions in each time step. The sensitivity test result for PSCF using multisite RGM measurements is shown in Figure 6. The top, middle, and bottom panels show the PSCF results with only backward dispersion (no-deposition), with both dispersion and deposition calculated using the three-layer model, and with both dispersion and deposition calculated assuming a constant dry deposition velocity, respectively. PSCF results when deposition was and was not included were very different from each other. However, similar results were obtained when deposition was calculated using the three-layer model and when a constant dry deposition velocity was assumed. When comparing the PSCF result based on backtrajectories (Figure 3) with results based on dispersion and deposition (Figure 5), several differences can be noted. To obtain a more quantitative comparison between modeling results, Spearman rank-order correlation coefficients (rs) were calculated (Table 2). Case 1 (middle panel of Figure 6) is the PSCF result incorporating dispersion, wet deposition, and dry deposition calculated using the three-layer resistance model, Case 2 (bottom panel of Figure 6) is the PSCF result incorporating dispersion, wet deposition, and dry deposition calculated using a constant dry deposition velocity, Case 3 (top panel of Figure 6) is the PSCF result incorporating only dispersion (without dry/wet deposition), and Case 4 (Figure 3) is the conventional PSCF result. All results were determined to be statistically significantly correlated with each other, even Cases 1 and 4. Correlation coefficients between PSCF values calculated by incorporating dispersion and deposition calculated using either the three-layer (Case 1) or a constant dry deposition velocity method (Case 2) is very high, indicating that the technique used to calculate dry deposition

TABLE 2. Spearman Rank-Order Correlation Coefficients (rs) between PSCF Results for RGM case 1 JP-PSCF

case 1 case 2 case 3 case 4

0.87 0.78 0.38

case 2

case 3

case 4

0.87

0.78 0.75

0.38 0.39 0.40

0.75 0.39

0.40

TABLE 3. Spearman Rank-Order Correlation Coefficients (rs) and the Spatial Correlation Index (r) between Modeling Results and Emissions Inventory JP-PSCF using two-day back-trajectories JP-PSCF incorporating back-dispersion JP-PSCF incorporating dispersion and deposition

FIGURE 6. Sensitivity analysis of JP-PSCF based on back-dispersion. Top, middle, and bottom panels indicate the PSCF results with only turbulence, with deposition calculated with the three-layer model, and with deposition calculated with a constant deposition velocity. did not change modeled results significantly. However, the correlation coefficients between Cases 1 and 3 and between Cases 2 and 3 are smaller, indicating that the consideration of deposition is important in the model result. PSCF results incorporating dispersion (Cases 1 and 2) are not as well

rs

r

0.34 0.51 0.39

45.87 55.54 47.69

correlated with the conventional PSCF result (Case 4), although they are still statistically correlated (at a significant level of 0.01). These lower correlation coefficients suggest that dispersion had a stronger influence than deposition in PSCF modeling of RGM. Statistical Comparisons between Emission Rates and PSCF Results for RGM. PSCF results for RGM were compared statistically with a recent RGM emissions inventory (3, 18, 19). However, it should be noted that this emissions inventory is likely to be even less accurate for Hg2+ than for Hg0. Some large mercury sources (for example, taconite facilities in Minnesota and Michigan) are not included in the current emission inventory, and Canadian sources may be particularly inaccurate (34). Also, the atmospheric reaction mechanisms were not considered in calculating source-receptor relationship at this time, and this omission of reaction can cause another uncertainty to the modeling result. Inclusion of reaction in receptor modeling will be investigated in future work. Additional problems in this comparison are the small number of samples and limited time period covered by the RGM sampling campaign. For these reasons, a good correlation between modeling results and the emissions inventory was not expected and higher correlations do not necessarily suggest better modeling results. Since the distribution of the estimated PSCF result values and emission rates is far from normal, the Spearman correlation coefficient rs was used (Table 3). Emission rates in the domain used in receptor modeling for RGM were ranked and zero emission fields were eliminated to avoid ties in ranking. Correlations were poor as expected, however all PSCF results were statistically correlated with the emissions inventory (at significant level of 0.05). Correlation between PSCF results incorporating dispersion and the emissions inventory was higher than when conventional PSCF result was used. The effect of including deposition is not clear and the correlation was smaller than for PSCF results incorporating only dispersion. Other possible reasons for these lower correlation coefficients in addition to those mentioned earlier include the possibility that the HYSPLIT model does a poor job of estimating deposition. The finding that PSCF incorporating dispersion is better correlated with the emissions inventory than conventional PSCF is consistent with single site results for the Potsdam, Stockton, and Sterling sites (12). Because the exact residence time of RGM in the atmosphere is unknown, 5-day back trajectories were also used in PSCF modeling. The 2-day and 5-day results were very similar although the results using 5-day back-trajectories suggested more regional sources (in Missouri, Iowa, Arkansas, and Nebraska) than did results using 2-day back-trajectories. When compared with the emissions inventory using Spearman rank-order correlation within the modeling domain, VOL. 39, NO. 6, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1721

2-day back-trajectory results were consistently better correlated with the emissions inventory for all sites than were the 5-day results (12). This result suggests that the extended backward simulation time of 120 h may not be necessary and that the regional source areas identified when one assumes a long-range transport for RGM may not be reliable. Another technique of spatial correlation index, r, was used to measure the association between two maps (39). In this analysis, a general spatial autocorrelation statistic (Γ) between two spatial data matrixes g and f is represented by two functions including a measure of the spatial (Gij) and a measure of nonspatial proximity (Cij).

Γ)

∑G C

(16)

ij ij

ij

For location i, xi, yi is the coordinate and gi and fi are the emission rate and modeling result, respectively. Gij and Cij, which are the spatial proximity of locations i and j, and the proximity of i and j, are defined as follows:

1 [(xi - xj)2 + (yi - yj)2 + 1]

(17)

[(gi - gj)(fj - hf ) + (gi - gj)(fi - hf )] 2

(18)

Gij )

Cij )

where gj and hf are the average values of g and f, respectively. The null hypothesis (H0) indicating that there is no spatial correlation between g and f was tested under a randomization assumption. The first two moments of Γ are

ER(Γ) ) S0T0/n(n - 1)

(19)

ER(Γ2) ) S1T1/2n(2) + (S2 - 2S1)(T2 - 2T1)/(4n(3) + S02 + S1 - S2)(T02 + T1 - T2)/n(4) (20) where,

S0

∑G ; ij

S1 )

ij

S2 )

1 2

∑(G

ij

+ Gji)2

ij

∑ ) (∑ G + ∑G ) ; 2

ij

i

i

ji

j

n(b) ) n(n - 1)(n - 2)‚‚‚(n - (b + 1)) Tk is the same as Sk with Gij replaced by Cij. The spatial correlation index, r, is obtained by the following equation under N(0,1) distribution.

r ) (Γ - E(Γ))/(E(Γ2) - (E(Γ))2)0.5

(21)

When the null-hypothesis is tested at a significance level of 0.05, the rejection region is [-∞, 1.6] for one tailed test. In this study, the spatial correlation index, r, has same trend as the Spearman correlation coefficient, showing that the correlation between PSCF result incorporating dispersion and the emissions inventory is higher than when the conventional PSCF result is used (Table 3). When deposition was included, the r became smaller than the PSCF result incorporating only dispersion. In summary, the identified source areas using conventional PSCF are considerably different from those incorporating transmission-corrected residence time. Using the HYSPLIT model, trajectories are estimated considering only advection. To consider atmospheric dispersion, threedimensional turbulence terms were added to the position estimated using advection. Unless HYSPLIT has serious problems in estimating meteorological turbulence fields, 1722

9

ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 39, NO. 6, 2005

PSCF result incorporating dispersion is expected to perform better than conventional PSCF. According to the Spearman rank-order correlation coefficients, PSCF results based on backward dispersion are better correlated with the emissions inventory than conventional PSCF. Since the emissions inventory of mercury is not likely to be accurate and there are only a limited number of samples in this study, the effect of dispersion and deposition on hybrid receptor modeling should be verified using a more accurate emission inventory and a sufficient number of measurements over a longer time period in future.

Acknowledgments We sincerely acknowledge the following colleagues for their help with sampling: Dr. Michael Milligan and Chris Andolina at SUNY Fredonia, James Pagano and Lauren Falanga at SUNY Oswego, Dr. Wei Liu currently at Georgia Tech, and Soon-Onn Lai and Keun-Wook Lee at Clarkson University. We also thank Dr. Mark Cohen at NOAA for sharing his knowledge and sparing his valuable time to discuss these results. This work is funded by New York State Energy Research and Development Authority under Contract 6083, U.S. EPA Regions 2 and 5, and the Great Lakes National Program Office (GLNPO).

Literature Cited (1) Slemr, F.; Schuster, G.; Seiler, W. Distribution, speciation and budget of atmospheric mercury. J. Atmos. Chem. 1985, 3, 407434. (2) Tokos, J. J.; Hall, B.; Calhoun, J.; Prestbo, E. Homogeneous gasphase reaction of Hg0 with H2O2, O3, CH3I, and (CH3)2S: Implications for atmospheric Hg cycling. Atmos. Environ. 1998, 32, 823-827. (3) U.S. EPA. Mercury study report to Congress. Office of Air Quality Planning and Standards and Office of Research and Development; EPA-452/R-97-005; U.S. Government Printing Office: Washington, DC, 1997. (4) U.S. EPA. Locating and estimating air emissions from sources of mercury and mercury compounds (Mercury L&E); Final draft report; EPA-454/R-97-0121997; Research Triangle Park, NC, 1997. (5) Hopke, P. K.; Gao, N.; Cheng M. D. Combining chemical and meteorological data to infer source areas of airborne pollutants. Chemom. Intell. Lab. Syst. 1993, 19, 187-199. (6) Stohl, A. Computation, accuracy and applications of trajectoriess A review and bibliography. Atmos. Environ. 1998, 32 (6), 947966. (7) Keeler, G. J. Ph.D. Thesis, University of Michigan, Ann Arbor, MI, 1987. (8) Seibert, P. Inverse modeling with a Lagrangian particle dispersion model: Application to point releases over limited time intervals. InSchiermeier, F. A., Gryning, S. E., Eds.; Air Pollution Modeling and its Application XIV; Kluwer Academic: Norwell, MA, 2001; pp 381-190. (9) Stohl, A.; Eckhardt, S.; Forster, C.; James, P.; Spichtinger, N.; Seibert, P. A replacement for simple back trajectory calculations in the interpretation of atmospheric trace substance measurements. Atmos. Environ. 2002, 36, 4635-4648. (10) Seibert, P.; Frank, A. Source-receptor matrix calculation with a Lagrangian particle dispersion model in backward model. Atmos. Chem. Phys. Discuss. 2003, 3, 4515-4548. (11) Han, Y. J.; Holsen, T.; Hopke, P.; Yi, S. M.; Pagano, J.; Milligan, M.; Lai, S. O.; Liu, W.; Falanga, L.; Andolina, C. Atmospheric gaseous mercury concentrations in New York State: relationships with meteorological data and other pollutants. Atmos. Environ. 2004, 38, 6431-6446. (12) Han, Y. J. Ph.D. Thesis, Clarkson University, Potsdam, NY, 2003. (13) Lindberg, S. E.; Stratton, W. J. Atmospheric mercury speciation: concentrations and behavior of reactive gaseous mercury in ambient air. Environ. Sci. Technol. 1998, 32, 49-57. (14) Stratton, W. J.; Lindberg, S. E.; Perry, C. J. Atmospheric mercury speciation: Laboratory and field evaluation of a mist chamber method for measuring reactive gaseous mercury. Environ. Sci. Technol. 2001, 35, 170-177. (15) Landis, M.; Keeler, G. Atmospheric mercury in the Lake Michigan basin: Influence of the Chicago/Gary urban area. Environ. Sci. Technol. 2002, 31, 4508-4517.

(16) Landis, M.; Keeler, G. Atmospheric mercury deposition to Lake Michigan during the Lake Michigan Mass Balance Study. Environ. Sci. Technol. 2002, 31, 4518-4524. (17) Malcolm, E. G.; Keeler, G. J.; Landis, M. S. The effects of the coastal environment on the atmospheric mercury cycle. J. Geophys. Res. 2003, 108, D12. doi: 2002JD003084. (18) Ryan, R. Point, area, and mobile sources of mercury prepared for U.S. EPA analysis of proposed multi-pollutant Power Generation Legislation. U.S. EPA, Office of Air Quality Planning and Standards, Emissions, Monitoring and Analysis Division, Emissions Factors and Inventory Group; U.S. Government Printing Office: Washington, DC, 2001. (19) Cohen, M.; Artz, R., Draxler, R.; Miller, P.; Niemi, D.; Ratte, D.; Deslauriers, M.; Duvar, R.; Laurin, R.; Slotnick, J.; Nettesheim, T.; McDonald, J. Modeling the atmospheric transport and deposition of mercury to the Great Lakes. Environ. Res. 2004, 95 (3), 247-265. (20) Stohl, A.; Forster, C.; Eckhardt, S.; Spichtinger, N.; Huntrieser, H.; Heland, J.; Schlager, H.; Wilhelm, S.; Arnold, F.; Cooper, O. A backward modeling study of intercontinental pollution transport using aircraft measurements. J. Geophys. Res. 2003, 108, D12, 4370. doi: 10.1029/2002JD002862. (21) Lin, J. C.; Gerbig, C.; Wofsy, S. C.; Andrews, A. E.; Daube, B. C.; Davis, K. J.; Grainger, C. A. A near-field tool for simulating the upstream influence of atmospheric observations: The Stochastic Time-Inverted Lagrangian Transport (STILT) model. J. Geophys. Res. 2003, 108, 16, 4493. doi: 10.1029/2002JD003161. (22) Draxler, R. R.; Hess, G. D. An overview of the HYSPLIT_4 modelling system for trajectories, dispersion, and deposition. Australian Met. Magzine 1998, 47, 295-308. (23) Ashbaugh, L. L.; Malm, W. C.; Sadeh, W. D. A Residence time probability analysis of sulfur concentrations at Grand Canyon National Park. Atmos. Environ. 1985, 19, 1263-1270. (24) Russell, T.; Cass, G. Acquisition of regional air quality model validation data for nitrate, sulfate, ammonium ion and their precursors. Atmos. Environ. 1984, 18, 1815-1827. (25) Zeng, Y.; Hopke, P. K. A study on the sources of acid precipitation in Ontario, Canada. Atmos. Environ. 1989, 23, 1499-1509. (26) Gao, N.; Cheng, M. D.; P. K. Hopke, Receptor modeling for airborne ionic species collected in SCAQS. Atmos. Environ. 1994, 28, 1447-1470.

(27) Draxler, R. R.; Hess, G. D. Description of the HYSPLIT_4 modeling system. NOAA Technical Memorandum ERL ARL-224; National Oceanic and Atmospheric Administration: Washington, DC, 1997. (28) Fay, B.; Glaab, H.; Jacobsen, I.; Schrodin, R. Evaluation of Eulerian and Lagrangian atmospheric transport models at the Deutscher Wetterdienst using ANATEX surface tracer data. Atmos. Environ. 1995, 29, 2485-2497. (29) Hicks, B. B. Aerosols: Research, Risk Assessment, and Control Strategies; Lewis Publishers: Chelsea, MI, 1986. (30) Seinfeld, J. H.; Pandis, S. N. Atmospheric Chemistry and Physics; John Wiley & Sons: New York, 1998. (31) Wesely, M. L.; Hicks, B. B. Some factors that affect the deposition rates of sulfur dioxide and similar gases on vegetation. J. Air. Pollut. Control Assoc. 1977, 27, 1110-1116. (32) Wesley, M. L. Parameterizations of surface resistances to gaseous dry deposition in regional-scale numerical models. Atmos. Environ. 1989, 23, 1293-1304. (33) U.S. EPA, Soil screening guidance: Technical background document and User’s guide, EPA/540/R-96/018; Office of Solid Waste and Emergency Response: Washington, DC, 1996. (34) Cohen, M. Personal communication, NOAA Air Resource Laboratory. (35) Hsu, Y. K., Ph.D. Thesis, Clarkson University, Potsdam, NY, 2001. (36) Polissar, A. V.; Hopke, P. K.; Harris, J. M. Source regions for atmospheric aerosol measured at Barrow, Alaska. Environ. Sci. Technol. 2001, 35, 4214-4226. (37) Polissar, A. V.; Hopke, P. K.; Poirot, R. L. Atmospheric aerosol over Vermont: Chemical composition and sources. Environ. Sci. Technol. 2001, 35, 4604-4621. (38) National Wildlife Federation. Cycle of harm: Mercury’s pathway from rain to fish in the environment; Great Lakes Natural Resource Center: Ann Arbor, MI, 2003. [email protected]. (39) Haining, R. Spatial Data Analysis in the Social and Environmental Sciences; Cambridge University Press: Cambridge, U. K., 1990.

Received for review January 28, 2004. Revised manuscript received November 22, 2004. Accepted December 3, 2004. ES0498540

VOL. 39, NO. 6, 2005 / ENVIRONMENTAL SCIENCE & TECHNOLOGY

9

1723