Luis A. de P. Vasconcelos

Department of Chemistry, Washington University, St. Louis, Missouri

Jonathan D. W. Kahl and Desong Liu

Edward S. Macias, and Warren H. White

Department of Chemistry, Washington University, St. Louis, Missouri

Abstract. The potential of routine back trajectory analyses to locate sources of contaminants in air at Grand Canyon is investigated with data on methylchloroform (CH₃CCl₃) concentrations collected during the Subregional Cooperative Electric Utility, Department of Defense, National Park Service, and Environmental Protection Agency Study(SCENES). Following a now-standard approach, the distribution of back trajectory segment endpoints over a latitude-longitude grid is examined as a function of measured concentrations at the fixed monitoring site. Grid cells in which segment endpoints are preferentially associated with high concentrations are then identified as candidate emissions sources. The method correctly identifies southern California as a source of CH₃CCl₃. An objective rationale is developed for screening out chance associations, attributable to statistical fluctuations, and the results are evaluated for real and hypothetical tracers with known, simple distributions.

The identification of sources (or source regions) and their contribution to the ambient concentrations of particles and gases is an essential part of visibility protection programs. It is a difficult task, and current source identification and apportionment methods carry with them a high degree of uncertainty [National Research Council, 1993]

One approach to source identification and apportionment uses statistical techniques to extract patterns of empirical association between the trajectories and concentrations of arriving air. Such methods are sometimes grouped generically as residence time analyses [White et al., 1994]. The type of analysis examined here was developed by Ashbaugh et al. [1985] to identify source areas influencing sulfur concentrations at Grand Canyon National Park. Bresch et al. [1987] and Gebhart and Malm [1991] extended the analysis to other national parks and also included carbonaceous species in their study. Their analysis was presented in terms of spatial probability distributions for air parcels arriving at varying particle sulfur and carbonaceous species concentrations; it will be described here in terms of time densities and conditional frequencies. Poirot and Wishinski [1985] used a similar approach to study sulfates and their impact on summertime visibility in northern Vermont. More recently, Zeng and Hopke [1989] and Gao et al. [1994] have used this approach in a variety of applications.

Copyright 1996 by the American Geophysical Union.

Paper number 95JD02609.

0148-0227/96/95JD-02609$05.00

Our study employs measurements of methylchloroform (CH₃CCl₃) near Grand Canyon to test the usefulness of con-ditional frequency analysis. Emissions inventories [Sheiman et al., 1990] and other analyses [White et al., 1990; Miller et al., 1990] show regional emissions of this long-lived halocarbon to be heavily concentrated in the Los Angeles basin. Our objective is to provide a "calibration" of conditional frequency analysis on this conservative species of known origin, which can then be used to evaluate future analysis with species that are nonconservative or have uncertain source distributions.

Chemical measurements were made from late 1984 to fall 1989 during the Subregional Cooperative Electric Utility, Department of Defense, National Park Service, and EPA Study (SCENES) on visibility [Mueller et al., 1986; Vasconcelos et al., 1994]. This paper focuses on methylchloroform (CH₃CCl₃) data taken at the Meadview observatory (Meadview, Ariz.; 36°01'N-114°04'W; 900 meters above mean sea level (msl)). These measurements were described by Miller et al.[1990]. Hourly concentrations of CH₃CCl₃ are averaged here according to the filter samples so that future analysis with the particle data can be comparable with this calibration. Several sampling schedules and protocols were followed during SCENES. Filter samples were collected during different epochs for 8-hour, 12-hour, 16-hour, and 24-hour periods. Only the SCENES summer quarter (June, July, and Aug.) is considered in this calibration. Future communications will extend the results for other seasons, variables, and locations.

The Air Resources Laboratory Atmospheric Transport and Dispersion (ARL-ATAD) model [Heffter, 1980] was used to estimate 3-day backward trajectories arriving at Meadview 4 times daily. Input to the model consists of vertical profiles of wind, temperature and pressure measured at National Weather Service rawinsonde stations. Air parcel transport is simulated as two-dimensional advection using horizontal winds averaged throughout the mixed layer, which is determined using the temperature and pressure profiles. The mixed layer is variable and may reach up to 3000 m in depth. The main output of the model is a set of 24 latitude-longitude coordinates that give the estimated position of an air parcel every 3 hours during the 72 hours prior to its arrival at the receptor. Occasionally, trajectories terminate prematurely or are missing entirely due to gaps in the input meteorological data set or to advection out of the computational domain. In this study, incomplete trajectories are considered invalid.

One could hope to locate the sources of a conservative contaminant by tracing back the trajectory of each infinitesimal volume of air that arrived containing above-background concentrations of that contaminant. If a large number of different back trajectories were observed, their intersections would indicate the positions of possible sources. However, this idealized approach would require high spatial resolution of the trajectories and high temporal resolution of the tracer concentrations; in reality, our information from wind data and filter samples is relatively coarse. As in many other problems that deal with spatial and temporal distributions, we have to characterize continuous variables through finite sets of discrete samples.

A grid with 0.5° longitude by 0.5° latitude cells is superimposed over the region defined by 130ºW-95ºW and 25ºN-50ºN. The ATAD trajectories used here describe parcels arriving in intervals of 6 hours, at 0000, 0600, 1200, and 1800 UT. They are based on spatially averaged winds, and need not track any actual volume element. The position of a parcel is logged every 3 hours by a "time stamp", which represents one 3-hour exposure of one parcel to a grid cell’s emissions. A parcel arriving at 1200 UT is assumed representative of air arriving between 0900 UT and 1500 UT, the hypothetical volume element following the mean transport layer wind is assumed representative of the entire diffusing puff, and the coordinates of this centerline element 36 hours upwind are assumed representative of its location between 34.5 and 37.5 hours. Similarly, the average tracer concentration measured during the sampling interval is assumed representative of the concentration during every instant of that time interval.

Let C denote the set of times during which contaminant concentrations were successfully measured, and let A denote the set of arrival times for which trajectory calculations were attempted. (In this paper, A is the discrete set {t: t = 0000, 0600, 1200, or 1800 UT, t Î June, July, or August of 1985 - 1989}.) Only the subset V Ì A of times for which full 72-hour trajectories were calculated will be considered. For any subset C* Ì C of concentration measurement times, simple bookkeeping yields the number G_ij(AÇ VÇ C*) of time stamps in each grid cell (i,j) from all complete trajectories arriving at the receptor during C*. The selected measurement times C* will generally be defined as those yielding extreme concentrations. For example, C^high will denote sampling periods yielding concentrations in the top quartile of an individual quarter in an individual year. (Extreme values are determined quarter by quarter and year by year in order to compensate for possible emission trends.)

Two distinct quantities will be derived from the time stamp counts: densities and conditional frequencies. A time stamp density, g_ij, is a count G_ij per cell area area(i,j), normalized by the total number of parcels. It has dimensions of stamps/(km²·parcel), or hour/km² since stamps have dimensions of hours·parcels. The product ( g [hours/km²] ´ q[kg/hours·km²] ´ D A[km²] ) of g with emissions density q and area element D A has dimensions kilogram per square kilometer of total vertical contaminant burden, and represents the contribution to the burden at the receptor from emissions in D A. A conditional frequency, CF_ij*, is the ratio of a special-case density g_ij* to the unrestricted density g_ij, and is accordingly dimensionless. It represents the proportion of time stamps contributed by selected parcels.

Our treatment departs slightly from the traditional one, which is presented in terms of probabilities and normalized distributions. Our conditional frequency can be interpreted as a conditional probability, but this requires additional and unnecessary assumptions. More important, it complicates the subsequent discussion of sampling uncertainty and statistical significance, which will require its own probabilistic framework. Retaining units in the time stamp densities helps to convey their physical significance, and is essential if results from different trajectory and grid formats are to be compared (e.g., 5-day versus 3-day trajectories, 12 trajectories/day versus 4 trajectories/day, 1-hour versus 3-hour segments, 1° ´ 1° versus 0.5° ´ 0.5° grid cells)

The special case time stamp density, g_ij*, is defined as follows:

The first factor in the numerator is the area density of time stamps on cell (i,j); the first factor in the denominator is the total number of parcels for which both trajectories and concentrations are available. Complete trajectories could not be calculated for some of the receptor measurement periods, and the fraction of all valid trajectories associated with days with concentration values in the selected quartile may thus differ from 25%. The rightmost terms in both the numerator and denominator are correction factors to remove the potential bias introduced by missing or incomplete trajectories; they are given by the ratio of the number of valid trajectories expected to arrive at the receptor during the sampling intervals (8-, 12-, 16- or 24-hour samples) to the actual number of trajectories successfully calculated,

and

These normalization factors assure that the number of special case time stamps are equal to approximately 25% of all time stamps.

Setting C* = C in 1 yields the unrestricted time stamp density g_ij, where all concentrations are considered. For this case

Figure 1. Time stamp density g_ij for all parcels in which methylchloroform was measured at Meadview in the summer. Results were obtained from 830 trajectories.

the correction terms f_{traj, selected} and f_traj,all cancel, and we have

(2)

The conditional frequency, CF_ij*, is defined as the ratio g_ij*/ g_ij, of time stamp densities,

Figure 2. Time stamp density g_ij^high for all parcels in which high methylchloroform concentrations were measured at Meadview in the summer. Results were obtained from 205 trajectories.

Figure 3. Time stamp density g_ij^rand for a random selection of parcels in which methylchloroform was measured at Meadview in the summer. Results were obtained from 205 trajectories.

Figure 4 shows the plot of CF_ij* obtained by taking the ratio, cell by cell, of results in Figures 2 and 1. The value at each grid cell is the conditional frequency, the proportion of parcels passing over the cell that subsequently arrive at Meadview with methylchloroform concentrations in the top

quartile. Figure 4 is "normalized" by the prevailing transport pattern, so this pattern no longer dominates the association between pathways and concentrations. Moreover, there is now a reference level against which empirical values can be compared: in the absence of any association between pathway and concentration, the expected conditional frequency would be 25% everywhere, since selected parcels make up 25% of the total. An example generated by random selection of parcels is given in Figure 5, which shows the conditional frequency CF_ij* corresponding to the density g_ij* exhibited in Figure 3.

The conditional frequencies in Figure 5 are not everywhere 25%, even though there is no underlying association between pathway and parcel selection in this case, but scatter to either side of this reference level. The observed scatter is strongest far from the receptor, and reflects the statistical variability inherent in small samples from a population. The issue of statistical significance must thus be addressed before conditional frequency maps can be reliably interpreted.

Figure 4. Conditional frequency CF_ij^high of parcels arriving at Meadview with high methylchloroform concentrations in the summer. CF values represent ratio of high-methylchloroform stamps (Figure 2) to all stamps (Figure 1) at each cell.

Figure 5. Conditional frequency CF_ij^rand for random parcel selection in Figure 3. High values of CF^rand arise by chance only.

Figure 5 shows that observed conditional frequencies can be expected to exceed the 25% expectation in some areas purely by chance. The magnitudes of such fluctuations vary with location, being least where high time stamp densities yield good counting statistics. Such spatial variations in precision complicate our interpretation of Figure 4 because they raise the possibility that the highest conditional frequencies, which are observed in the Pacific Ocean, central Nevada, and the Four Corners region, could be statistical accidents with no physical meaning. It is necessary to evaluate conditional frequency maps against the null hypothesis that there is no association between the concentrations and pathways of arriving air.

Without making any assumption about the statistical nature of the distribution of time stamps over the grid, one can use resampling techniques (alternatively known as "bootstrap" or "randomization" methods) to examine the statistical significance of conditional frequency values. This approach, straightforward but computationally intensive, consists of repeatedly computing CF_ij* for sample subsets containing 25%´ n(AÇ VÇ C) trajectories drawn randomly (without replacement) from the set of all trajectories. After a sufficient number (~1000) of repetitions, the "random" conditional frequency’s mean approaches 25% everywhere and their standard deviation s _ij stabilizes. It can then be estimated that 95% of the CF_ij* arising from the null hypothesis of no association

Figure 6. Conditional frequency significance thresholds estimated by random resampling; values plotted correspond to .

Figure 7. Conditional frequency significance thresholds calculated from binomial distribution; values plotted correspond to the least empirical CF_ij whose 95% confidence interval lies above 25%.

will lie within 2s _ij of 25%. The value 25%+2s _ij thus becomes a threshold below which CF_ij* is statistically not significant at the 95% (two-sided) confidence level. Figure 6 plots this 25%+2s _ij threshold for summer at Meadview. The surface in Figure 6 is bowl-shaped, with values that rise from near 30% along frequently traveled parcel pathways to near 100% at the periphery. An observation of CF_ij^high, when methylchloroform is used to select trajectories, needs to exceed the values plotted in order to be considered significant.

An alternative approach [Bresch et al., 1987] is to treat the stamp count G_ij*= G_ij(AÇ VÇ C*) as the number of "successes" in G_ij= G_ij(AÇ VÇ C) Bernoulli trials each with probability CF_ij*. This model ignores the interdependence of stamps left by the same trajectory in different grid cells. The confidence interval for G_ij*/G_ij as an estimate of CF_ij* can then be derived from the binomial distribution [IMSL/STAT Library, 1989]. An empirical value is considered statistically significant if the null level n(AÇ VÇ C*) / n(AÇ VÇ C) is below the 95% (two-sided) confidence interval for G_ij*/G_ij. To each G_ij corresponds a minimum significant value of G_ij*, and this in turn yields a threshold below which CF_ij* is not significant. Figure 7 plots this binomial test threshold for summer at Meadview. The surface in Figure 7 resembles that in Figure 6, with values a little lower in the middle and a little higher at the periphery.

The binomial test seems to be the more conservative in the most uncertain region, far from the monitoring site. If a conditional frequency estimate is significantly higher than the random expectation according to the resampling test, it is also likely to be significant by the binomial test. Because the binomial test is so commonly used by others, and is computationally easier, it will be adopted in the remainder of this paper despite the resampling test’s sounder theoretical foundation. It should be noted that the assumed independence of different time stamps becomes harder to justify if trajectory segments are interpolated to shorter intervals [Green and Gebhart, 1994].

Figure 8 distinguishes the grid cells whose CF_ij* values in Figure 4 are above and below the thresholds in Figure 7 for statistical significance. The dark region extending from Meadview west-southwest to Los Angeles and Point Conception is empirically associated with high methylchloroform concentrations at Meadview, to a degree not explained by chance. This result contrasts with that in Figure 9, the analo-

Figure 8. Significance of conditional frequencies CF_ij^high (Figure 4) for parcels arriving at Meadview with high methylchloroform concentrations in the summer.

gous plot for the randomly generated CF_ij* values in Figure 5. Many fewer of the randomly generated values are significantly high, and these are scattered about the map.

Figure 10 provides an alternative check on physical meaningfulness, showing grid cells significantly associated with low methylchloroform concentrations at Meadview, those in the bottom quartile. Even though constrained by the same prevailing southwesterly flow (Figure 1) from which the high-concentration associations in Figure 8 were derived, the area associated with low concentrations is wholly distinct, extending to the south-southwest.

Because all trajectories converge toward the monitoring site, certain grid cells may exhibit high conditional frequencies merely because they are located immediately upwind or downwind of genuine source cells. This coupling places limits on the spatial resolution that can be achieved by any analysis of back trajectories. The nature of these limits can be studied by applying the methods developed above to synthetic data generated by known, ideal sources.

Figure 11 shows the conditional frequency map corresponding to an imaginary source located in the 0.5° ´ 0.5° cell (a,b) centered at 34° N and 118° W (Los Angeles). This map was derived in the same manner as the methylchloroform map in Figure 4, but with a different parcel selection criterion. The selected parcels in Figure 11 are those which arrived at Meadview during a valid determination of methylchloroform and

Figure 9. Significance of conditional frequencies CF_ij^rand (Figure 4) for random parcel selection in Figure 3.

Figure 10. Significance of conditional frequencies CF_ij^low for parcels arriving at Meadview with low methylchloroform concentrations in the summer.

left at least one time stamp on cell (a,b): C^ab = {t(t ): G_ab(t ) > 0}. The conditional frequency is thus the proportion CF_ij^ab = G_ij(AÇ VÇ C^ab) / G_ij(AÇ VÇ C) of parcels passing over cell i,j that also passed over cell a,b. This is the CF map that would be generated by an ideal tracer with no background level and no sources outside cell a,b.

Some 148 trajectories out of 830 arriving at Meadview during a valid determination of methylchloroform over five summers passed over cell (a,b). This proportion, n(AÇ VÇ C^ab) / n(AÇ VÇ C) = 17.8%, constitutes the standard against which observed conditional frequencies should be judged. Figure 12, which distinguishes the conditional frequency values in Figure 11 that are significantly higher than this standard, shows the limited resolution of conditional frequency analysis, particularly in the radial direction since trajectories passing through cell (a,b) are likely to pass through its neighbors. In this example, a source confined to a single grid cell is smeared throughout a radial sector by the data inversion process. Our imaginary source could be more closely located by choosing higher reference levels for the observed conditional frequency, as Figure 11 shows, but any choice other than n(AÇ VÇ C^ab) / n(AÇ VÇ C) would be ad hoc and arbitrary.

A statistical sampling framework has been developed for analyzing spatial distributions of back trajectory segment endpoints as areal time densities; exposure hours per air parcel

Figure 11. Conditional frequency CF_ij^ab of parcels arriving at Meadview after passing over a cell centered at a = 34° N and b = 118° W (Los Angeles) in the summer. Results from 148 trajectories are depicted.

Figure 12. Significance of conditional frequencies CF_ij^ab (Figure 11) for parcels arriving at Meadview after passing over a cell centered at Los Angeles.

per square kilometer. Ratios of these densities for parcels arriving with high concentrations and all concentrations yield conditional frequencies; proportions of parcel exposures yielding high measured concentrations. High conditional frequencies indicate empirical associations between exposures to specific locations and subsequent observations of elevated concentrations. Such associations may arise from the presence of major emissions sources, and analyses of conditional frequency maps may thus shed light on possible locations of such sources.

Because empirical associations arise by chance as well as from transport patterns, observed conditional frequencies must be tested for their statistical significance. An objective rationale has been developed for such screening, and the idealized Bernoulli model for statistical fluctuations in conditional frequency has been validated against an exact model that accounts for spatial correlations. Limits on spatial resolution imposed by spatial correlations have been illustrated by applications to synthesized data from hypothetical point source emissions.

The method has been evaluated with data on concentrations near the Grand Canyon of methylchloroform, a major industrial solvent, and correctly identifies the Los Angeles basin as a candidate area of emissions. As might be expected, the directional resolution of conditional frequency maps is considerably better than the radial resolution. Future papers will treat spatial resolution as a function of distance, direction, and season, and apply conditional frequency analyses to particle light scattering, particle sulfur, and other aerosol fractions. The method examined here is universal, and could be applied with appropriate modifications to any set of back trajectories. Previous applications have employed the primitive but widely used ATAD model, and that tradition is continued here for consistency. Our results support the meaningfulness of ATAD trajectories in statistical aggregations, but do not address the meaningfulness of individual trajectories.

Acknowledgments. The authors would like to thank the Electric Power Research Institute (EPRI) for its support of this work. In particular, we wish to acknowledge Peter K. Mueller and Pradeep Saxena for their comments and suggestions. We would also like to thank Rudolph Husar and Bret Schichtel at the Center for Air Pollution Impact and Trend Analysis (CAPITA) for their useful discussions and suggestions. Aleksandar Juric and Kari Höijärvi at CAPITA provided essential support during use of the CAPITA developed software MapEdit.

E.S. Macias, L.A. de P. Vasconcelos, and W.H. White, Department of Chemistry, Box 1134, Washington University, St. Louis, MO 63130. (e-mail: macias@hilltop.wustl.edu; vasconcelos@wuchem.wustl.edu; white@ wuchem.wustl.edu)

J.D.W. Kahl and D. Liu, Department of Geosciences, University of Wisconsin, Milwaukee, P.O. Box 413, Milwaukee, WI 53201. (e-mail: kahl@csd.uwm.edu; desong@csd4.csd.uwm.edu)

(Received March 30, 1995; revised July 22, 1995;
accepted July 23, 1995.)