Chemosphere 66 (2007) 1396–1407 www.elsevier.com/locate/chemosphere
A multimedia environmental model of chemical distribution: Fate, transport, and uncertainty analysis Yuzhou Luo, Xiusheng Yang
*
Department of Natural Resources Management and Engineering, University of Connecticut, Storrs, CT 06269, USA Received 17 April 2006; received in revised form 12 September 2006; accepted 17 September 2006 Available online 13 November 2006
Abstract This paper presented a framework for analysis of chemical concentration in the environment and evaluation of variance propagation within the model. This framework was illustrated through a case study of selected organic compounds of benzo[a]pyrene (BAP) and hexachlorobenzene (HCB) in the Great Lakes region. A multimedia environmental fate model was applied to perform stochastic simulations of chemical concentrations in various media. Both uncertainty in chemical properties and variability in hydrometeorological parameters were included in the Monte Carlo simulation, resulting in a distribution of concentrations in each medium. Parameters of compartmental dimensions, densities, emissions, and background concentrations were assumed to be constant in this study. The predicted concentrations in air, surface water and sediment were compared to reported data for validation purpose. Based on rank correlations, a sensitivity analysis was conducted to determine the influence of individual input parameters on the output variance for concentration in each environmental medium and for the basin-wide total mass inventory. Results of model validation indicated that the model predictions were in reasonable agreement with spatial distribution patterns, among the five lake basins, of reported data in the literature. For the chemical and environmental parameters given in this study, parameters associated to air–ground partitioning (such as moisture in surface soil, vapor pressure, and deposition velocity) and chemical distribution in soil solid (such as organic carbon partition coefficient and organic carbon content in root-zone soil) were targeted to reduce the uncertainty in basin-wide mass inventory. This results of sensitivity analysis in this study also indicated that the model sensitivity to an input parameter might be affected by the magnitudes of input parameters defined by the parameter settings in the simulation scenario. Therefore, uncertainty and sensitivity analyses for environmental fate models was suggested to be conducted after the model output was validated based on an appropriate input parameter settings. 2006 Elsevier Ltd. All rights reserved. Keywords: Great Lakes; Monte Carlo method; Multimedia environmental model; Sensitivity analysis; Uncertainty analysis
1. Introduction Multimedia environmental fate models are increasingly recognized as useful tools to predict screening level chemical distributions in multiple environmental media simultaneously. These models facilitate management of chemicals by providing approximations of the actual and more complex fate of chemical transport and transformation within
*
Corresponding author. Tel.: +1 860 486 0135; fax: +1 860 486 5408. E-mail address:
[email protected] (X. Yang).
0045-6535/$ - see front matter 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.chemosphere.2006.09.026
and among media. Many of the quantitative parameter values used in the multimedia chemical fate assessments are either uncertain or variable or both. True uncertainty in parameters can result from measurement or estimation error. In contrast, variability refers to quantities that are distributed empirically (McKone, 1996a). For example, chemical properties may be highly uncertain due to difficulty in measurement, while landscape data are highly variable at spatial scales (Bennett et al., 1999). Therefore, it is very important to present clearly the propagation of parameter variances into model outputs that reflect environmental persistence of chemicals. A major objective of
Y. Luo, X. Yang / Chemosphere 66 (2007) 1396–1407
parameter uncertainty and sensitivity analyses is to identify individual model inputs that most significantly contribute to uncertainties in model outputs. These results provide guidances for additional research efforts to increase model accuracy, and can also be useful for model users in applying the existing models to other chemical species or landscape conditions. With the increasing computational power of personal computers, probability distributions are used in place of discrete values, and appropriate Monte Carlo analysis is currently the major technique for quantifying uncertainty in environmental assessments (McKone, 1996a; Bennett et al., 1999; Liu et al., 1999). Most of the existing probability analyses in the multimedia environmental fate models consider only the uncertainties in the chemical properties (Kuhne et al., 1997; Citra, 2004; Fenner et al., 2004). There is a small, but growing, number of multimedia environmental fate models that perform stochastic simulations by including both the uncertainty in chemical parameters and spatiotemporal variability within the environment (McKone, 1996b; Bennett et al., 1999; Liu et al., 1999; MacLeod et al., 2002). Few efforts, however, are found in the literature to the uncertainty and sensitivity analyses for a spatially explicit simulation of multimedia environmental fate. In such models, the meteorological and hydrological conditions associated with air or water connectivity, such as wind speed, wind direction, and stream flow rate, were proven to be very important parameters in determining the pattern and trend of chemical distributions (Cohen and Cooter, 2002a,b; Luo et al., in press-b). The goal of this paper was to extended Monte Carlo analysis into a fugacity-based multimedia environmental fate model by analyzing chemical transport and fate across multiple environmental media. This approach enabled holistic integration of modeling capabilities that have not been fully achieved in existing models, in terms of multimedia mass transport simulation, spatial resolution at basin scale, and uncertainty/sensitivity analyses. The propagation of the uncertainty and variability associated with physiochemical, geophysical and landscape parameters were demonstrated by performing stochastic calculation of the multimedia simulation models. In this study, the primary probability analyses included uncertainty analysis, sensitivity analysis, and comparison of input variances in terms of their contribution to the variances in outputs. The methodology was demonstrated through a case study involving two representative chemicals (benzo[a]pyrene and hexachlorobenzene) in the Great Lakes region. 2. Approach 2.1. Multimedia environmental fate model A multimedia environmental fate model has been developed at the University of Connecticut for analyzing unsteady-state dispersion and distribution of chemicals in multimedia environmental systems (Luo et al., in
1397
press-a). Chemical transport processes were formulated in seven environmental compartments of air, canopy, surface soil, root-zone soil, vadose-zone soil, surface water, and sediment. The model assumed that all compartments were completely mixed and chemical equilibrium was established instantaneously between the phases within each compartment. A fugacity approach was utilized to formulate the mechanisms of diffusion, advection, physical interfacial transport, and transformation reactions. The governing equations of chemical mass balances in the environmental compartments were solved simultaneously to reflect the interactions between the compartments. Temporal variations in model parameters were identified by considering hydrologic flows, temperature, canopy growth, chemical emissions, and distant inputs. The outputs of the model included time-dependent chemical concentrations in each compartment and its sub-compartments, and inter-media mass fluxes between adjacent compartments at daily time steps. Contaminants are moved among and lost from each compartment through a series of transport processes. These processes were represented mathematically by first order equations in the model based on the fugacity concept Dij Qij ¼ M ij N i ¼ ð1Þ Ni V iZi where Qij (mol s1) is the unidirectional chemical flux from compartment i to j, Mij (s1) is the transfer rate constant from i to j, Ni (mol) is the chemical inventory in i, Dij (mol Pa1 s1) is the overall transport coefficient (Mackay-type D value) from i to j, Vi (m3) is the total volume of the compartment, and Zi (mol Pa1 m3) is the fugacity of chemical in the compartment. Similar formulations were used for other dilution mechanisms such as degradation, sediment burial, and transport to external environment. Mass balance for the chemical inventory in a compartment was described by a set of differential equations ! m m X X dN i ¼ Si þ ðM ji N j Þ M ij þ M ix þ M Ri N i dt j¼1 j¼1 ¼ Si þ
m X
ðM ji N j Þ M Oi N i
ð2Þ
j¼1
where m is the total number of defined compartments, Si (mol s1) is a total source term as the sum of distant sources, local emission, and mass gain from transformations, MRi (s1) is the reaction rate constant of chemical degradation, Mix (s1) is the transfer rate constant of chemical loss to a hypothetical receiver x in external environment, and MOi (s1) is the total loss rate constant of the chemical in compartment i. For the whole simulation domain, Eq. (2) was written for m compartments to solve for the time-dependent chemical inventories. Each of these equations was rearranged in implicit forms and resolved by a finite difference method using MATLAB code in IBM PC platform.
1398
Y. Luo, X. Yang / Chemosphere 66 (2007) 1396–1407
Geographic information system (GIS) database and geospatial analysis were integrated into the chemical transport simulation to provide spatially explicit estimations of model parameters at watershed scale. With watershed segmentation, soil and canopy compartments were considered to be isolated from the neighborhood. Water connectivity was based on the channel network and the actual river flow direction. The connectivity of air compartments was estimated from wind direction in a grid system based on similar studies at University of Connecticut (Xu et al., 2000a,b; Luo et al., in press-b). The transport fluxes of chemicals in air were calculated by applying an advective-dispersion equation in the grid system. The inter-media transport fluxes of chemicals between air and underlying compartments (canopy, surface soil, and surface water) were computed based on the projective area of air columns on the ground. Using trichloroethylene (TCE) as a test agent, the model was used to explore the behavior of contaminant in the Connecticut River Basin (Luo and Yang, 2005; Luo et al., in press-b). Driven by actual weather conditions from 2000 to 2002, the model results were reported as time series of concentrations and inter-media transport fluxes in the environmental compartments. The predicted concentrations were compared to field measurements or predictions from well-validated models. The monthly variations of the TCE concentrations were evaluated by comparing the temporal trends in chemical inputs and inter-media transport fluxes. The results of comparison indicated that the model simulation yielded reasonable agreement with reported data in the literature. 2.2. Stochastic simulation The performance of stochastic multimedia environmental fate simulation was demonstrated by applying the
model into the field conditions of the Great Lakes region (Fig. 1). Benzo[a]pyrene (BAP) and hexachlorobenzene (HCB) were selected as chemical tracers in the stochastic simulation. Both the true uncertainty in the chemical properties and the spatiotemporal variability in the environmental parameters, such as landscape properties, meteorological condition, and hydrological flows, were considered in the stochastic simulation. To simplify the calculation, some model input parameters were assumed to be constant, including the environmental compartmental dimensions, densities, chemical emissions, and background concentrations. Model input parameters were assumed to follow independent lognormal distributions (Slob, 1994; McKone, 1996a,b; Liu et al., 1999; MacLeod et al., 2002). The lognormal distribution is parameterized by the mean (l) and standard deviation (r, or r2 as variance) of the corresponding normal distribution on a log scale. In the uncertainty analysis, variance in lognormal distribution is often expressed by a coefficient of variance (CV), defined as the ratio of the standard deviation to the mean of the original samples (not on a log scale) (McKone, 1993, 1996a; Liu et al., 1999; MacLeod et al., 2002). The theoretical relationship between the standard deviation of a lognormal distribution and the corresponding coefficient of variance is r2 ¼ lnðCV2 þ 1Þ
ð3Þ
In a practical random data generation, however, model parameters are sampled within a specific interval of cumulative frequency. Given CV value in random sample generation, the resultant standard deviation may be less than the theoretical values given by Eq. (3). For example, based on a Latin Hypercube Sampling process given CV = 1 and 95% of cumulative frequency, the actual standard deviation is 0.725 vs. calculated value of 0.832 from Eq. (3).
CANADA Ontario
Quebec N
Lake Superior W
E S
Wisconsin
Lake Huron
Minnesota
Lake Ontario
Lake Michigan
New York
Lake Erie
Illinois
Indiana
UNITED STATES
Ohio
Pennsylvania 200
0
200
400
km
Fig. 1. The Great Lakes region modeled in this study. The shaded area shows the drainage basins of the Great Lakes.
Y. Luo, X. Yang / Chemosphere 66 (2007) 1396–1407
Each model input parameter was characterized by its mean and CV value, while a zero CV indicted a model input with fixed values assumed in the model simulation. Model simulation was conducted to predict chemical concentrations in all compartments for a ten-year period from 1991 to 2000 at monthly steps. The chemical emissions prior to the simulation period were accounted in the background concentrations. The driving forces were average emission rates and 10-year means of monthly average meteorological and hydrological conditions. Twelve datasets of input parameter were generated for each month of a year, and used repeatedly for the whole simulation period. For each simulation year, the model outputs were reported as monthly time series. The model simulation achieved dynamic equilibrium between compartments in 2 to 3 simulation years. For the sequent simulation years, the model outputs were the same time series at monthly step and repeated year by year, due to the monthly variations in hydrological and meteorological conditions. Therefore, model outputs in the last simulation year (2000), which reflected the average chemical distribution pattern for the whole simulation period during 1991–2000, were used for model evaluation and uncertainty/sensitivity analyses. The means and CV values of predicted chemical concentrations during simulation year 2000 in compartments of air,
1399
surface water, and sediment were compared with reported data in the literature. Observed chemical concentrations of benzo[a]pyrene and hexachlorobenzene in the air during 1990–2000 were retrieved from the Integrated Atmospheric Deposition Network (IADN, 2006). Measured concentrations of chemicals in surface water and sediment were obtained from a large number of sources complied by Zhang et al. (2003) and Booty et al. (2005). The chemical portioning patterns were compared with the CHEMGL and Mackay Level III models results (Zhang et al., 2003), in which steady-state mass balance models were applied to the field condition of the Great Lakes regions, without considerations of parameter uncertainties. Part of the sensitivity results were compared with the study by MacLeod et al. (2002) in which an analytical approach was used for uncertainty and sensitivity analyses in the ChemCAN model with a case study of benzo[a]pyrene fate in Southern Ontario. 2.3. Uncertainty in input parameters The input parameters required for the model simulation comprised environmental properties, chemical properties, chemical releases, and background concentrations. Environmental properties included landscape parameters,
Table 1 Great Lakes landscape properties and flow rates as annual means and coefficients of variance (Mackay, 1989; Holtschlag and Nicholas, 1998; Scurlock et al., 2001; GLERL, 2006; GLIN, 2006) Assumed CVa
Parameter
Annual means of the parameter in the basin of Lake Superior
Lake Michigan
Lake Huron
Lake Erie
Lake Ontario
Total basin area (km2) Land area (km2) Vegetated area (km2) Water and sediment area (km2) Air boundary layer height (m) Equivalent foliage depth (mm) Depth of surface soil (m) Depth of root-zone soil (m) Depth of vadose-zone soil (m) Water depth (m) Sediment active depth (m) Wind speed (m/s) Wind direction ()b Air temperature over basin (C) Water temperature (C) Outflow rate (m3/s) Transpiration rate (m3/s) Phloem flow (m3/s) Precipitation to land (mm/year)c Evapotranspiration (mm/year)c Surface runoff (mm/year)c Infiltration (mm/year)c Precipitation to lake (mm/year)d Lake evaporation (mm/year)d Soil erosion (mm/year) Sediment burial (mm/year)
210 000 128 000 64 000 82 100 700 2.6 0.01 0.75 5.5 152 0.01 5.7 183 3.1 4.8 2250 1509 75 757 372 97 288 774 594 0.2 0.06
176 000 118 000 59 000 57 800 700 2.6 0.01 0.75 5.5 85 0.01 6.0 160 7.4 5.4 1550 1893 95 812 506 65 241 802 643 0.2 0.60
193 000 134 000 67 000 59 600 700 2.6 0.01 0.75 5.5 59 0.01 3.6 179 5.9 6.7 5710 1359 68 826 320 143 363 831 631 0.2 0.60
104 000 78 100 39 050 25 700 700 2.6 0.01 0.75 5.5 19 0.01 6.6 175 9.3 5.3 6660 1504 75 862 607 133 122 890 906 0.2 0.60
83 000 64 000 32 000 19 000 700 2.6 0.01 0.75 5.5 86 0.01 5.1 177 7.5 6.3 7930 893 45 915 440 166 309 854 661 0.2 0.60
a b c d
Zero CV indicates a constant value was assumed for the corresponding parameters in this study. Wind direction degree: 0 = from North, 90 = from East, etc. Flow rate expressed in basin land-area equivalent value. Flow rate expressed in lake-area equivalent value.
0 0 0 0 0.2 0.2 0 0 0 0 0 0.1 0.5 0.2 0.2 1.0 1.0 1.0 0.1 1.0 0.4 1.0 0.1 0.1 1.2 5.0
1400
Y. Luo, X. Yang / Chemosphere 66 (2007) 1396–1407
hydrological flow rates, and Meteorological data. Shown in Table 1 are dimension parameters and flow rates as annual averages for each of the five lake basins. The height of air boundary layer was estimated as 700 m for all basins (Hanna et al., 1982), while other parameters for basin dimensions were obtained from the spatial data supplied by the Great Lakes Information Network (GLIN), or followed suggested settings for the Great Lakes in the literature (Mackay, 1989; McKone, 1993; MacLeod et al., 2002; GLIN, 2006). The plant foliage compartment was characterized by 50% vegetation coverage and a leaf area index of 8.72 (Scurlock et al., 2001). Meteorological and hydrological data of wind, temperature, rainfall, runoff, evaporation, and outflow rates were obtained from the Great Lakes Environmental Research Laboratory in the National Oceanic and Atmospheric Administration (NOAA) (GLERL, 2006). There data are derived from daily or monthly station measurements, and reported monthly averages for over-land and over-lake datasets for each of the Great Lakes (Croley and Hunter, 2001). Evapotranspiration and subsurface runoff rates were calculated based on the results of the basin water supply for the Great Lakes (Holtschlag and Nicholas, 1998). Transpiration and phloem flow rates were assumed to be proportional to foliage area (Paterson et al., 1994). Landscape parameters of compartment composition, organic carbon content and density are presented in Table 2. To simplify the simulation, annual mean values of these parameters
Table 2 Other environmental properties as annual means and coefficients of variance of the Great Lakes drainage basins (Mackay, 1989; McKone, 1993; Zhang et al., 2003; Booty et al., 2005; USDA, 2006) Parameter
Annual means
Assumed CVa
Volumetric content Fraction of particles in air Fraction of particles in water Fraction of biota in water Fraction of air in surface soil Fraction of water in surface soil Fraction of air in root-zone soil Fraction of water in root-zone soil Fraction of air in vadose-zone soil Fraction of water in vadose-zone soil Fraction of pore water in sediment
2.0 · 1011 2.08 · 107 1.0 · 106 0.20 0.45 0.2 0.3 0.1 0.3 0.85
0.2 1 0.2 0.2 0.4 0.3 0.5 0.4 0.5 0.2
Organic carbon (OC) content in particles OC content in surface soil 0.02 OC content in root-zone soil 0.02 OC content in vadose-zone soil 0.001 OC content in water particles 0.08 OC content in sediment particles 0.03
2.0 2.0 0.14 0.5 1.0
Density (kg/m3) Plant fresh mass density Soil solids Suspended solids Aquatic biota Sediment solids
0.2 0 0 0 0
830 2400 2000 1000 2000
a Zero CV indicates a constant value was assumed for the corresponding parameters in this study.
were assumed to be invariant over the five basins. The aerosol loading was estimated from total suspended particles (TSP) in the Integrated Atmospheric Deposition Network (IADN), by assuming air particle density of 1500 kg/m3 (Mackay, 1989). Organic carbon content in soil layers were adopted in the literature or obtained from the State Soil Geographic (STATSGO) database (USDA, 2006). For volumetric contents and densities which were not available in the literature, the values were taken from the Equilibrium Criterion or EQC model (Mackay et al., 1996; Mackay et al., 2006). In addition to the annual or monthly mean values, coefficients of variance for input parameters are also shown in Tables 2 and 3. The CV values for hydrological and meteorological data were processed from the measurements in the NOAA Great Lakes Environmental Research Laboratory (GLERL, 2006). The variabilities in other properties were estimated from the uncertainty analysis studies in the literature (McKone, 1996a,b; Liu et al., 1999; MacLeod et al., 2002). To simplify the calculation, all five basins shared an average CV value for each input parameter. The parameters for compartmental dimensions, except for air upper boundary and plant foliage, were assumed fixed values and marked with zero CVs. Benzo[a]pyrene and hexachlorobenzene were selected due to their physical–chemical properties, and as important contaminants in the Great Lakes region. More over, the availability of environmental data allows validating their environmental multimedia fate simulation. These two chemicals are recognized by the Canada-Ontario Agreement (COA) as toxic, persistent, and bioaccumulative chemicals in the Great Lakes. There are great uncertainties associated with their chemical properties, especially the half lives for sediment. Table 3 shows the physicochemical properties and degradation rate constants with coefficients of variance for these two chemicals, which were obtained from the supporting database of CalTOX model (McKone and Enoch, 2002). The transformation rate in foliage or root was assumed to be the same as the growth rate due to the lack of available data (Trapp and Matthies, 1998). The coefficients of variance were set as 1.0 and 1.1 for the half-life in foliage and root, respectively. The annual mean emissions of benzo[a]pyrene and hexachlorobenzene are also shown in Table 3. Chemical air emissions were obtained from the Great Lakes Regional Air Toxic Emissions Inventory (GLC, 2004), while chemical discharge to soils were estimated from published data (Zhang et al., 2003; Booty et al., 2005). To highlight the outcome variance by chemical and landscape properties, uncertainty in emission data was not considered in the stochastic simulation in this study. Uncertainty in emission or discharge data can be very high and is often difficult to estimate (Maddalena et al., 2001). Studies with estimated emissions found that uncertainty in emission contributed less to outcome variance than chemical properties (Bennett et al., 1999; Hertwich et al., 1999; Huijbregts et al., 2000). Apart from direct emission, the chemical transport from adjacent
Y. Luo, X. Yang / Chemosphere 66 (2007) 1396–1407
1401
Table 3 The physiochemical property data and daily emission rates for benzene[a]pyrene and hexachlorobenzene in this study Parameter
Benzo[a]pyrene
Hexachlorobenzene
Air Surface soil Root-zone soil Vadose-zone soil Surface water Sediment
252.30 (0.01) 7.13 · 107 (0.07) 0.09 (1.00) 2.49 · 10+6 (0.91) 0.44 (0.08) 5.26 · 105 (0.25) 450.65 (0.03) 329.00 (0.41) 0.06 (1.00) 228.6 (1.10) 228.6 (1.20) 880.0 (1.00) 2.3 (1.20) 1172.5 (1.40)
284.79 (0.01) 1.73 · 103 (0.32) 107.40 (0.46) 45809 (0.56) 0.47 (0.08) 5.83 · 105 (0.25) 510.15 (0.03) 13490 (0.60) 844.7 (1.00) 466.2 (1.10) 466.2 (1.20) 5147.0 (1.00) 1533.0 (1.20) 4672.5 (1.40)
Lake Lake Lake Lake Lake Lake Lake Lake Lake Lake
21.0 40.0 28.0 38.0 40.0 N/A N/A N/A N/A N/A
10.5 11.8 0.15 10.7 0.25 130.0 130.0 0.56 130.0 0.56
a
Physiochemical property Molecular weight (g/mol) Vapor pressure (Pa) Henry’s law constant (Pa-m3/mol) Organic carbon partition coefficient (l/kg) Diffusivity in pure air (m2/d) Diffusivity in pure water (m2/d) Melting point (K) Bioconcentration factor () Reaction half lives (day)
Daily emission ratesb Emission to air (kg/day)
Emission to soil (kg/day)
a b
Superior basin Michigan basin Huron basin Erie basin Ontario basin Superior basin Michigan basin Huron basin Erie basin Ontario basin
Coefficients of variance in parentheses. Emission rates were assumed to be fixed values, i.e., the corresponding coefficients of variance were zero.
basins via air and water flows were recognized as important sources. Variations in these processes were evaluated by assigning CV values for wind speed, wind direction, and lake outflows (Table 1).
2.4. Uncertainty and sensitivity analysis Latin Hypercube Sampling (LHS) was conducted to generate random input datasets for each simulation run. In LHS, the data range is divided into equiprobable subintervals, and stratified random samples are generated for each of the subinterval. In this study, model input parameters were taken from the 95% of cumulative frequency of the corresponding lognormal distribution as defined in Tables 1–3. Propagation of variances in chemical properties and landscape parameters into a model output was evaluated by a Monte Carlo analysis. Uncertainty analysis in this study was aimed to estimate total variance associated with the model outcome, and to evaluate the variability of model predictions when some uncertainty parameters were fixed. The stochastic simulations of chemical fate and transport were based on the random input dataset generated by LHS. Firstly, the distributions of the predicted range of model outputs were determined by allowing all parameter values to vary simultaneously. Furthermore, the environmental property dataset and the chemical prop-
erty dataset were varied independently. The results of the stochastic simulation were reported as means and standard deviations of chemical concentration in each of the environmental compartments. For each basin, the moments of total chemical inventories were also reported. At least 5000 discrete Monte Carlo events were included in each simulation set. The number of repetitions was determined by requiring consistent values of the first and second moments of all output distributions, i.e., the moments were invariant from one simulation to another when a different random input dataset was used. A sensitivity analysis was performed to understand which input parameters were important to the output by characterizing the degree of monotonic relationship between the model prediction and uncertain inputs. In a conventional sensitivity analysis, the sensitivity (S I i ) of output to changes in input is expressed as oO I i ð4Þ jS I i j ¼ oI i O where Ii is the ith model input, and O is the model output. Eq. (4) implied that partial derivatives of input parameters were assumed to be constant over the range of the uncertainties. The sensitivity analysis could be considered as a mathematical process, in which the complicated relationship between model output and inputs are approximated by using a power function
1402
Y. Luo, X. Yang / Chemosphere 66 (2007) 1396–1407 S
O ¼ a I i Ii
ð5Þ
where a is a constant. For example, unit sensitivity (S = 1) indicates that a linear relationship, with zero interception, could be established between the model output and the ith input parameter. It’s noteworthy that the results of sensitivity analysis may be only valid for the given model input settings (i.e., the means and CV values) for chemical properties, environmental parameters, and emission rates. If the input parameter settings change significantly, the sensitivity results may also change by applying different power functions in approximating the relationship between output and inputs. When fitting with a log-normal distribution, standard deviation of the input or output distribution is a measure of the magnitude of variation in the input or output variables (MacLeod et al., 2002). Therefore, the sensitivity can be calculated as the ratio between the standard deviations of output and input S 2I i ¼
r2O;I i r2I i
ri ti ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð1 ri Þ=ðN 2Þ
ð9Þ
has a Student’s t-distribution in the null case of zero correlation. For 5000 random samples in this study, the critical value of significant rank correlation coefficient is 0.05 at significant level of 0.1%. The correlation coefficients were squared and normalized to calculate contribution factor (k) to output variance from each individual input parameter (MacLeod et al., 2002) r2 k i ¼ Pi 2 ri
ð10Þ
i
As results of sensitivity analysis, this paper reported significant rank correlation coefficients between the input parameters and simulated chemical concentration in compartments of air, root-zone soil, surface water, and sediment. Model sensitivity was evaluated by presenting the sensitivity of total chemical inventories in the Lake Superior basin to all input parameters.
ð6Þ 3. Results and discussion
where r2I i is the variance in the lognormal distribution of the ith input parameter, and r2O;I i is the variance of output contributed by the uncertainty in the ith input parameter. Based on the assumption that uncertainty in all input parameters were independent, total variance in the output (r2O ) can be expressed as sum of variances in the output contributed by uncertainty of individual model input (r2O;I i ). This linear relationship allows to define a contribution factor k (0 < k < 1) to present the contribution of individual input to total variance in output r2O;I i ¼ k i r2O
ð7Þ
As shown in Eqs. (6) and (7), three parameters (input variance, total output variance, and contribution factor) were required in calculating the model sensitivity (S) of a specific model input. Variance of mode input (r2I i ) was defined in Tables 1–3; while total output variance (r2O ) was obtained based on the Monte Carlo simulation. The contribution factors (k’s) were determined from Spearman rank correlation analysis between model inputs and output. By assigning a rank to each observation in the model input and output groups, the Spearman’s method calculates the sums of the squares of the differences in paired ranks !, N X 2 ri ¼ 1 6 dj ð8Þ ðN 3 N Þ
3.1. Simulated concentrations and comparison with field data The simulated concentrations in the bulk compartments during 2000 are presented in Table 4(a) as means and CV Table 4 Annual means of predicted concentrations (g m3) and coefficients of variance (in parentheses) of benzo[a]pyrene and hexachlorobenzene in the Lake Superior basin (a) with variances of all parameters, and with only variations of (b) chemical parameters or (c) environmental parameters Bulk compartment
Benzo[a]pyrene
Hexachlorobenzene
Panel (a) Air Foliage Surface soil Root-zone soil Vadose-zone soil Surface water Sediment
1.23 · 1011 (0.72) 2.64 · 105 (1.14) 1.17 · 104 (1.39) 1.32 · 104 (3.24) 4.01 · 106 (2.79) 4.03 · 109 (1.67) 2.13 · 104 (2.62)
8.51 · 1011 (0.14) 1.84 · 106 (0.44) 2.90 · 105 (1.56) 2.45 · 105 (1.93) 1.32 · 106 (1.94) 7.87 · 109 (0.50) 1.05 · 105 (1.25)
Panel (b) Air Foliage Surface soil Root-zone soil Vadose-zone soil Surface water Sediment
1.21 · 1011 (0.72) 2.52 · 105 (1.07) 1.60 · 104 (1.11) 9.26 · 105 (1.17) 2.21 · 106 (1.28) 4.19 · 109 (1.61) 2.12 · 104 (1.83)
8.32 · 1011 (0.05) 1.74 · 106 (0.29) 3.42 · 105 (0.63) 2.31 · 105 (0.68) 1.06 · 106 (0.70) 7.62 · 109 (0.34) 1.12 · 105 (0.61)
Panel (c) Air Foliage Surface soil Root-zone soil Vadose-zone soil Surface water Sediment
1.37 · 1011 (0.02) 2.17 · 105 (0.25) 1.43 · 104 (0.66) 1.39 · 104 (1.63) 4.69 · 106 (1.10) 4.22 · 109 (0.16) 2.39 · 104 (0.79)
8.65 · 1011 (0.13) 1.67 · 106 (0.31) 2.74 · 105 (1.29) 2.35 · 105 (1.60) 1.40 · 106 (1.18) 8.69 · 109 (0.33) 1.19 · 105 (0.91)
j¼1
where ri is the rank correlation coefficient between the model output and the ith input, d’s are differences between ranks of corresponding model input and output, and N is the sample size in the Monte Carlo simulation. The significance of the rank correlation coefficient was determined by the t-test. For the sample size above about 20, the variable
Y. Luo, X. Yang / Chemosphere 66 (2007) 1396–1407
1403
Table 5 Comparison of model predictions (mean and CV in parentheses) and reported ranges for chemical concentrations Basins
Concentration in air (pg m3) Predicted
Concentration in dissolved water (lg m3)
Concentration in sediment particle (lg kg1)
Reporteda
Predicted
Reportedb
Predicted
Reportedb
Benzo[a]pyrene Superior 12.3 Michigan 25.7 Huron 17.6 Erie 35.4 Ontario 46.6
(0.72) (0.70) (0.75) (0.66) (0.67)
6.0–17.2 8.9–55.1 12.2–28.6 35.5–121.7 5.0–60.0
0.002 0.002 0.003 0.009 0.009
(1.72) (1.55) (1.61) (1.42) (1.56)
0–0.39 NA NA 0–1.1 0–1.48
0.307 1.110 1.970 7.080 6.610
(2.54) (2.22) (2.61) (3.47) (3.44)
NA ND–9.6 NA NA 30–119
Hexachlorobenzene Superior 85.1 Michigan 61.7 Huron 58.3 Erie 64.4 Ontario 12.3
(0.14) (0.13) (0.14) (0.15) (0.14)
62.7–118.2 59.7–128.2 26.8–33.2 62.3–131.7 26.6–40.4
0.008 0.002 0.006 0.004 0.002
(0.50) (0.39) (0.40) (0.35) (0.41)
0.010–0.026 0.007–0.014 0.007–0.033 0.014–0.078 0.045–0.063
0.015 0.011 0.028 0.019 0.010
(1.09) (1.74) (1.46) (1.43) (1.65)
NA NA NA NA 0.04–31.2
NA = not available, ND = non-detect. a Concentrations estimated from measurements in the integrated atmospheric deposition network (IADN, 2006). b Field measurements complied by Zhang et al. (2003) and Booty et al. (2005).
Due to the data limitation, only annual means of chemical concentrations with variations were compared with reported data in this study. The comparison between the model predicted and reported concentrations in the literature during 1991–2000 is shown in Table 5. Predicted concentrations in other compartments or sub-compartments were not tabulated due to a lack of measured data. The predicted and measured concentrations in the air with 5% and 95% percentiles were also plotted in Fig. 2. The model predicted mean air concentration fell within the observed concentration ranges between 5% and 95% percentiles for both chemicals. In this study, the air emission rates were
150
1. B[α]P in the Superior LB 2. B[α]P in the Michigan LB 3. B[ α]P in the Huron LB 4. B[α]P in the Erie LB 5. B[ α]P in the Ontario LB 6. HCB in the Superior LB 7. HCB in the Michigan LB 8. HCB in the Huron LB 9. HCB in the Erie LB 10. HCB in the Ontario LB
120 3
Observed concentration (pg/m )
values. For each chemical compound, the distribution patterns in the compartments were similar over the five basins. About 94% of benzo[a]pyrene was accumulated in the soil compartments. This finding was consistent with the simulation results of CHEMGL and Mackay Level III models (Zhang et al., 2003). With low vapor pressure and high organic carbon partition coefficient, most of benzo[a]pyrene (about 70% in surface water and more than 98% for air, soil, and sediment compartments) was bound to the solid phase. Similarly, more than 90% of hexachlorobenzene was predicted in soil compartments, and about 98% of hexachlorobenzene in soils and sediments were bound to the solid phase. As a volatile chemical, however, hexachlorobenzene in air compartments exists primarily in the vapor phase, and only 3% of hexachlorobenzene was found in aerosol particles. Compared to the air and surface water, porous compartments of soil and sediment showed larger uncertainties in the predicted chemical concentrations. CV values for the predicted concentrations in soil and sediment were 2–8 times lager than those in the air and surface water. As shown in Tables 1–3, high variabilities were assumed for the parameters related to soil and sediment compartments, especially for organic carbon contents and reaction halflives. It was also instructive to compare the uncertainty in the concentration of the chemical tracers in corresponding compartments. When the same set of variabilities in environmental parameters was applied in stochastic simulations, benzo[a]pyrene, compared to hexachlorobenzene, showed higher uncertainty in the predicted concentrations in all compartments except surface soil. In the chemical properties illustrated in Table 3, only Henry’s law constant (H) and organic carbon partition coefficient (KOC) were assumed with higher uncertainties for benzo[a]pyrene than those for hexachlorobenzene. When applied to the input parameter settings in this study, therefore, the multimedia environmental fate model might be highly sensitive to these chemical properties.
90
7 9 6 4
60
10
30
3
5
8
2
1 0 0
30
60
90
120
150
Simulated concentration (pg/m 3 )
Fig. 2. Comparison of simulated and observed air concentration of benzo[a]pyrene and hexachlorobenzene in the Lake Basins (LB) during 1991–2000. The error bars indicate the 5% and 95% percentiles of the data. Observed air concentration were compiled from the dataset of Integrated Atmospheric Deposition Network (IADN, 2006).
1404
Y. Luo, X. Yang / Chemosphere 66 (2007) 1396–1407
quantified and the assumption of well-mixed compartment might be more reliable for air compartment. Therefore, the predicted concentrations in the air were judged to be satisfactory for the present screening purpose. The predicted concentrations in surface water and sediment were lower than observations by a factor of 10–100. This might be caused by assuming zero chemical direct discharge to surface water due to limitation of available data. Furthermore, multimedia environmental fate simulation for basin-average contamination levels in well-mixed compartments may report lower concentrations compared to the mean values of site-specific measurements (Zhang et al., 2003). Even with inaccurate emission data, the predicted concentration also showed consistent distribution pattern as with the reported measurements, especially for benzo[a]pyrene. For example, high concentrations of benzo[a]pyrene in surface water and sediment were predicted in Lake Erie and Lake Ontario. 3.2. Sensitivity and key fate processes The datasets of chemical and environmental parameters were varied independently to gain further insight to separate variability and true uncertainty, as shown in Table 4(b) and (c). The outcome variance resulting from variation of one input group of parameters indicates the influence of uncertainty of this parameters group on the total outcome variance. In Table 4(b), the predicted variations for benzo[a]pyrene resulting from the chemical properties were lager than those for hexachlorobenzene as shown in Table 3. As mentioned previously, the high variations in predicted concentration of benzo[a]pyrene were attributed to the high uncertainties in its chemical properties of Henry’s law constant and organic carbon partition coefficient. Since the same set of landscape parameters were applied in the Monte Carlo simulation (Tables 1 and 2), the variations for the two chemicals resulting from the landscape parameters were generally comparable in corresponding compartments, as shown in Table 4(c). In the comparison of output variations for the two different chemicals, it was noteworthy that magnitudes of the chemical properties had signi-
ficant effects on the nonlinearity of the model simulation. For example, a lager half-life in the air would result in a lager variation in predicted chemical concentrations in air given datasets for other inputs were invariant over simulations. For benzo[a]pyrene, the output variances resulting from the chemical parameters was much greater than those from the environmental parameters for all compartments except for root-zone soil. For hexachlorobenzene, in contrast, the environmental parameters contributed more to the output variances compared to the chemical properties. For both chemicals, higher variations of predicted concentrations in the soil layers were observed relative to other compartments. This indicated that if soil was the most important exposure pathway, the values of landscape parameters were the most influential. According to the significance of rank correlation defined in Eq. (9), all rank correlation coefficients larger than 0.05 were significant at 0.1% level in this study. To show the parameters which were most strongly correlated with model outputs, an arbitrary value of 0.3 was chosen as critical value of the coefficients. Table 6 shows the rank correlation coefficients of primary input parameters (with absolute value of the coefficient larger than 0.30) on the predicted concentrations in the Lake Superior basin. To simplify the stochastic simulation considerably, all parameters were assumed to be independent or uncorrelated. In reality, some inputs were correlated, especially for the degradation rate constants and hydrological flow rates. This assumption and simplification were consistent with techniques previously applied for contaminant fate models (Kuhne et al., 1997; Liu et al., 1999). Introducing covariance would not affect the results of sensitivity analysis if some parameters, under independent assumption, were identified as important contributor to output variance. The predicted concentrations in the air and surface water were influenced primarily by the chemical properties of corresponding half-lives (HLa and HLw), indicating that the chemical degradation was the most important fate process in these compartments. On the other hand, chemical distributions in porous compartments of root-zone soil
Table 6 The input parameters with Spearman rank correlation coefficients (in parentheses) with absolute values greater than 0.3 for predicted concentrations in air, root zone, surface water, and sediment in the Lake Superior basin Chemical
Primary influential parameters in bulk compartment Air
Root-zone soil
Surface water
Sediment
Benzo[a]pyrene
HLa (1.00)
fOCsp (0.65) HLa (0.39) KOC (0.32)
HLa (0.89) HLw (0.42)
HLa (0.74) HLw (0.35) KOC (0.35) fOCdp (0.32)
Hexachlorobenzene
LAI (0.56) vw (0.50) Dair (0.49)
fOCsp (0.83) KOC (0.31)
HLw (0.64) vw (0.44) VP (0.38)
fOCdp (0.61) KOC (0.44) fdw (0.40)
Note: HLa = half live in air, HLw = half live in surface water, KOC = organic carbon partition coefficient, fOCsp = organic carbon content in root-zone soil solid, and fOCdp = organic carbon content in sediment particle, Dair = diffusivity in air, vw = wind speed, LAI = leaf area index, VP = vapor pressure, and fdw = pore water content in sediment.
Y. Luo, X. Yang / Chemosphere 66 (2007) 1396–1407
S= 1.
0
1.2
Standard deviation of output
1
0.8 HLa 0.6
Koc
focsp 0 .4 S=
H transpiration
0.4
fgw
0.2
0 0
0.2
0.4 0.6 0.8 Standard deviation of inputs
1
1.2
S= 1. 0
1.2
1 Standard deviation of output
and sediment were primarily affected by organic carbon partition coefficient and organic carbon contents in the particles (fOCsp and fOCdp). This finding was consistent with the results in the separation of variability and true uncertainty in Table 4(b) and (c). For the plant foliage (not shown in Table 4), the outcome variances were contributed mainly by the same parameters which were strongly associated with predicted concentrations in either air or rootzone soil, because these compartments interact with each another closely. Similarly, the predicted concentrations in sediment shared the same input parameters which have significant effects on the variance of concentrations in surface water. Without direct source input, water-sediment diffusion and settling from water column were the controlling transport processes of these chemicals in sediment compartment. Rank correlation results supported the identification of the key fate and transport processes for chemicals. Being released to air and removed rapidly by degradation, for example, benzo[a]pyrene concentrations in all compartments were sensitive to its half-life in air (HLa). For hexachlorobenzene, the predicted air concentration was sensitive to wind speed (vw) and leaf area index (LAI) as shown in Table 6. This indicated that the mass loss via advective air flow and mass exchange with plant foliage were the dominant transport mechanisms in the air. It was noteworthy that the parameter of wind speed reflected air advection and chemical inputs from upwind regions as well. With explicit spatial simulation of chemical transport over basins, the parameter uncertainties in one compartment were allowed to be propagated to other compartments within the basin, and to neighbor basins as well following the air and water connectivites. In soil compartments, the predicted concentrations of benzo[a]pyrene and hexachlorobenzene were sensitive to organic carbon content in soil solids and organic carbon partition coefficient (fOCsp and KOC). Diffusive fluxes between soil compartments were the predominant processes over advections. Further calculations showed that the flux ratio of advection to diffusion was about 1:9 and 1:17, for benzo[a]pyrene and hexachlorobenzene, respectively. A graphic analysis of contribution to uncertainty in predicted total mass in the Lake Superior basin is shown in Fig. 3. The output variances contributed by an individual parameter (r2O;I i ) were calculated in Eq. (7). For a prescribed set of input parameters, the sensitivities were defined properties of the model. If it was desired to reduce the output standard deviation, therefore, the standard deviations of input should be decreased by traveling down lines of constant S values. To reduce output uncertainty more efficiently, the input parameters with high sensitivity and large assumed CV should be targeted. The input parameters or intermediate variables were labeled if their sensitivity values were larger than 0.4. High sensitivities were also observed for some input parameters which had small standard deviations (
