Article pubs.acs.org/est
Temporal Trends of Persistent Organic Pollutants: A Comparison of Different Time Series Models Marta Venier,† Hayley Hung,‡ Włodzimierz Tych,§ and Ronald A. Hites†,* †
School of Public and Environmental Affairs, Indiana University, Bloomington, Indiana 47405, United States Air Quality Research Division, Environment Canada, Toronto, ON, M5P 1W4, Canada § Lancaster Environment Centre, Lancaster University, Lancaster, LA1 4YQ, U.K. ‡
S Supporting Information *
ABSTRACT: We analyzed the vapor phase atmospheric concentrations of representative persistent chemicals (i.e., α- and γ-hexachlorocyclohexane, phenanthrene, PCB-18 and PCB-52) in samples collected at a remote site near Eagle Harbor, Michigan, and at an urban site in Chicago, Illinois, using four time series models: a modified Clausius−Clapeyron equation, a multiple linear regression that includes both a linear and an harmonic dependence on time, digital filtration (DF), and dynamic harmonic regression (DHR). The results of these different models were evaluated in terms of goodnessof-fit, long-term trends, and halving times. The four approaches all provided highly significant descriptions of the data, with coefficients of determination (R2) ranging from 0.33 to 0.96. In general, the DF and DHR methods fit the data better, capturing not only the seasonal variations of the atmospheric concentrations but also smaller scale interannual variations in the long term trends. The halving times calculated using the four methods were generally similar to one another, and they ranged from about 4 years for γ-HCH at Chicago to about 60 years for PCB-52 at Chicago. This analysis showed that each of these four statistical methods for evaluating long-term time series has advantages and disadvantages. The choice of the appropriate method should depend on the output needed, the type of audience, and the availability and usability of the necessary software.
■
INTRODUCTION The Integrated Atmospheric Deposition Network (IADN) is a joint project between the U.S. Environmental Protection Agency and Environment Canada. IADN has been in operation since 1990 as a response to Annex 15 of the Great Lakes Water Quality Agreement that required the United States and Canada to provide regular updates on the atmospheric deposition of toxic chemicals to the lakes. Several persistent organic pollutants are measured every 12 days at five sites around the U.S. shores of the lakes. The compounds measured include 22 organochlorine pesticides (such as DDTs, chlordanes, and hexachlorocyclohexanes), 84 polychlorinated biphenyl congeners or groups of congeners (PCBs), 17 polycyclic aromatic hydrocarbons (PAHs), and 43 flame retardants. Given the number of samples collected since the beginning of this project and given the number of chemicals analyzed in each sample, this project has gathered about 700 000 atmospheric concentration measurements. As these data have been acquired, we have used progressively larger data sets in several publications, but there have been several common threads of these studies: First, IADN data have been used to determine when the atmospheric concentrations of these chemicals maximize during the year in each phase. For example, lower molecular weight compounds are relatively © 2012 American Chemical Society
abundant in the vapor phase, and thus, their concentrations are influenced by atmospheric temperature and maximize during the warmer summer months. Second, the IADN data have been used to study the spatial distribution of these chemicals in the Great Lakes region. In general, a significant correlation between local human population density and atmospheric concentrations has been found,1,2 but in some instances, localized sources affect the spatial distribution of chemicals. For example, the atmospheric concentrations of Dechlorane Plus maximize near its production source. 3 Third and perhaps most importantly, IADN data have been used to determine the rate at which atmospheric concentrations have decreased (if at all) as a function of time. For example, PCBs were banned in 1977, but their atmospheric concentrations are now declining slowly.2 Conversely, the concentrations of some other compounds, notably α- and γ-hexachlorocyclohexane (HCH) have decreased rapidly.1 This information, which is the focus of this paper, is helpful for predicting the behavior of atmospheric concentrations in the future, especially for regulatory purposes. Received: Revised: Accepted: Published: 3928
September 2, 2011 February 28, 2012 March 7, 2012 March 20, 2012 dx.doi.org/10.1021/es204527q | Environ. Sci. Technol. 2012, 46, 3928−3934
Environmental Science & Technology
Article
the IADN Quality Control Project Plan.7 Data below the detection limit were considered as nondetects (empty cells). Data Analysis. Because environmental concentration data are log-normally distributed, the gas phase concentration data (in pg/m3) for each chemical were first converted into their natural logarithms. These data were then fitted using four different approaches: (a) The first was a modified Clausius−Clapeyron equation (in this case, the partial pressure of the compound was replaced by its concentration)8
In the scientific literature, several different statistical approaches have been used to determine the rates at which atmospheric concentrations of persistent organic pollutants (POPs) change. While these approaches all share the common purpose of separating seasonal signals from long-term trends, they are quite different. They range from the relatively simple Clausius−Clapeyron correction of the vapor phase partial pressures4 to the relatively advanced autocorrelation of residuals.5 In this paper, we concentrated on four widely used statistical approaches for the analysis of POPs time trends: (a) a modified Clausius−Clapeyron equation, (b) a multiple linear regression that includes both a linear and a harmonic dependence on time, (c) digital filtration, and (d) dynamic harmonic regression. We have applied them to selected compounds representative of their chemical categories: α-HCH, γ-HCH, phenanthrene, PCB-18, and PCB-52. In each case, we have analyzed the gas phase concentrations of these five compounds at two representative IADN sites (one urban site at Chicago, Illinois (IL) and one remote site at Eagle Harbor, Michigan (MI)). The ultimate goal of this exercise was to compare these four approaches in terms of the information they provide and their ease of use.
ln C = a0 − a1/T + a2t
(1) 3
where C is the gas phase concentration in pg/m , T is the average atmospheric temperature at the sampling site during the sampling period in degrees Kelvin, t is time expressed in Julian days starting from 1 January 1990, a0 is an intercept that rectifies the units, and a1 and a2 are the coefficients that describe the dependence of the concentration on temperature and time, respectively. Note that a2 is a first order rate constant with units of days−1, which can be used to calculate a halving time in years (t1/2 = −ln(2)/ 365.25 a2). Among the four methods applied here, this approach is the one most closely related to a physical process: The a1 coefficient is directly related to the enthalpy of surface-air partitioning which describes the partitioning of the chemical between air and an unknown component of the earth’s surface (i.e., soil, vegetation, and water).9 While this approach provides a direct relationship with a physical process, it requires knowledge of the average atmospheric temperature during the sampling event and whether equilibrium has been attained, knowledge which is not always available. Equation 1 was fitted to the data using Minitab 15, which gave the coefficients and their standard errors. The long-term trend was directly calculated using coefficient a2, which is equivalent to calculating a linear regression of the model fit with time. (b) The second approach was a multiple regression equation2
■
EXPERIMENTAL SECTION Sampling and Analytical Methodology. In this study, air samples collected at two of the six United States IADN sites were used: Eagle Harbor, MI, and Chicago, IL. The first is representative of a very remote site (only about 500 people live within a 25 km radius of this site), and the second is representative of an urban site (about 3.6 million people live within a 25 km radius of this site). Data from the entire period of operations were used. The dates were 1990 to 2009 (inclusive) for the Eagle Harbor site and 1996 to 2009 (inclusive) for the Chicago site. Details of the sample collection, extraction, and analysis procedures can be found elsewhere,6 and only a brief description will be presented here. A modified, Anderson high-volume air sampler (General Metal Works) was used to collect air samples for 24 h every 12 days at a flow rate giving a total sample volume of ∼820 m3. The gas phase was collected on Amberlite XAD-2 resin (Supelco, 20−60 mesh) held in a 9 cm o.d., stainless steel cartridge. After sampling, the XAD was Soxhlet extracted for 24 h with a 1:1 acetone hexane mixture. Prior to extraction, surrogate standards were spiked into the sample. The extract was reduced in volume by rotary evaporation, the solvent was exchanged to hexane, and this solution was fractionated on a column containing 3.5% water deactivated silica gel. The column was eluted with 25 mL of hexane and then with 25 mL of a 1:1 hexane dichloromethane mixture. After blow down with a gentle stream of nitrogen, the sample was spiked with the internal standards. PCBs and pesticides were analyzed by gas chromatography on Hewlett-Packard 5890 and 6890 instruments equipped with 63Ni electron capture detectors and DB-5 and DB-1701 (J&W Scientific) 60 m columns (250 μm i.d. and 0.1 μm film thickness). Quantitation was done using the internal standard method. Surrogate standards were used to estimate recoveries of each compound in each sample. PAHs were analyzed by gas chromatographic mass spectrometry (GC/MS) on an Agilent 6890 GC with a 5973 mass spectrometer, using a DB-5 (J & W Scientific) 30 m column (250 μm i.d. and 0.25 μm film thickness).6 Quality control and quality assurance procedures were followed to ensure data accuracy. The detailed procedures are described in the IADN Quality Assurance Program Plan and in
ln C = b0 + b1sin(zt ) + b2cos(zt ) + b3t
(2)
where z = 2π/365.25 (this fixes the periodicity at one year), b0 is an intercept that rectifies the units, b1 and b2 are harmonic coefficients that describe the seasonal variations of the concentrations with time, and b3 is a first-order rate constant. Equation 2 was fitted to the data using Minitab 15, which gave the coefficients and their standard errors. This method is a simple harmonic regression with a seasonally varying component (represented by b1 and b2) and a long-term exponential component (represented by b3). As above, the parameter b3 can be converted to a halving time. This regression has recently been applied to a superset of the data reported here.1 This approach does not require any assumptions or knowledge of atmospheric temperature. Its main advantage is that it allows an exact determination of the specific date when the concentrations maximize or minimize. The long-term trend was calculated as described above for method (a) using instead coefficient b3. (c) The third approach used the digital filtration (DF) method.10 For each compound, an approximate longterm trend and an average seasonal cycle were determined by fitting a smoothing Reinsch-type cubic spline and Fourier components to the data, respectively, in an 3929
dx.doi.org/10.1021/es204527q | Environ. Sci. Technol. 2012, 46, 3928−3934
Environmental Science & Technology
Article
Captain toolbox for MATLAB.13 The DHR model allows for these variations by estimating the time-varying amplitudes, ai,t and bi,t and αi,t , and βi,t described in eqs 4 and 5, in the data and then removing these components to display the underlying trend. The observation disturbance noise component, et , also interpreted as the irregular component, includes a variety of other, individually unobserved, processes such as the erratic incursion of polluted air masses, the impact of local sources, or even variations in removal processes, such as precipitation events. The harmonics identified by the DHR were 12, 6, and 4 months, roughly equivalent to 30, 15, and 10 samples/annual cycle. The first base annual frequency is derived from IADN’s sampling frequency (sampling once every 12 days results in about 30 samples per year), whereas the other two represent 6- and 4-month harmonic cycles, respectively. The harmonics occurring over shorter periods ( 0.05). At the other extreme, the plots in Figure 2 for γ-HCH at Eagle Harbor indicate that these coefficients are highly significant (p < 0.01). This suggests that the PCB-52 concentrations at Chicago are changing slowly, if at all, and that the γ-HCH concentrations at Eagle Harbor are decreasing rapidly. These same long-term features are also seen in the DF and DHR results from these data, but in these cases, more information is available on subtle changes over the years. Nevertheless, the approaches encapsulated by eqs 1 and 2 are useful when a relatively simple representation of the time-series data is needed; both approaches are particularly useful when the goal of the analysis is to easily determine if the concentrations have significantly increased or decreased at a given site. Both the DF and DHR approaches give a more nuanced depiction of the long-term trends, capturing also shorter time scale interannual variations. For example, both the DF and DHR approaches indicate that the atmospheric concentrations of PCB-52 at Chicago (see Figure 1), have not been steadily decreasing since 1996; rather, after reaching a minimum around 2003−2004, the concentrations started to increase again, reaching a maximum in 2007, followed by another decrease.
Figure 3. Top: Seasonally adjusted concentration measurements for PCB-52 in Chicago together with the trend estimate (red line) and its uncertainty band equivalent to two standard errors (dashed black lines) calculated using DHR. Bottom: Rate of change or slope of the trend depicted in the top panel with its uncertainty band. The trend is significant only when the slope and its uncertainty do not include the zero. The gray bands identify the time periods where the trend was statistically significant.
adjusted data (i.e., the raw data after removing the seasonality), which included the long-term trend. One particularly useful aspect of both the DF and the DHR approaches is that they also provide the rate of change or slope (i.e., the derivative of the trend), see Figure 3 and SI Figure S4. In addition, DHR inherently provides the uncertainty of the changing slope, which shows when a concentration increase or decrease became statistically significant (see Figure 3, bottom). Only when the slope’s confidence interval (shown as a gray band) did not include zero was the trend considered to be statistically significant. For PCB-52 at Chicago, the trend increased significantly from May 2004 until November 2005 (see the first gray band in Figure 3). The concentrations then decreased significantly after August 2007 (see the second gray band in Figure 3). The end of this last interval should be considered cautiously because, at the extremities of the data distribution, the uncertainties of the estimates are higher than in the middle of the distribution (as in any regression analysis). This abrupt change in concentration was not observed for PCB-52 at Eagle Harbor or for PCB18 at either site (SI Figure S3). After excluding possible experimental changes or unusual deviations from the protocols, we speculate that the 2007 maximum in PCB-52 concentrations in Chicago may be related to the demolition of an old building and the construction of a new condominium complex near the IADN sampling site. Such construction would have disturbed the local soil, releasing a bolus of PCBs to the atmosphere. Work on the interior and exterior of the new building continued until December 2006. In addition, in the summer of 2005, the exterior finish from another building 200 m west of the sampling site was removed using sand blasting, which could have released PCBs from the 3932
dx.doi.org/10.1021/es204527q | Environ. Sci. Technol. 2012, 46, 3928−3934
Environmental Science & Technology
Article
caulk used in the exterior façade of this older building. A similar increasing trend was observed also for other PCB congeners (i.e., 101 and 180) and for total PCBs but not for PCB 18. This hypothesis is strengthened by the observation of somewhat higher levels of PAHs in Chicago, perhaps caused by an increased number of heavy duty construction vehicles in the area, around the same time (see SI Figures S2 and S5). The DF and DHR approaches also identified an interesting feature in the data for γ-HCH concentrations in Chicago (Figure 2). Both approaches showed that the concentrations of γ-HCH did not steadily decline. The concentrations of γ-HCH in Chicago were somewhat stable from 1996 to the beginning of 2002, at which time, they declined rapidly until 2007, when they tapered off. In the DHR approach a significantly negative slope of the trend indicated that this decrease was statistically significant between January 2002 and October 2005. γ-HCH, also known as lindane, was introduced as a replacement for technical HCH (which contained 60−70% α-HCH, 5−12% β-HCH, and 10−15% γ-HCH) in the early 1980s. The use of γ-HCH was gradually restricted, and by 2002, it was used only for the treatment of barley, corn, oats, rye, sorghum, and wheat. γ-HCH was finally deregistered in the U.S. on October 4, 2006, and by October 1, 2009, it was not used anymore.15 In Canada, the use of γ-HCH followed a similar timeline; its use for the treatment of various agricultural seeds was phased out between 2002 and 2004. The sharp overall declines in the atmospheric concentrations of γ-HCH in Chicago and Eagle Harbor (Figure 2) are most likely a reflection of these restrictions. A similar acceleration in the declining trend was also observed in Arctic air at Alert, Canada, and at Zeppelin, Norway.16 These rapid decreases are examples of how restrictions on the use of a compound can have a positive outcome in the environment in a relatively short time. Comparison of Halving Times. The halving time is defined as the time required for the atmospheric concentration of a given substance to drop by a factor of 2. This time provides useful insights concerning the effectiveness of legislative measures and control actions. Given that the purpose of this study was a comparison of the four time trend analysis approaches, we will not discuss in depth the environmental relevance of the calculated halving times nor will we compare them to previously calculated ones. Rather, we will discuss the differences and/or similarities among the values obtained using these four approaches. In general, there is good agreement among the halving times calculated with the four different methods; see Table 1. The differences were around 1−3%, with a few exceptions. The halving times of PCB-18 were around 10 years at Eagle Harbor and 8 years at Chicago, while those of PCB-52 were around 18 years at Eagle Harbor and >50 years at Chicago. In fact, for PCB-52 at Chicago, the halving times calculated with eqs 1 and 2 were not statistically significant. As described earlier, the high halving times for PCB-52 at Chicago were likely driven by the increase in concentration centered in 2007; clearly, more data are needed to determine if the PCB-52 concentrations at Chicago have started to decrease again at a similar rate as before. In any case, we suggest that a halving time in excess of 50 years represents no real decrease in these atmospheric concentrations. There were no significant differences in the halving times of α- and γ-HCH calculated using the four different methods. In fact, these two compounds were the ones for which we observed the smallest differences among the four halving times (
