Environ. Sci. Technol. 1998, 32, 115-123
Chemical Dynamics of Persistent Organic Pollutants: A Sensitivity Analysis Relating Soil Concentration Levels to Atmospheric Emissions JOSEPH N. S. EISENBERG* School of Public Health, University of California, Berkeley, California 94720 DEBORAH H. BENNETT Department of Mechanical Engineering, University of California, Berkeley, California 94720 THOMAS E. MCKONE Lawrence Berkeley National Laboratory and School of Public Health, University of California, Berkeley, California 94720
A dynamic modeling framework is presented that provides a method to estimate deposition of air emissions in soil, based on either single or multiple time-point soil measurements. These dynamic estimates, made in the context of the uncertainty and variability associated with environmental systems, differ significantly from those based on steady-state single-valued assumptions. A case study on the global mass balance for polychlorinated dibenzo-p-dioxins (PCDD) and dibenzofurans (PCDF) illustrates that, by explicitly representing the uncertainty of the parameters used to assess mass balance, there exists a discrepancy between deposition and source emissions of between 6- and 20-fold. The uncertainty in these soil deposition estimates was most strongly dependent on two landscape properties (the deposition velocity of air particles and the thickness of the atmosphere) and three chemical properties (the organic-carbon partition coefficient, chemical degradation within the atmosphere, and vapor pressure). Using lower-bound values of the deposition velocity of air particles, lower-bound values of the organic-carbon coefficient, and upper-bound values of atmospheric thickness, this discrepancy was estimated to be over 1 order of magnitude; whereas using upperbound values of the deposition velocity combined with lowerbound values of vapor pressure, this discrepancy was estimated to be around 6-fold. Constraints on these five properties decreased the uncertainty measure, as expressed by the coefficient of variation, from 0.8 to as low as 0.17.
Introduction Accumulations of environmentally persistent chemical pollutants are increasing societal concerns about environmental sustainability and the health of ecosystems and humans (1* Corresponding author fax: 510-642-5815; e-mail address:
[email protected]. S0013-936X(97)00337-4 CCC: $14.00 Published on Web 01/01/1998
1997 American Chemical Society
5). These contaminants originate from human activities, such as combustion related to energy production, industrial processes, and agriculture, and accumulate in various environmental compartments, e.g., soils, vegetation, animals, and humans. Of particular concern are the persistent pollutants that may have reproductive and endocrine disrupting effects. Examples of persistent pollutants include dioxins, polychlorinated biphenyls (PCBs), pesticides (e.g., toxaphene, carbaryl, lindane), certain heavy metals (e.g., cadmium, lead, and mercury), and phthalates (6). Many of these chemicals share similar physical-chemical properties (7). They are semivolatile organic compounds (SVOCs) with vapor pressures typically between about 0.1 and 10-5 Pa; they are resistant to degradation; and they are lipophilic, partitioning into carbon and lipid phases of the biosphere. As a result, many of these compounds can persist for years, and thus have been labeled “persistent organic pollutants” (POPs) (8). There is increasing recognition that many of these air pollutants can be dispersed over large regions well beyond the communities and even the countries from which the sources first emanate (4, 9, 10). This recognition has led many to attempt to determine the global distribution of organic pollutants and how they are transported from their sources to various environmental media (7, 9, 11-13), and have led others to assess the regional and global mass balances of these chemicals, including dioxins (14-17). While much has been gained from research on the global distribution of POPs, many of these studies did not account for the dynamics of transport or for the uncertainty of these processes when analyzing data. Much of the data collected from the environment are taken at one point in time providing no temporal information and have been interpreted as information obtained from a steady-state system. Other data taken from sediment cores are also interpreted as if the system dynamics are instantaneous, again resulting in a steady-state analysis. Given that releases of many of these chemicals began in the 1930s (1-5), have rate processes on the order of 10-100 years, and have been emitted from sources at variable rates over time, these environmental systems are most likely in a transient rather than steady-state phase. The kinetics of these systems depend on both landscape and chemical properties. The uncertainty and variability of landscape and chemical properties place limits on the precision of chemical distribution estimates. Without the explicit treatment of these uncertainties, results from analytical studies produce point-value estimates, as opposed to distributions reflecting the full range of possible values. Furthermore, depending on the complexity of the model, these point values may or may not accurately represent central values. To adequately incorporate the dynamic and uncertain elements of the system, there is a need for models and field sampling programs that provide an appropriate regional/ global scale framework for assessing the dispersion, persistence, and potential long-term impacts on human health and ecosystems. Also, models and field data should be integrated to best inform control strategies at different levels of political hierarchyslocal, national, regional, and global. In this paper, we illustrate how simple dynamic models combined with uncertainty and sensitivity analyses can be used to put bounding interpretations on limited data. We also introduce a scheme for making this process iterative by systematically selecting parameters for which lack of information (high variance) is significant to outcome predictions. VOL. 32, NO. 1, 1998 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
115
This makes possible a “value-of-information” strategy in the evolution of global models. Specifically, to address these issues we have developed a simple but dynamic modeling framework from which both cumulative soil content measurements and multiple time-point measurements, such as sediment core measurements, can be interpreted and a simulation approach which explicitly acknowledges the uncertainty and variability associated with system parameters. A regional sensitivity analysis (18) is used to assess the relative contribution of each parameter on the spread of the output distribution. To illustrate this approach, we conduct a case study on the global mass balance for polychlorinated dibenzo-p-dioxins (PCDD) and dibenzofurans (PCDF) (throughout the remainder of the text, dioxin will be used to refer to PCDD/PCDF).
Mass Balance Before dioxin emissions can be managed effectively, the sources need to be identified and quantified. If all major sources cannot be identified, we have little confidence that managing known sources will reduce potential human and ecosystem exposure, especially if the largest sources are not yet known. To assess whether or not sources have been well characterized, a mass balance study can match these sources against global deposition patterns. In recent years, researchers have tried various approaches to match soil deposition with estimated emissions into the atmosphere. One approach is to make direct measurements of the deposition from the atmosphere to the soil (19). These measurements include not only deposition from present day sources but also that resulting from resuspension of old sources, which might account for up to 95% of the total deposition (20). Alternatively, measured soil concentrations have been used to estimate emissions (14, 21). There are a variety of approaches to interpret the concentrations in the soil in relation to emission estimates, all rely on either an implicit or explicit model. For example, Brzuzy and Hites (14) interpreted deposition data as a static system, with the soil as a conservative matrix. This result bounds the lower end of emissions because it assumes there has been no degradation in air or soil and that the system is in steady state. At the other extreme, the emissions are bounded by the compartment with the fastest dynamics, corresponding to all the chemical being in the air (22). A dynamic model provides a more flexible representation of the range of possible source emission levels by acknowledging that chemical distributions do not attain steady state instantaneously. The temporal characteristics must be accounted for in any system that, within the time frame of interest, has continual input, finite decay rates, and noninstantaneous transition rates between media. Our premise here is that a dynamic representation is necessary to evaluate the source/deposition ratios for dioxin. In defining the dynamic model, a two-compartment system, atmosphere and soil, was used to represent the global dioxin transport system. Two primary issues justify the use of a two-compartment model. First, although a third ocean compartment is usually included in a global distribution model, this compartment can be incorporated into the atmospheric half-life parameter since (a) the ocean is a sink and relatively unimportant as a re-emission compartment and (b) very little chemical will partition into the ocean due to the lipophilic nature of the chemical. The deposition to oceans is uncertain and has been estimated by others to be on the order of 5% of total global deposition (14). Second, the uncertainty and variability of the limited available information on the spatial relationship between sources and 116
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 32, NO. 1, 1998
sinks precludes the identification of more than a few aggregate parameters.
Approach and Methods The two-compartment model, shown in Figure 1, describes the time variation of the quantity of chemical in the atmosphere, Na(t), and in the soil, Ns(t). The movement of the chemical between the atmosphere and the soil is governed by the rate constants uas and usa, and the degradation is governed by the parameters µa and µs. There is one variable, S, that describes the input of source emissions into the atmosphere. All rate constants are presented in units of day-1, and the two state variables are presented in an arbitrary unit, which we call mass units (mu). Subsequently, the input variable is in units of mu/day. The equations for this model are shown in the Appendix (see Supporting Information). To provide a framework to assess the implications of either explicit or implicit assumptions, we present closed form solutions under five different sets of limiting assumptions of transport and degradation rates. These assumptions are as follows: (a) chemical transport is unidirectional, from atmosphere to soil, and there is no soil degradation; (b) condition a is relaxed to allow for soil degradation; (c) condition a is relaxed to allow for chemical transport from the soil to the atmosphere, but with no degradation in either the atmosphere or soil; (d) condition c is relaxed to allow for degradation in the air; and (e) condition d is relaxed to allow for degradation in both the air and soil. For each of these five conditions, both the early transient phase and steady-state phase are presented in the Appendix (see Supporting Information). These analytical solutions provided guidance for our numerical simulation studies. The results from the numerical simulation studies were used to interpret data from cumulative soil samples. Running simulations required parametric estimates of the two transport rate constants, which in turn required values for six chemical properties and 13 landscape properties. The equations for these rate constants are shown in Figure 1 and taken from CalTOX (23). These equations are based on diffusion as well as two advective processes: rain and dry particle exchange between soil and air. For example, the first equation shows that air to soil transport is calculated as the sum of three terms: dry deposition of the vapor phase, wet deposition of the vapor phase, and wet and dry deposition of the particulate phase. For further details, see McKone et al. (23). The values for the defined chemical and landscape properties are shown in Table 1. Included are nominal values for the chemical properties that are representative of 2,3,7,8-TCCD from U.S. EPA estimates (24), and for the landscape properties, reasonable average global values are given. In addition to these nominal values, uniform distributions were assigned to each property that account for the estimated uncertainty and variability of each value. Uniform distributions were used for a priori estimates so as to not bias one study over another; i.e., given that the available data are inadequate to statistically identify a given distribution, the uniform is considered the least biased estimator (25). Also listed in Table 1 are the degradation parameters in the atmosphere and soil. The atmospheric degradation halflife is assumed to be between 7 and 46 days, which takes into account both the more rapid vapor phase and the slower particle phase transformation (26-29). The soil degradation half-life is assumed to be 70 years. It should be noted here that we distinguish between chemical transformation and volatilization; volatilization is a loss term taken into account by the transport rate usa, whereas degradation refers only to chemical transformation. Chemical transformation is variable depending on the soil type, chemical form, depth of contamination, and co-contamination (30-36). Since most
FIGURE 1. Compartmental diagram of the model used to interpret deposition data. Na(t) and Ns(t) are quantities of chemical in the atmosphere and soil, respectively (mu); uas and usa are transfer rate constants from air to soil and soil to air, respectively; µa and µs are degradation rates in the atmosphere and soil, respectively (day-1); S is the input variable (mu/day); Yas is the fugacity mass-transfer coefficient (mol/Pa-m2); Da is the diffusion coefficient (m2/d) for ground surface soil; the Z values are the fugacity capacities for the different media; and da is the diffusion length in soil. Further definitions of parameters are shown in Table 1.
studies on chemical half-life in soil do not distinguish between these mechanisms, we selected 70 years for most simulations as an upper-end estimate from the literature. In addition, we ran one simulation using 700 years to assess the sensitivity of this parameter. Monte Carlo simulations were conducted using Crystal Ball (37) to generate the model parameter values and Matlab (38) to obtain soil concentration values over the simulation time period. For each simulation, the transport rates were calculated using eq A16 (see Supporting Information) and used a randomly selected set of parameter values from Table 1. (Four chemical parameters (Dair, Dwtr, H, and Tm) were assigned constant values since preliminary simulations suggested that the uncertainty of these factors did not appreciably affect the output.) These rates were then used in the dynamic model to estimate current dioxin levels within the soil, assuming that emissions of dioxins began 60 years ago (1-5) and that emission rates have been constant throughout this period. A subsequent set of simulations was conducted in which the assumption of constant emissions was relaxed. For these simulations, source emissions were assumed to be zero at time zero, increasing to 1 mu/day by year 40, and remaining at 1 mu/day for the next 20 years. This input profile corresponds to a linearly increasing source emission from 1935 to 1975, leveling to a constant value after 1975.
Using the soil concentration values, we estimated the following percentage, F:
F)
s × 1000 u
where s is the average net soil deposition over the 60-year period (mu/d), and u is the average source emission (mu/d). Therefore, F is an output variable that represents a measure of comparison between the source emission and a cumulative soil measurement; i.e., F is the percentage of source emissions that can be accounted for by a cumulative soil sample. The implications of a given F value are as follows: If, after 60 years, 50% of what was emitted ultimately accumulates in soil, the global source term required to account for the measured soil levels is 1/F or twice the amount in soil. If 10% is in soil, the source term must be 10 times what was in soil. Our Monte Carlo simulations produce a distribution for F, characterizing the range of possible values due to parametric uncertainty. Each single value of F is associated with a specific set of chemical and landscape properties. To analyze the output distribution of F and how this corresponds with the multivariate distribution of chemical and landscape properties, we use a combination of a Classification and Regression Tree (CART) algorithm and more traditional VOL. 32, NO. 1, 1998 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
117
TABLE 1. List of Chemical and Landscape Properties Used to Calculate Soil Concentrations nominal octanol-water partition coefficient vapor pressure (Pa) diffusion coeff; pure air (m2/d) diffusion coeff; pure water (m2/d) melting point (K) Henry’s law constant (Pa-m3/mol)
organic carbon partition coeff surface area of earth (m2) annual av precipitation (m/d) atmospheric dust load (kg/m3) deposition velocity of air particles (m/d) thickness of the ground soil layer (m) soil particle density (kg/m3) water content in surface soil (vol fraction) air content in surface soil (vol fraction) fraction of land area in surface water ambient environmental temperature (K) organic carbon fraction in upper soil zone boundary layer thickness in air above soil (m) atmospheric mixing height (m) reaction half-life in air (d) ) 0.693/µa reaction half-life in surface soil (y) ) 0.693/µs
Chemical Properties Kow 4.62 × 106 VP 1.6 × 10-6 Dair 0.421 Dwtr 5.11 × 10-5 Tm 578 H 2.47
Landscape Property A 5.1 × 1014 rain 2.4 × 10-3 Fba 5.0 × 10-8 vd 5.0 × 102 ds 1.0 × 10-2 Fss 2.6 × 103 βs 0.22 Rs 0.35 farw 0.7 T 288 focs 0.015 δas 5.0 × 10-3 da 1000 Degradation Parameters τa τs 70
Results The closed form solutions under the five conditions listed in the Methods section are shown in Table A1 (see Supporting Information). These equations suggest that assumptions that soil concentration levels track the source emissions with little or no delay are met in two of the five cases: the unidirectional case in which only atmospheric degradation is considered and uas . ma and the bi-directional case in which no 9
refs
[106, 109] [10-12, 10-5]
24, 42, 43 24 24 24 24 24, 42
[0.3 × 10-3, 3.0 × 10-3] [2.0 × 10-8, 10 × 10-8] [1.0 × 102, 1.0 × 103] [0.3 × 10-2, 3.0 × 10-2] [2.0 × 103, 3.0 × 103] [0.15, 0.25] [0.3, 0.4] [0.6, 0.8] [273, 301] [0-0.05] [1 × 10-3, 9 × 10-3] [300, 3000]
44 45, 46 47, 20, 19,46
Calculated Chemical Property Kds ) focsKoc Koc ) 0.48Kow
multivariate regression analysis (39). The goal in this analysis is to characterize the importance of variance in chemical and landscape properties, identifying the regions in this property space that result in the greatest reduction of estimated uncertainty for F as well as the regions that correspond to large and small values of F. The CART analysis produces a tree structure using a parametric decision at each node based on an inequality, pi < c, where pi corresponds to property i and c is a constant value. Observations that satisfy this condition are sent to the left node while the others are sent to the right node. Each node is characterized by a mean value and standard deviation associated with the observations satisfying the above inequality. A specific split is chosen that maximizes the reduction in variance between the parent node and the resultant two child nodes. Splitting continues until a stopping rule is satisfied: either the node has no variance or the number of observations is small. Higher nodes are revisited and if necessary readjusted to improve the lower level errors. This revisitation is especially important when particularly strong interactions among parameters exist. From this decision tree, an exhaustive list of subtrees is made. For each tree size, defined by the number of terminal nodes, the tree is chosen that maximizes the decrease in variance. In general, as a tree size increases, the variance is reduced, and as the tree size decreases, the variance is increased. The objective is to find an optimal combination of tree size and variance, somewhat analogous to choosing the correct bin size for a histogram (see Breiman (39) for details).
118
range
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 32, NO. 1, 1998
7-46
44 44 44 44 47
24, 26-29 24, 29-33
degradation is considered and uas . usa. Otherwise, soil concentrations are always going to be a fraction of the amount emitted by the source; i.e., F < 100. To ascertain the dynamic properties of this system under nonlimiting conditions, we conducted numerical simulations. Five thousand simulations were run. The transport rate constants varied over many orders of magnitude: uas ranged over 3 orders of magnitude, from 10-3 to 1.0 d-1; whereas usa ranged 6 orders of magnitude, from 10-8 to 10-2 d-1. However, for the majority (90%) of simulations, usa ranged over only 4 orders of magnitude, between 10-7 and 10-3. Figure 2 is a scatter plot showing the transport rate parameter values for each of the 5000 simulations. For reference, four line types are marked on the graph to distinguish regimes with different F values and attributes. The upper horizontal line designates a 60-year time period, τ. All points either near or below this line are still in the transient phase of the simulation after 60 years. The two vertical lines designate the upper and lower bounds for atmospheric degradation. Points near or to the left of these lines are parameter sets that represent appreciable degradation within the atmosphere. The lower horizontal line designates the soil degradation corresponding to a 70-year half-life. Above the line are parameter sets that represent an appreciable amount of soil degradation. The diagonal lines designate the ratio of rate constants, usa/uas. The closer usa is to uas, the greater the relative residence time within the air compartment as compared with the soil compartment and, therefore, the more air degradation that occurs. Thus, moving toward the upper left-hand corner results in higher overall degradation rates. Since µa is orders of magnitude higher than µs, overall degradation is very sensitive to this ratio shift. To illustrate the implications of the parameter space data points in Figure 2, Figure 3A shows time trajectories for two example parameter sets. The dashed line is a trajectory using transfer rate values from region 3, a typical point in the lower right-hand region that is characterized by fast transfer from air to soil and slow transfer from soil to air. This results in
FIGURE 2. Scatter plot of the soil to atmosphere transition rate parameter (usa) vs the atmosphere to soil transition rate parameter (uas) for the 5000 Monte Carlo simulations. Various line types and numbered regions are drawn for reference (see text). a short residence time in air and thus little decay in air, allowing a continuous buildup in the soil. The solid line trajectory is from region 9 that quickly equilibrates and is characterized by fast transfer from soil to air. The concentration in the air remains relatively constant, quickly coming to equilibrium with the soil as dioxin is degraded in air, driving the system toward steady state. Points residing between regions 3 and 9 exhibit intermediate transfer rate values, whereas points within regions to the left of regions 3, 6, and 9 are characterized by slower movement from the air to soil, resulting in an increased residence time in the air compartment and thus increased degradation. Therefore, each region is characterized by a trajectory with a similar shape as other regions with similar values of usa but attaining a lower level of soil concentration for smaller values of uas. The percentage of the emitted source that resides in the soil after 60 years, F, was calculated for all 5000 data points. For the constant emission case, the histogram of percentages can be seen in Figure 3B. The mean value for the 5000 simulations was 22% with a standard deviation of 18% and a range of 0-73%. An additional simulation was run that varied the source term, assuming zero emission at time equal zero, increasing linearly to a unit mass per day after 40 years, and maintaining that source term for the next 20 years. This input profile corresponds to a linearly increasing source emission from 1935 to 1975, leveling to a constant value after 1975. The ratio of F resulting from constant emissions over F resulting from non-constant emissions was 89 ( 0.07%, suggesting that, though there was a quantitative difference shifting toward larger percentage values, this difference was not large. Results from the CART analysis are summarized in Figure 4. The numbers within the oval are mean values, and the numbers below are standard deviations. Therefore, the top node shows that the overall mean of F for the 5000 simulations was 22%. Lower level nodes show means and standard deviations of regions defined by the splitting rules. For example, node 7 represented the simulations in which the factors vd (the deposition velocity), da (the atmospheric thickness), and koc (the organic-carbon partition coefficient) were constrained to the ranges [10, 103] m/d, [714, 3000] m, and [1.0e5, 3.2e8], respectively; this subregion contains a mean F of 6 ( 6%. The CART analysis suggests that the landscape properties most sensitive to the classification of
FIGURE 3. Simulations of chemical accumulation into the soil given a constant emission into the atmosphere. (A) Trajectory for two parameterizations, one leading to a steady-state system while the other leads to a transient state: the dashed line, uas ) 0.2 d-1, usa ) 2e-6 d-1 represents relatively rapid transfer from air to soil and slow transfer from soil to air (region 3); the solid line, uas ) 0.2 d-1, usa ) 2e-4 d-1, represents more rapid transfer from soil to air (region 8). For both plots ma ) 0.04 (half-life ) 17.25 days) and ms ) 2.74e-5 (half-life ) 70 years). (B) Density plot of all 5000 G values, the ratio of the average net soil deposition to the average source emission over a 60-year period, based on uncertainty ranges for landscape and chemical properties. output values with a minimal variance were vd (deposition velocity) and secondarily da (effective atmospheric mixing height). The three chemical properties, koc, VP (vapor pressure), and µa (atmospheric degradation rate) were also found to be sensitive parameters. This classification had an equivalent R 2 of 0.73, suggesting that 73% of the output data could be explained by the decision tree shown in Figure 4. To verify this result, we looked at various combinations of these three parameters, with and without a log transformation. We found that the best fit, using multiple regression analysis, was a linear model of the following form:
F ) -232µa - 0.009da + 15 log(koc) + 7 log(VP) 16 log(νd) - 4 log(νd) log(VP) (R 2 ) 0.78). Table 2 shows the statistics associated with this fit. The only significant interaction was between vd and VP. The results of the CART analysis suggest that the factor vd was the single-most sensitive parameter in determining whether the output value of F was high or low; i.e., for vd < 103 m/d, all nodes had mean values of less than 20, with the exception of node 4; for vd > 103 m/d, all nodes had mean VOL. 32, NO. 1, 1998 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
119
FIGURE 4. CART description of a multiple regression in which the 12 landscape factors were the independent variables and G was the dependent variable. Inside the ovals are mean G values, and below each oval are standard deviations for the subregion defined by the branching rules.
TABLE 2. Statistics for Multivariate Regression of G (p < 0.0001) intercept µa da log(vd) log(VP) log(koc) log(vd) × log(VP)
coeff
standard error
t value
52 -232 -0.009 -15 7.0 2.0 -4
2.1 5.0 1.6e-4 0.9 0.23 0.055 0.11
25 -47 -56 -16 30.0 37 -37
values of greater than 22 with the exception of node 24. In general, a small vd (714 m), and a small koc (307 m/d) and either a small VP (
