Dynamic modeling of chemical fate and transport in multimedia environments at watershed scale—I: Theoretical considerations and model implementation

ARTICLE IN PRESS

Journal of Environmental Management 83 (2007) 44–55 www.elsevier.com/locate/jenvman

Dynamic modeling of chemical fate and transport in multimedia environments at watershed scale—I: Theoretical considerations and model implementation Yuzhou Luoa, Qiong Gaob, Xiusheng Yanga, a

Department of Natural Resources Management and Engineering, University of Connecticut, Storrs, CT 06269-4087, USA b Institute of Resources Science, Beijing Normal University, Beijing 100875, China Received 7 April 2005; received in revised form 19 January 2006; accepted 26 January 2006 Available online 11 May 2006

Abstract A geo-referenced environmental fate model was developed for analyzing unsteady-state dispersion and distribution of chemicals in multimedia environmental systems. 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 the compartments were completely mixed and chemical equilibrium was established instantaneously between the sub-compartments 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. A 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. Temporal variations of the environmental properties and source emissions were also considered in the parameter estimations. 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. r 2006 Elsevier Ltd. All rights reserved. Keywords: Chemical dynamics; Fugacity; GIS; Multimedia environment; Watershed scale

1. Introduction Effective environmental management requires an understanding of and capability to quantitatively analyze the transport mechanism and health risks associated with chemical contaminations. The level of contamination resulting from accidental or continuous chemical releases is usually assessed with the use of environmental transport models. These models are traditionally based on a single environmental pathway, such as the transport of a contaminant through water (Coulibaly et al., 1998). However, chemicals released to the environment are cycled through various environmental media as the consequence of complex physical, chemical, and biological processes. Corresponding author. Tel.: +860 486 0135; fax: +860 486 5408.

E-mail address: [email protected] (X. Yang). 0301-4797/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.jenvman.2006.01.017

Many real-world environmental problems involve processes that occur both within and between environmental media (e.g., air, soil, surface water, and biota). The traditional approach in simulating transport processes across environmental media is to couple existing singlepathway models (Charnock et al., 1996). According to Brandmeyer and Karimi (2000), there are five types of methods in model coupling, i.e., one-way data transfer, loose coupling, shared coupling, joined coupling, and tool coupling. Typically, coupled environmental models are built upon or expanded from single-medium simulations. Outputs from one medium are used as inputs for another adjacent medium, often termed as ‘‘boundary conditions’’ or ‘‘driving functions’’. Inter-media mass transport fluxes are thus simplified as fixed values or prescribed time-series determined by processes within a single medium (Cohen and Cooter, 2002a). Therefore, these coupled single-media

ARTICLE IN PRESS Y. Luo et al. / Journal of Environmental Management 83 (2007) 44–55

models are not capable of analyzing the interactions and mechanisms of chemical transport in real multimedia environments. Multimedia modeling has been introduced as a new tool to track the distribution of a contaminant in all environmental components connected together, by incorporating not only the single-medium fate and transport but the inter-media transfers as well. The most frequently used multimedia models are fugacity-based mass balance models (or Mackay-type models) that focus on the partitioning of chemicals among environmental media. The fugacity concept was firstly used in the QWASI model (Mackay and Yeun, 1983) to assess the distribution of contaminants in the Great Lakes, and has been widely adopted in environmental modeling over the last decade (McKone, 1993; Cowan et al., 1995; Devillers et al., 1995; Diamond et al., 2001; ChemCAN, 2003; Zhang et al., 2003). Most of these applications were designed to describe the steadystate response to chemical inputs (so-called level III models), and treat the entire simulation domain as a unit environment without spatial variations. There is a small, but growing, number of multimedia environmental fate models that consider more accurate descriptions of the environment and analyze chemical dynamics (MacLeod et al., 2001; Woodfine et al., 2001; Sweetman et al., 2002, 2004; Coulibaly et al., 2004a,b). However, these models do not support the utilization of time-dependent environmental parameters such as temperature, wind, precipitation, and water flows. Instead, long-term averages of environmental parameters are used in these models to simulate chemical dynamics. The lack of spatial and temporal variations in the environmental properties thus limits their ability to quantify chemical distribution and persistence in the multimedia environment (Pennington, 2001; Luo and Yang, 2004). This study aimed to develop a multimedia environmental fate model for chemical dynamics with spatial resolution at watershed scale. In this first part of a two-paper series, theoretical considerations and model implementations of a general multimedia modeling framework are presented. This model relied on fugacity to describe the mass potential of chemical species in an Eulerian system, in which the region of interest was divided into a number of connected boxes, or compartments, representing the connected environmental media. The chemical transport and fate were simulated by a set of fugacity-based transport equations to mimic inter-media chemical fluxes and dynamic mass balance in various environmental media, including air, plant, soil, water, and sediment. The model was formulated and parameterized with the following requirements: (1) inclusion of all major mass transport and transformation processes at watershed scale, (2) integration of all environmental media to reflect intermedia interactions, and (3) estimations of model parameters with spatial and temporal variations. This approach enabled holistic integration of modeling capabilities that have not been fully achieved in existing models, in terms of

45

inter-media mass transport simulation, spatial resolution at watershed scale, and time-dependent estimation of chemical distribution. 2. Mathematic formulations 2.1. Domain and compartments The simulation domain of this study was an entire river basin. It included basic elements of atmosphere, terrestrial and aquatic biota, unsaturated soil, surface water, and sediment, with the top of the troposphere as the upper boundary, and the bottom of the vadose-zone soil as the lower boundary. The simulation domain was horizontally segmented into n small regions depending on the size of the domain. The segmentation of regions followed the delineation of watersheds to minimize technical complexity in handling water flows between adjacent watersheds. The threshold areas in the watershed delineation were in the range of 102–104 km2. This spatial resolution was used to ascertain that a chemical was likely to spend enough time in each watershed to allow reactions and inter-media transport to occur. The environment in each watershed was divided into a number of boxes or compartments linked by a variety of inter-media transport processes. Seven major compartments, including the atmosphere, plant canopy, surface soil, root-zone soil, vadose-zone soil, surface water, and sediment, were considered in each watershed (Fig. 1). The total number of compartments was m ¼ 7n where n is the number of watersheds delineated in the simulation domain. These environmental compartments were considered well mixed and homogeneous, in terms of environmental properties and chemical concentration. Each compartment included different sub-compartments characterized by their physical properties. For example, the surface water compartment has three sub-compartments of pure water, suspended particles, and aquatic biota. The air compartment horizontally covered the area of the corresponding watershed and vertically extended from the ground surface to the top of the troposphere. This compartment consisted of pure air and suspended aerosol particles. The plant canopy compartment was separately set over the soil and comprised forest, cropland and pasture. Characterized by the coverage and the mass of plant canopy, this compartment included only the aboveground portion of plants. The portion of plants below the ground was considered in the root-zone soil compartment. The archetypal structure of soil layers was aggregated into three well-mixed soil compartments. In the absence of tilling, particles deposited from the atmosphere were accumulated in and resuspended from a thin surface soil layer (0.1–1 cm). The root-zone soil was below the surface and encompassed the region capturing the plant-rooting zone and the maximum diffusion depth. The vadose-zone soil compartment was defined from the bottom of the rootzone soil to the top of the groundwater table. All soil compartments consisted of air, water, and particles in soil.

ARTICLE IN PRESS Y. Luo et al. / Journal of Environmental Management 83 (2007) 44–55

46

Watershed #1

Watershed #2

Air

Air

Surface Soil Root-zone Soil

Canopy Root

Vadose-zone Soil

Water

Water

Groundwater Sediment Advection and physical interfacial transport, e.g. transport via inflow/outflow deposition runoff uptake... Diffusive transport processes Fig. 1. Environmental description with primary inter-media mass transfer processes (details in Table 1).

The root-zone soil contained plant roots as well. Transport parameters in soil were derived by matching compartment inventories to those obtained from analytical solution (McKone, 1993, 1996; Bennett et al., 1998). The surface water compartment was the surface water bodies in a watershed. This compartment comprised pure water, suspended solids, and aquatic biota. The sediment compartment, with an area extent equal to that of the surface water compartment, consisted of pure water and particles in the top active layer of sediment (2–5 cm) where active contaminant exchange occurred with the overlying water column. The current version of the model did not consider the upper layer atmosphere (the stratosphere) and coastal aquatic environment. The chemical flux across the troposphere–stratosphere was neglected by assuming that the bidirectional fluxes are balanced by each other. The migration and dilution of chemicals in groundwater was not explicitly simulated. Instead, the contaminant leaching from the vadose-zone soil was considered as an input to the groundwater. Due to the complexity in describing the structure of fractured permeable media, we neglected the bulk transport from air to soil resulting from barometric fluctuations (Nilson et al., 1991; William et al., 1997). Snow melting and snow scavenging were considered in the model. The influence of snow pack on the air–ground interactions, however, was neglected in the current phase of the study.

2.2. Overview of mathematical formulations Under the assumption of homogenous composition of the environmental compartments, the following transport and transformation mechanisms were considered in

developing mass balance equations for multimedia chemical dynamics. (1) mass exchange between species, i.e., chemical degradation or mass gain from parent compounds, (2) mass exchange between compartments within a watershed, i.e., inter-media mass transport, (3) mass exchange between watersheds resulting from the advective flows of air and water, (4) mass exchange between the simulation domain and the external environment, including distant chemical inputs and outputs by advective flows of air and water, sediment burial, and chemical loss with groundwater recharge, and (5) source emission of chemicals. By considering all these processes, a general unsteadystate mass balance can be described with the following differential equation: m X dN i ¼ Si þ ðQji  Qij Þ  ðQix þ QRi Þ, dt j¼1

i ¼ 1; 2; . . . ; m,

ð1Þ

where Ni (mol) is the chemical inventory in the compartment i at time t, m is the total number of compartments defined in the simulation domain, Si (mol s1) is the total chemical source in this compartment, Qij and Qji (mol s1) are the rates of unidirectional chemical flux from compartment i to j, and vice versa, and QRi (mol s1) is the chemical degradation rate in i. For a compartment located on the boundary of the simulation domain, Qix (mol s1) is the chemical loss rate from i to a hypothetical receptor compartment, x, in the external environment outside of the simulation domain. All these transport and

ARTICLE IN PRESS Y. Luo et al. / Journal of Environmental Management 83 (2007) 44–55

transformation processes were represented mathematically as first-order equations based on the fugacity concept discussed in the following section. 2.3. Fugacity concept

Table 1 Inter-media mass transport processes recognized and formulated in this study Interface

Transport processes

Air–canopy Air–surface soil

Diffusion Dry deposition Wet deposition by rain and snow Wind resuspension

Air–water

Dry depositiona Dry deposition Wet deposition by rain and snow

Surface soil–vadose zone

Diffusion Infiltration

Surface soil–surface water

Overland flow Soil erosion

Root zone–vegetation

Plant uptake Phloem flow

Root zone–vadose zone

Infiltration Flow form vadose zone to root zone

Vadose zone–groundwater

Recharge

Vadose zone–surface water

Interflow

Surface water–bottom sediment

Diffusion Sedimentation Resuspension

Mathematically, the fugacity is defined as (Mackay et al., 1983) f ik ¼

C ik , Zik

(2)

where fik (Pa) is the fugacity of a chemical in a subcompartment k of compartment i, Cik (mol m3) and Zik (mol Pa1 m3) are the corresponding concentration and fugacity capacity, respectively. Formulations for fugacity capacities in various sub-compartments can be found in Coulibaly (2000) and Mackay (2001). In each compartment, a local equilibrium was assumed to exist between all the sub-compartments. Therefore, the fugacity of a chemical in a homogenous compartment can be calculated as (Hemond and Fechner-Levy, 2000) f i ¼ f i1 ¼ f i2 ¼    ¼ f ik P ðC ik V ik Þ Ni k ¼ fi ¼ P , ðZ ik V ik Þ Z i V i

ð3Þ

k

47

a

where fi is the fugacity of chemical in compartment i, N (mol) is the total inventory of the chemical in i, Vk (m3) is the volume of the sub-compartment k, Vi (m3) is the total volume of compartment i, and Zi (mol Pa1 m3) is the bulk fugacity capacity in compartment i. 2.4. Inter-media mass transport fluxes Transport equations for various media were coupled through the formulation of inter-media transport processes, including advective, diffusive, and physical interfacial processes (Table 1). Based on the fugacity concept, the overall inter-media flux from compartment i to j was quantified by the D values that were described in detail by Mackay (2001) Qij ¼ Dij f i ¼ ðDDij þ DAij þ DPij Þ f i , 1 1

(4)

where Dij (mol Pa s ) is the overall transport coefficient from compartment i to j, and DDij, DAij, and DPij (mol Pa1 s1) are the transport coefficients for advective, diffusive, and physical interfacial transport processes from i to j, respectively. In addition to the D values defined by Mackay (2001), this model also included transport mechanisms in plant canopy and soil layers. It was noteworthy that Qij presented only the unidirectional flux rate from i to j, while the net mass exchange by the inter-media diffusion was given by the algebraic sum of Qij and Qji. The D value of the inter-media diffusion was formulated based on the two-film theory (Whitman, 1923; Lewis and

Dry deposition was defined for aerosol particles only. The dry deposition of gases was formulated as a diffusion process.

Whitman, 1924):   diji dijj 1 þ , DDij ¼ Aij Z i Dti Zj Dtj

(5)

where Aij (m2) is the interface area between compartment i and j, d’s (m) are the boundary layer depths of the two compartments at the interface, and Dt’s (m2 s1) are the bulk diffusivities of the chemical in the respective compartments. The boundary layer depth of each compartment in the diffusive transport was estimated using empirical or semi-theoretical equations from the literature (Jury, 1983; McKone, 1993; Thibodeaux, 1995; Bennett et al., 1998; Meyers et al., 1998; Mackay, 2001; Cohen and Cooter, 2002a, b). Based on the volume fractions of air and water in a compartment, the bulk diffusivity can be calculated from the diffusivities of the chemical in pure air and water (Coulibaly, 2000). The inter-media advection was driven by hydrologic flows and plant–soil interactions (plant uptake and phloem flow) as shown in Table 1. The D value of the inter-media advection was calculated as DAij ¼ Aij Z c uij ,

(6) 1

3

where Zc (mol Pa m ) is the fugacity capacity of the chemical in the carrying media (water or phloem fluid), and uij (m s1) is the flow velocity. Eq. (6) was also valid in simulating the advective chemical fluxes between two connected watersheds by air/water flows. For calculating

ARTICLE IN PRESS Y. Luo et al. / Journal of Environmental Management 83 (2007) 44–55

48

the advective chemical fluxes, Zc denotes the fugacity capacity of the chemical in air or water, and uij is the advective flow velocity of air or water across the watershed boundaries. The physical interfacial processes involve chemical transport by particles in air or in water, and cannot be categorized as either advection or diffusion. In this study, these processes included chemical transports by atmospheric dry deposition of aerosol particles, wind resuspension of particles from surface soil and plant foliage, soil erosion, sedimentation and resuspension of solids in water bodies (Cohen et al., 1990; Cohen and Cooter, 2002a). The general formulation of the D values for these transport processes can be expressed as DPij ¼ Aij ðgip Z ip Þ u0ij ,

(7)

where the index i represents the originating compartment from which chemicals are removed, including air (for atmospheric deposition), surface soil (for wind resuspension and soil erosion), or sediment (for sediment resuspension and burial), gip (dimensionless) is the volume fraction of particles in compartment i, Zip (mol Pa1 m3) is the fugacity capacity of the chemical in particles in compartment i, and u0 ij (m s1) is the velocity of the corresponding interfacial transport. These velocities were estimated based upon empirical and semi-theoretical methods from the literature (Ackers and White, 1975; Cowherd et al., 1985; Hicks et al., 1987; Meyers and Baldocchi, 1988; Meyers et al., 1998). For computational convenience, the mass balance was formulated in term of chemical inventory by introducing an overall transfer rate constant. Eq. (4) was rearranged as Qij ¼ M ij N i ; with N i ¼ Z i V i f i ; and M ij ¼ Dij =ðZi V i Þ,

ð8Þ

1

where Mij (s ) is the overall transfer rate constant of the inter-media chemical flows from compartment i to j. Similarly, Mji (s1) was defined as the overall transfer rate constant of Qji. The transfer rate constant for an individual transport process was given by M Dij ¼ DDij =ðZ i V i Þ; or M Pij ¼ DPij =ðZ i V i Þ,

M Aij ¼ DAij =ðZ i V i Þ, ð9Þ

where MDij, MAij, and MPij are the transfer rate constants for diffusive, advective, and physical interfacial processes, respectively. The transfer rate constant indicated the fraction of the chemicals in the originating compartment that was removed by the corresponding transport process per unit time. For example, if a transfer rate constant of a dry particulate deposition was found to have the value of 0.1 s1, then 10% of the chemical in aerosol particles would be deposited per second. The transfer rate constants of inter-media transport processes were summarized in Table 2. The reciprocal of the transfer rate constant (M1) was the characteristic time of the corresponding transport process, i.e., the time required for the transport process to reduce the fugacity to e1 of the original fugacity of a chemical in a compartment.

Chemical losses to the external environment outside of the simulation domain were formulated in a similar way. These processes included chemical burial in sediment, and advective/diffusive transport in air or water compartments located on the boundaries of the simulation domain. The overall rate constant (Mix) of the chemical loss from compartment i to an external sink was defined as M ix ¼

DDix þ DAix þ DPix , Zi V i

(10)

where the subscript x denotes a hypothetical receptor compartment in the external environment connected to i, and DDix, DAix and DPix are the corresponding D values for diffusion, advection, and sediment burial, respectively. 2.5. Other sources and sinks The rates of transformation or degradation of the chemical were estimated via first-order kinetics, which implied that the substrate concentration was the primary factor in affecting the decomposition rate, and the microbial biomass was always present in high enough concentrations not to be rate limiting. As for reactions that do not follow first order kinetics, a pseudo first order reaction was assumed. The degradation flux of a chemical in specific compartment i, QRi (mol s1), was calculated by QRi ¼ M Ri Zi V i f i ,

(11)

1

where MRi (s ) is the degradation rate of a chemical in compartment i. If transformations between inter-converting species were applicable for the modeled species, the mass gain (QTi, mol s1) from a parent compound was QTi ¼ M 0Ri Z0i V i f 0i , f 0i

(12)

0

where, and Z i are the fugacity and fugacity capacity of the parent species in compartment i, respectively, and M0 Ri (s1) is the degradation rate of the parent species. For watersheds on the boundary of the simulation domain, chemical inputs by advective flows of air or water from external environment were taken as distant sources (QIi, mol s1). In a compartment, a total source term (Si, mol s1) was defined as the sum of the chemical inputs from transformation (QTi), distant transport (QIi), and emission source (QSi), Si ¼ QTi þ QSi þ QIi .

(13)

2.6. Mass balance equations and numerical method Substituting Q’s in terms of N and M, the mass balance equation in (1) was rearranged as ! 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 , j¼1

ð14Þ

ARTICLE IN PRESS Y. Luo et al. / Journal of Environmental Management 83 (2007) 44–55

49

Table 2 Overall transfer rate constants for inter-media transport processes Inter-media transport Air to canopy

Air to surface soil

Air to water

Canopy to air

Overall transport rate constant M ap ¼

LAI DDap þ Ap I w urain ðZwater þ gap Zap QÞ þ Ap I d gap udepo Zap Za V a

M ag ¼

DDag þ Ag urain ðZwater þ gap Z ap QÞ þ Ag gap udepo Zap Za V a

M aw ¼

DDaw þ Aw urain ðZwater þ gap Zap QÞ þ Aw gap udepo Zap Za V a

M pa ¼

LAI DDap þ Ap gap uwindres Zap Zp V p

Canopy to root zone

M ps ¼ Ap uphlm Zphlm =ðZp V p Þ

Canopy to surface soil

M pg ¼ Ap ulitter Zp =ðZp V p Þ

Surface soil to air

Surface soil to root zone

Surface soil to water

M ga ¼

DDag þ Ag ggp uwindres Zgp Zg V g

M gs ¼

DDgs þ Ag uinfil Zwater Zg V g

M gw ¼

Ag ðurunoff Zwater þ ggp uerosion Zgp Þ Zg V g

Root zone to canopy

M sp ¼ Ag uuptake Zwater =ðZs V s Þ

Root zone to surface water

M sg ¼ DDsg =ðZs V s Þ

Root zone to vadose zone

M sv ¼

Vadose zone to root zone

M vs ¼ DDsv =ðZv V v Þ

Vadose zone to water

M vw ¼ Ag uinter Zwater =ðZv V v Þ

Vadose zone to groundwater

M vq ¼ Ag urech Zwater =ðZv V v Þ

Water to air

M wd ¼ DDwd =ðZw V w Þ

Water to sediment

Sediment to water

Burial to deep sediment

DDsv þ Ag uperc Zwater Zs V s

M wd ¼

DDwd þ Aw gwp usedmt Zwp Zw V w

M dw ¼

DDwd þ Aw gdp uresus Z dp Zd V d

M dx ¼ Aw gdp uburial Zdp =ðZd V d Þ

Where LAI (dimensionless), leaf area index; Iw and Id (dimensionless), canopy interception fractions of wet and dry deposition; Zwater (mol Pa1 m3), fugacity capacity of the chemical in pure water; Zap, Zgp, Zwp, and Zdp (mol Pa1 m3), fugacity capacity of the chemical in aerosol particles, in surface soil solids, in suspended particles of water, and in sediment solids, respectively; gap, ggp, gwp, and gdp (dimensionless), volume fraction of particles/soilds in air, in surface soil, in water, and in sediment, respectively; Q (dimensionless), particle scavenging ratio for rain or snow; urain (m s1), precipitation rate; udepo (m s1), dry deposition velocity of aerosol particles; uwindres (m s1), wind resuspension velocity of particles from canopy and surface soil; uphlm (m s1), velocity of fluid that moves from canopy tissues down into the roots through the phloem tubes; ulitter (m s1), velocity of leaf litterfall; uinfil (m s1), infiltration rate; urunoff (m s1), surface runoff rate; uuptake (m s1), velocity of water that moves from soil into the roots and up through the plant as a result of transpiration; uperc (m s1), percolation rate, from root zone to vadose zone; uinter (m s1), lateral inter-flow rate, from vadose zone to water; urech (m s1), groundwater recharge rate; usedmt (m s1), velocity of particle sedimentation; uresus (m s1), velocity of resuspension of sediment solids; ubural (m s1), velocity of sediment burial. Indices a, p, g, s, v, w, and d denote the compartment of air, plant canopy, ground surface soil, root-zone soil, vadose-zone soil, surface water, and sediment, respectively, M’s (s1) and D’s (mol Pa1 s1) are the transfer rate constants and Mackay D values of transport processes, respectively, A’s (m2) are the horizontal projective areas of compartments, V’s (m3) are the compartment volumes, and Z’s are the fugacity capacities of chemical in compartments or sub-compartments. The areas of surface soil, root-zone soil, and vadose-zone soil were assumed to be the same and denoted as Ap. Similarly, the areas of surface water and sediment were assumed to taken the same value of Aw.

ARTICLE IN PRESS Y. Luo et al. / Journal of Environmental Management 83 (2007) 44–55

50

where MOi (s1) is the total loss rate of the chemical in compartment i. For the whole simulation domain, Eq. (14) was written for m compartments to solve for the timedependent chemical inventories. Each of these equations was re-written in implicit forms with the finite difference method and solved for N’s by ð1=Dt þ M Oi ÞN i jtþ1 

m X

  ðM ji N j Þtþ1 ¼ ðS i þ N i =DtÞt ,

j¼1

i ¼ 1; 2; . . . ; m.

ð15Þ

The model was implemented with both capabilities for providing steady-state and time-dependent solutions. At steady state, sources were in balance with sinks and there was no chemical accumulation in each compartment. Consequently, the set of ordinary differential equations were transformed to a set of linear equations expressed in matrix form. The numerical solution was implemented with MATLAB codes in IBM PC platform. 3. Model implementation 3.1. Model structure and input requirements In general, the model consisted of three functional layers: spatial characterization, parameter estimation, and multimedia fate simulation. Spatial analysis was performed by geographic information system (GIS)-powered functions, which were utilized in the landscape characterization to assign environmental properties and hydrologic simulation results to each watershed. The function layer of parameter estimation was designed to quantify the spatial

Henry’s law constant Temperature, Vapor pressure, Densities, Partition coefficients, OC contents

and temporal variations of the environmental properties. In the multimedia fate and transport simulation, the D values of the inter-media transport processes were estimated based on chemical and environmental properties (Fig. 2). The resultant D values were converted to transfer rate constants by considering interface area and compartment dimensions (Table 2). Then the transfer rates were combined with inventories estimated in the prior time step to compute inter-media transport fluxes and loss rates, and used to solve for the inventories in the current time step from the dynamic mass balance equation (Eq. (15)). The model was run at daily time increments with data inputs for all watersheds at each simulation step. The input information required for the model simulation comprised environmental properties, chemical properties, chemical releases, and background concentrations. Environmental properties included landscape parameters and meteorological data. The simulation domain and its watersheds were characterized by the landscape characteristics, such as connection relations, compartment dimensions, volume fractions of sub-compartments, and organic carbon contents in particles (Table 3). Meteorological data included temperature, wind speed, wind direction, precipitation, relative humidity, and solar radiation. These meteorological data were used to drive a hydrological module to simulate the fluxes of evapotranspiration, infiltration, percolation, surface runoff, interflow, groundwater recharge, and stream flow. For a chemical species, the required properties included the molecular weight, vapor pressure, octanol–water partition coefficient, sorption coefficient, diffusivities in air and water, and Henry’s law constant. The degradation rate constant or half-life of the

Sub-compartment fugacity capacity Volumetric fractions of air, water, particle Bulk fugacity capacity

Diffusivity in air, Diffusivity in water Bulk diffusivity

Reaction half-lives

Diffusive D values

Air/water flows

Advective D values Compartment dimensions

Transfer rates

Emission rates, Initial conditions Chemical properties, Source emissions, Background concentrations

Mass balance equation

Air/water connectivities

Model simulation

Landscape characteristics, Meteorological and Hydrological data

Fig. 2. Schematic of multimedia transport simulation with input requirements.

ARTICLE IN PRESS Y. Luo et al. / Journal of Environmental Management 83 (2007) 44–55 Table 3 Major environmental characteristics used by the multimedia model Environmental compartment

Dimensional and physical parameters

All compartments

Surface area Depth or thickness Volume fractions of sub-compartments Densities of subcompartments Organic carbon content

Air

Particles size distribution Atmospheric scavenging ratio

Interfacial zones

Air–soil Air–surface water Air–canopy

Canopy

Canopy type Forest inventory Leaf area index

Canopy–air Canopy–surface soil Canopy–root-zone soil

Soils

Land slope, land use Porosity, permeability Initial water content

Soil–air Soil–surface water Soil–canopy

Surface water and sediment

Suspended solid size and sphericity Bio-accumulation factor

Surface water–air

51

With watershed segmentation, soil and plant canopy compartments were considered to be isolated from the neighborhood, i.e., there was no direct chemical transport between the soil or canopy compartments in any two adjacent watersheds. The connectivity structures in water, soil, and canopy compartments were considered to be invariant with time. However, the connectivity of air compartments is time dependent, and changes with wind direction. The current segmentation of air compartments based on watershed delineation (Fig. 1) was not suitable for horizontal transport of chemicals in the atmosphere. Based on wind direction, advective air outflows from one watershed may affect more than one downwind watershed. Therefore, the atmosphere was re-segmented as interconnected grid cells (Fig. 3a). The grid size can be set to user-defined values, and different cell size can be used for a nesting area. GIS

Surface water–sediment

chemical species in a compartment was an overall parameter used for a variety of processes, such as photolysis, hydrolysis, and microbial breakdown. Source emission was considered to be uniformly distributed in a compartment. The implementation of the model also required the environmental chemical concentrations adjacent to the simulation domain, for the evaluation of the distant chemical inputs by advective flows of air or water entering the domain (Cohen and Cooter, 2002b).

(a)

Sub-watershed

Air grid cell

3.2. Landscape description and compartment connectivities The connectivities of environmental compartments were defined in the model simulation design. Watershed delineation and water compartment connectivity was based on the surface hydrologic analysis developed by Jenson and Domingue (1988). Using the National Hydrography Dataset (NHD) as a reference, watersheds were delineated from digital elevation data (DEM) in the simulation domain. The connection relations from upper to lower river segments were given by calculating the flow accumulation in the watersheds. Lakes were incorporated into the river-networking system as a special river segment. Digital maps of ground and soil characteristics were overlapped in this function to aggregate hydrologic parameters on a watershed basis. The resultant topological relations and derived environmental properties enabled the evaluations of spatial and temporal variations in the simulations of hydrology and chemical transport.

A D

B C Air Column

(b)

Interactive surface of soil, canopy, or water

Fig. 3. Schematics of (a) air connectivity design, and (b) inter-media transport calculation based on projective areas between air grid cell and compartments at ground.

ARTICLE IN PRESS 52

Y. Luo et al. / Journal of Environmental Management 83 (2007) 44–55

functions were used to generate the grid system, and to calculate the projective area between the grid cells and the underlying compartments of soil, canopy, and surface water. Calculations of the horizontal transport and air–ground interactions were conducted using the following procedures. (1) input meteorological data were interpolated and assigned to each grid cell by Ordinary Kriging (Jarvis and Stuart, 2001; Bai and Feng, 2003); (2) daily averages of the meteorological data were used to estimate the advective and diffusive transport parameters using the effective velocity strategy (Strand and Hov, 1993). The transport parameters then were used to calculate the advective and diffusive transport between adjacent air grid cells; (3) for each of the patches (e.g., A, B, C, and D in Fig. 3b) intercepted by air grid cell and ground compartments, inter-media transport fluxes were calculated by the Eqs. (5)–(7); (4) the fluxes from (2) and (3) were incorporated with chemical transformation and emission into the mass balance equation, Eq. (14), to calculate the chemical inventories in each air grid cell; (5) chemical inventories in air and inter-media fluxes between air and ground were reported on a watershed basis. These values are calculated based on the projective area of air grid cells on the ground compartments in a watershed (Fig. 3b).

3.3. Model parameter estimation at watershed scale At every time step, fate and transport of the chemical species were simulated with representative environmental parameters for the compartments defined in the simulation domain. A number of spatial functions were utilized based on GIS technology for parameter estimations in the multimedia environmental system. These functions were developed in ArcGIS 9.0 with its extensions of Spatial Analyst, Network Analyst, and Geostatistical Analyst (ESRI, 2001). For parameters available from GIS digital maps, areaweighted or length-weighted values were calculated in each watershed. These parameters included compartment areas, land use, river hydraulic properties, leaf area index, and soil properties. Aerosol loads measured at monitoring sites were interpolated into the air grid cells defined for the air connectivity. Chemical emission data were mainly provided in the form of ‘‘facility street address’’. A geo-coding technique was used to locate the emission site and sum the emission rates into each watershed. Similarly, particle load and organic carbon content in the surface water from USGS gages were rearranged for each watershed. Arithmetic averages of these measurements within a watershed were used as representative values of the corresponding parameters.

The Hydrological Simulation Program—Fortran (HSPF) version 12 was adopted to simulate the hydrologic processes at watershed scale (Bicknell et al., 1997; Duda et al., 2001). HSPF was parameterized with the environmental properties resulting from the watershed delineation function, and calibrated in terms of river flows by comparing the simulation results with observed streamflow rates from the USGS NWIS sites (Bicknell et al., 1997). The major hydrological outputs for each watershed included actual evapotranspiration, stream flow, surface runoff, infiltration, interflow, percolation, and groundwater recharge. Flow velocities were determined from these hydrological flows and used in estimating the advective D values and transfer rate constants (Table 2). Plant transpiration (uuptake) and phloem flow (uphlm) were estimated from the calculated evapotranspiration rates (McKone, 1993). If snow was significant in a watershed, the SNOW module in HSPF was activated. The aggregated environmental properties from the geographic database or monitoring measurements were used to calculate the D values and other transport parameters for model simulation in each watershed (Fig. 2). The implied assumption was that the spatial heterogeneity would be propagated linearly in the simulation. This assumption was valid for most of the equations in this study, including the D value formulations for advective and physical interfacial processes, transport and transformation fluxes, and dynamic mass balances. For the non-linear formulation of the D values for diffusive transport in Eq. (5), however, this simple lumping of environmental properties could introduce a bias. In this study, this error was evaluated based on the specifics of the non-linear equations and the heterogeneity of the environmental properties (King, 1990). 3.4. Temporal variations in model parameters Model simulations were conducted at daily time steps, and input parameters were required for each simulation step. The dynamics of metrological variables directly controlled the transfer rate constants of horizontal diffusion and advection, atmospheric deposition, and wind resuspension in the air compartment. The temporal variation of the transfer rate constants of inter-media advections were mainly determined by hydrologic flows. In addition, the time-dependent changes in temperature, canopy growth, chemical emission, and background chemical concentrations were also considered in estimating the chemical and environmental properties. The temperature of the plant canopy, soil, and surface water was assumed equal to the ambient air temperature, under the conditions that the water temperature does not drop below 2 1C. The relationship defining changes of partition coefficients and half-life of a chemical with temperature suggested by Anderson and Hites (1996) and

ARTICLE IN PRESS Y. Luo et al. / Journal of Environmental Management 83 (2007) 44–55

Toose et al. (2004) was used in the study,    DE 1 1  kt ¼ kr exp , R Tr Tt

(16)

where kt is the temperature-dependent partition coefficient or half-life at temperature Tt (K), kr is the parameter at reference temperature Tr (298.15 K), DE (J mol1) is the activation energy for the chemical in a medium, and R is the gas law constant (8.134 J K1 mol1). At temperatures below 0 1C, precipitation becomes snow and soil water freezes. In this case, we assumed that the partition coefficients of a chemical between air and snow or ice can be approximated by extrapolating the air–water partition coefficients (Henry’s law constants) below the freezing point of water (Wania and Mackay, 1995). During a snow accumulating season (starting when the temperature drops below 0 1C, and ending when the temperature rises above 0 1C), wet deposition flux was formulated based on a scavenging ratio for snow. The current version of the model did not simulate canopy growth in the true sense, but merely calculated seasonal variations of leaf area index (LAI) for forests. By identifying the canopy as deciduous and coniferous forests, the seasonal variations of LAI described by Wania et al. (2000) were used in this study. The LAI of the coniferous canopy was assumed to increase slightly throughout spring and summer, and decline during the rest of the year. The deciduous LAI was assumed to increase quickly during spring, to be constant during summer, and then decrease rapidly in autumn. The LAI of a mixed canopy was calculated based on the area fractions of the forest types. The resultant LAI was used to calculate the timedependent parameters of (1) canopy volume, (2) litter fall rates, (3) atmospheric deposition velocities, and (4) interfacial areas in air–canopy diffusion. If chemical emissions and background concentrations from monitoring sites were not available for every time step, linear interpolation was employed to generate any missing data in the measurement time series. The model also allowed users to create temporal scenarios by specifying time series input data of emission rates and background concentrations. 3.5. Model capability and evaluation The outputs of the model provided comprehensive information of the contaminant distribution in a multimedia environment at watershed scale. The primary numerical results included (1) inter-media transport fluxes between any two adjacent compartments, (2) chemical inventory, fugacity, and concentration in all the compartments and sub-compartments, (3) chemical loss or gain rates by the mass exchange with the external environment for the compartments on the domain boundary, and (4) chemical loss rates by degradation. All the model outputs were computed simultaneously from the mass balance equations and reported as time series on a daily basis for

53

each compartment. This model was designed for use in a limited range of spatiotemporal scales, geographic conditions, and chemical classes. This model did not simulate chemical fate processes in the stratosphere and groundwater, and should not be applied for coastal and marine environments. The fugacity approach in the model was best suited to non-ionic organic chemicals. This model was no longer valid when chemical concentration exceeded the solubility limit in any compartment. As for time scale, the current model was intended for simulation over long time scales of several months to decades at daily steps. A full validation of the multimedia environmental fate model was impractical due to the lake of accurate and complete data of chemical emissions and concentrations. Comparisons of predicted and observed time trends or spatial patterns in terms of normalized concentrations, rather than absolute values, were conducted in this study to determine if the spatial and temporal variations in chemical distribution were well captured by the model even if the emission rates were inaccurate (Sweetman et al., 2002; Lee et al., 2004). The model was evaluated using trichloroethylene (TCE) as a test agent subjected to the field conditions in the Connecticut River Basin. Details of the model initialization and simulation are provided in the accompanying paper (Luo and Yang, 2006). 4. Conclusion A multimedia environmental fate model was developed for studying the temporal dynamics and spatial distributions of chemical pollutants at watershed scale. The core of the modeling approach was a set of fugacity-based mass balance equations to mimic the inter-media chemical fluxes in the multimedia environment, including air, canopy, soil, surface water, and sediment. Spatial and temporal variations in the chemical and environmental properties were considered in the simulation of chemical transport and transformation processes. With detailed environmental descriptions, this newly developed model was able to produce a better description of chemical transport and distribution in a multimedia environment compared to existing simulation tools. This model facilitated in-depth analysis of inter-media transports and multimedia system behaviors under dynamic conditions while preserving the requirements of modest data input and rapid scenario analysis. Therefore, this model is expected to provide a greater degree of realism relative to simpler environmental models, and lead to a more sophisticated and yet practical multimedia environmental simulation and analysis. Acknowledgement This study was conducted with financial support from the University of Connecticut Research Foundation, the Storrs Agricultural Experiment Station, and the Connecticut River Airshed–Watershed Consortium.

ARTICLE IN PRESS 54

Y. Luo et al. / Journal of Environmental Management 83 (2007) 44–55

References Ackers, P., White, W.R., 1975. Sediment transport: new approach and analysis. ASCE Journal of the Hydraulics Division 101, 621–625. Anderson, P.N., Hites, R.A., 1996. OH radical reactions: the major removal pathway for polychlorinated biphenyls from the atmosphere. Environmental Science and Technology 30, 1756–1763. Bai, Z., Feng, Y., 2003. GPS water vapor estimation using interpolated surface meteorological data from Australian Automatic Weather Stations. Journal of Global Positioning Systems 2 (2), 83–89. Bennett, D.H., Mckone, T.E., Matthies, M., Kastenberg, W.E., 1998. General formulation of characteristic travel distance for semivolative organic chemicals in a multimedia environment. Environmental Science and Technology 32, 4023–4030. Bicknell, B.R., Imhoff, J.C., Kittle, J.L.J., Donigian, A.S.J., Johanson, R.C., 1997. Hydrological Simulation Program—Fortran, User’s Manual for Version 11. EPA/600/R-97/080. National Exposure Research Laboratory, USEPA, Athens, GA. Brandmeyer, J.E., Karimi, H.A., 2000. Coupling methodologies for environmental models. Environmental Modelling and Software 15 (2000), 479–488. Charnock, T.W., Elgy, J., Hedges, P.D., 1996. Application of GIS linked environment models over a large area. In: Third International Conference/Workshop in Integrating GIS and Environmental Modeling, Santa Barbara, CA. ChemCAN, 2003. Level III fugacity model of regional fate of chemicals version 6.00, September 2003, developed by Mackay, D., Di Guardo, A., Paterson, S., Tam, D.D. Canadian Environmental Modelling Centre, Trent University, Peterborough, Ont., Canada. Cohen, Y., Cooter, E.J., 2002a. Multimedia environmental distribution of toxics (Mend-Tox). I: hybrid compartmental-spatial modeling framework. Practice Periodical of Hazardous, Toxic, and Radioactive Waste Management 6 (2), 70–86. Cohen, Y., Cooter, E.J., 2002b. Multimedia environmental distribution of toxics (Mend-Tox). II: software implementation and case studies. Practice Periodical of Hazardous, Toxic, and Radioactive Waste Management 6 (2), 87–101. Cohen, Y., Tsai, W., Chetty, S.L., Mayer, G.J., 1990. Dynamic partitioning of organic chemicals in regional environments: a multimedia screening-level modeling approach. Environmental Science and Technology 24 (10), 1549–1558. Coulibaly, L., 2000. Multimedia modeling of organic contaminants in the passaic river watershed in New Jersey, Thesis. Department of Environmental Engineering. New Jersey Institute of Technology. Coulibaly, L., Labib, M.E., Hazen, R., 1998. Multimedia modeling of environmental contamination: a case study of the New York harbor. In: The Fourth International Symposium on Environmental Geotechnology and Global Sustainable Development, Boston, MA. Coulibaly, L., Labib, M.E., Hazen, R., 2004a. A GIS-based multimedia watershed model: development and application. Chemosphere 55, 1067–1080. Coulibaly, L., Labib, M.E., Meegoda, J.N., 2004b. Multimedia model for analysis of contaminant releases in Passaic River Watershed. Practice Periodical of Hazardous, Toxic, and Radioactive Waste Management 8 (4), 220–227. Cowan, E.C., Mackay, D.M., Feijtel, T.C.J., van de Meent, A.D.G.D., Davies, J., Mackay, N., 1995. The Multi-Media Fate Model: A Vital Tool for Predicting the Fate of Chemicals. Society of Environmental Toxicology and Chemistry (SETAC) Press, Pensacola, FL. Cowherd, C., Muleski, G., Englehart, P.J., Gillette, D.A., 1985. Rapid assessment of exposure to particulate emissions from surface contamination sites. EPA/600/8-85/002, U.S. Environmental Protection Agency. Devillers, J., Bintein, S., Karcher, W., 1995. ChemFRANCE: a regional level III fugacity model applied to France. Chemosphere 30, 457–476. Diamond, M.L., Priemer, D.A., Law, N.L., 2001. Developing a Multimedia model of chemical dynamics in an urban area. Chemosphere 44, 1655–1667.

Duda, P., Kittle, J.L.J., Gray, M., Hummel, P., Dusenbury, R., 2001. WinHSPF version 2.0, an interactive Windows Interface to HSPF (WinHSPF). 68-C-98-010. Office of Water, USEPA, Washington, DC. ESRI, 2001. ArcGIS Extension Guides. ESRI Press. Hemond, H.F., Fechner-Levy, E.J., 2000. Chemical Fate and Transport in the Environment. Academic Press, New York. Hicks, B.B., Baldocchi, D.D., Meyers, T.P., Hosker Jr., R.P., Matt, D.R., 1987. A preliminary multiple resistance routine for deriving dry deposition velocities from measured quantities. Water, Air and Soil Pollution 36, 311–330. Jarvis, C.H., Stuart, N., 2001. A comparison among strategies for interpolating maximum and minimum daily air temperatures, part ii: the interaction between number of guiding variables and type of interpolation method. Journal of Applied Meteorology 40, 1075–1084. Jenson, S.K., Domingue, J.O., 1988. Extracting topographic structure from digital elevation data for geographic information system analysis. Photogrammetric Engineering and Remote Sensing 54 (11), 1593–1600. Jury, W.A., 1983. Behavior assessment model of trace organics in soil: I. Model description. Journal of Environmental Quality 12 (4), 558–564. King, A.W., 1990. Translating models across scales in the landscape. In: Turner, M.G., Gardner, R.H. (Eds.), Quantitative Methods in Landscape Ecology. Springer Publishers, New York, NY. Lee, Y., Lee, D.S., Kim, S.-K., Kim, Y.K., Kim, D.W., 2004. Use of the relative concentration to evaluate a multimedia model for PAHs in the absence of emission estimates. Environmental Science and Technology 2004 (38), 1079–1088. Lewis, W.K., Whitman, W.G., 1924. The two-film theory of gas absorption. Industrial Engineering Chemistry 16, 1215–1220. Luo, Y., Yang, X., 2004. A High-resolution multimedia model for chemical transport using GIS, ASAE Paper No. 041144, ASAE, Ottawa, Canada, August 1–4, 2004. Luo, Y., Yang, X., 2006. Dynamic modeling of chemical fate and transport in multimedia environment at watershed scale, Part II. Trichloroethylene test case. Journal of Environmental Management 83 (1), 56–65. Mackay, D., 2001. Multimedia Environmental Models: The Fugacity Approach, second ed. Lewis Publishers, Boca Raton, FL. Mackay, D., Yeun, A.T.K., 1983. Mass transfer coefficients correlations of organic solutes from water. Environmental Science and Technology 17, 211–233. Mackay, D., Joy, M., Paterson, S., 1983. A quantitative water, air, sediment interaction (QWASI) fugacity model for describing the fate of chemicals in lakes. Chemosphere 12, 981–997. MacLeod, M., Woodfine, D., Mackay, D., McKone, T., Bennett, D., Maddalena, R., 2001. BETR North America: a regionally segmented multimedia contaminant fate model for North America. Environmental Science Pollution Research 8 (3), 156–163. McKone, T.E., 1993. CalTOX, A Multi-media Total-Exposure Model for Hazardous Wastes Sites Part II: The Dynamic Multi-media Transport and Transformation Model. UCRL-CR_111456. The Lawrence Livermore National Laboratory, Livemore, CA. McKone, T.E., 1996. Alternative modeling approaches for contaminant fate in soils: uncertainty, variability, and reliability. Reliability Engineering and System Safety 54, 165–181. Meyers, T.P., Baldocchi, D.D., 1988. A comparison of models for deriving dry deposition fluxes of O3 and SO2 to a forest canopy. Tellus Series B 40 (1988), 270–284. Meyers, T.P., Flinkelstein, P., Clarke, J., Ellestad, T.G., Sims, P.F., 1998. A multilayer model for inferring dry deposition using standard meteorological measurements. Journal of Geophysical Research 103 (D17), 22,645–622,661. Nilson, R.H., Peterson, E.W., Lie, K.H., 1991. Atmospheric pumping: a mechanism causing vertical transport of contaminated gases through fractured permeable media. Journal of Geophysical Research 96 (B13), 21,933–921,948. Pennington, D.W., 2001. An evaluation of chemical persistence screening approaches. Chemosphere 44 (7), 1589–1601.

ARTICLE IN PRESS Y. Luo et al. / Journal of Environmental Management 83 (2007) 44–55 Strand, A., Hov, 1993. A two-dimensional zonally averaged transport model including convective motions and a new strategy for the numerical solution. Journal of Geophysical Research 98, 9023–9037. Sweetman, A.J., Cousins, I.T., Seth, R., Jones, K.C., Mackay, D., 2002. A dynamic level IV multimedia environmental model: application to the fate of Polychlorinated Biphenyls in the United Kingdom over a 60year period. Environmental Toxicology and Chemistry 21 (5), 930–940. Thibodeaux, L.J., 1995. Environmental Chemodynamics—Movement of Chemicals in Air, Water, and Soil. Wiley, New York. Toose, L., Woodfine, D.G., MacLeod, M., Mackay, D., Gouin, J., 2004. BETR-World: a geographically explicit model of chemical fate: application to transport of a-HCH to the Arctic. Environmental Science and Technology 128, 223–240. Wania, F., Mackay, D., 1995. A global distribution model for persistent organic chemicals. The Science of the Total Environment 160/161, 211–232.

55

Wania, F., Persson, J., Guardo, A.D., McLachlan, M.S., 2000. CoZMoPOP, a Fugacity-based Multi-compartmental Mass Balance Model of the Fate of Persistent Organic Pollutants in the Coastal Zone. WECC Report 1/2000. WECC Wania Environmental Chemists Corp., Toronto, Ont., Canada. Whitman, W.G., 1923. The two-film theory of gas absorption. Chemical and Metallurgical Engineering 29 (1923), 146–150. William, L., Mason, N., Chipman, V., 1997. Barometric Pumping with a Twist: VOC Containment and Remediation without Boreholes. Science and Engineering Associations, Inc., Santa Fe, NM. Woodfine, D.G., MacLeod, M., Mackay, D., Brimacombe, J.R., 2001. Development of continental scale multimedia contaminant fate models: integrating GIS. Environmental Science Pollution Research 8 (3), 164–172. Zhang, Q., Crittenden, J.C., Shonnard, D., Mihelcic, J.R., 2003. Development and evaluation of an environmental multimedia fate model CHEMGL for the Great Lakes region. Chemosphere 50, 1377–1397.