Evaluation of health risks for contaminated aquifers.

This review focuses on progress in the development of transport models for heterogeneous contaminated aquifers, the use of predicted contaminant concentrations in groundwater for risk assessment for heterogeneous human populations, and the evaluation of aquifer remediation technologies. Major limitations and areas for continuing research for all methods presented in this review are identified. ImagesFigure 2.


Introduction
Important objectives of an aquifer remediation technology are to reduce exposure concentrations of aqueous and nonaqueous phase contaminants below prescribed health standards for the least cost. Information required for meeting these objectives includes: determination of exposure concentrations in pumping or extraction wells prior to, during, and after remediation; health impacts of exposure on heterogeneous human populations; and evaluation of technology performance or effectiveness in achieving these goals. In the evaluation of remediation technologies, the variability of exposure doses and the variability of response in heterogeneous human populations are required. Figure 1 identifies the necessary information required, and the central role of the contaminant transport model for both risk assessment and for developing and evaluating reliable remediation technologies for contaminated aquifers. Predicted exposure doses from these models are required for risk assessments and they are required to evaluate the effectiveness of remediation technologies. This review focuses on progress in developing models for contaminant transport in the subsurface Groundwater quality data, variability of transport parameters, and aquifer characteristics data I Transport model development fore contaminated aquifer Modify initial values of transport parameters under conditions of uncertainty, in developing cancer risk assessment models for heterogeneous human populations, and in developing approaches for evaluating aquifer remediation technology performance.
Contaminant transport in the subsurface is transport through a porous structure that is neither uniform nor homogeneous. Subsurface contaminant transport models must incorporate and account for the impacts of spatial and compositional variabilities on contaminant transport processes and model parameters. In a given aquifer, the boundaries and flow channels of a groundwater aquifer cannot be determined with great precision. There are often zones that are completely saturated adjacent to zones that are partially saturated or unsaturated. The description of hydrodynamics in the subsurface is further complicated because the extent to which ground and surface waters are connected and the extent to which contaminants flow between ground and surface waters are difficult to characterize. The subsurface can also contain fractures, high permeability layers, discontinuous layers, impermeable lenses of soil and rock, and dislocations that can drastically alter local transport dynamics. In some situations, rapid convective transport channels can be located adjacent to slower convective-dispersive transport channels.  1. Development and uses of contaminant transport models. The development of contaminant transport models is central to developing and evaluating waste site remediation technologies and providing exposure assessment information for quantitative risk assessment for heterogeneous human populations.
Environmental Health Perspectives -Vol 105, Supplement -February 1997 I Structural and spatial heterogeneity of the subsurface affects the convective and dispersive transport parameters of groundwater and contaminant transport. Compositional variability affects interphase sorption-desorption processes of different chemical contaminants to organic matter and clay mineral constituents of the subsurface. Populations of microorganisms in saturated and unsaturated zones of the subsurface that transform contaminants into other chemical forms can vary in their capacities to transform groundwater contaminants by oxidative or reductive mechanisms. Transport in the unsaturated zones of the subsurface may require the inclusion of multiphase, multicontaminant, and energy transport equations to describe completely groundwater and contaminant transport. Taken together, all these factors contribute additional uncertainties to contaminant transport model predictions. A conceptual picture of the variety and complexity of the subsurface and the processes describing contaminant transport is shown in Figure 2.
For a contaminated aquifer, there is usually a paucity of data on the spatial variability of contaminant transport parameters, a major source of uncertainity in model predictions. Optimization methods are used to provide better estimates of these parameters because initial estimates of the statistical characteristics of these parameters must be inferred from information collected from a limited number of sparsely distributed observation wells or from textbook information. This process, referred to as the inverse problem, provides methods for developing better estimates of model parameters for a particular aquifer and makes it possible to predict spatial-temporal changes in contaminant concentration with less uncertainty.
As shown in Figure 1, the output from contaminant transport models provides important exposure information for estimating health risks in human populations. The responses of populations exposed to contaminants in groundwater have been studied and extensive assessments made of such populations' risks of cancer. The assessment of cancer risk as a result of exposure to environmental chemicals is a four-step process. The first step is hazard identification, in which a preliminary causative link is established between exposure and adverse health impacts. This is followed by an assessment to determine how and to what extent people are exposed. The third step is to establish a dose-response relationship that correlates exposure dose with cancer incidence. In this step, laboratory animal experiments are used to determine the cancer incidence as a function of exposure dose. Mathematical correlations determined in this step are used to estimate cancer risks at the lower exposure concentrations that are measured in contaminated aquifers. A safe dose or narrow range of safe doses are determined for a very low or negligible cancer incidence. Epidemiological studies on human populations also attempt to establish doseresponse relationships. Compared with the controlled environment of laboratory studies, epidemiology investigations must develop and use techniques to identify exposure, socioeconomic, and lifestyle factors from a very noisy group of candidate factors that may be primary contributors to disease incidence. It should be noted that the same process used in cancer risk assessment is being adapted and modified for risk assessment analyses of other health effects such as neurotoxicity, immunotoxicity, or reproductive toxicity that have environmental chemical exposure components.
Because of uncertainty in measuring and predicting exposure concentrations for contaminants in drinking water, cancer risk assessments that rely on these measurements and predictions are also uncertain and may not offer protection for sensitive populations. In heterogeneous human populations, there are variabilities in important physiological and pharmacological parameters for individual members that must be incorporated into the risk assessment process. Inclusion of these variabilities represents another source of uncertainty in risk assessment. In this review, heterogeneous human populations are defined as males and females of different ages and states of health as well as different measurements for physiological parameters (such as respiratory volume and heart rate) and pharmacological parameters (e.g., metabolic rates). Because of this variability, the same level of exposure to chemical contaminants in drinking water from contaminated aquifers will cause individuals to respond differently. Figure 3 identifies the interactions and intersections of individual physiologic, pharmacologic, and genetic characteristics, environmental factors, and lifestyle characteristics that contribute to the variability of responses observed in heterogeneous human populations as a result of exposure to contaminant chemicals via ingestion, inhalation, or dermal contact (1).
Once the extent of aquifer contamination has been estimated, selecting the most effective remediation technology or process and evaluating its performance at reducing health risks is the next major consideration. One method for evaluating the reduction of health risks and remediation costs for a specific remediation technology or process is the second moment method. For each level of contaminant reduction, a distribution function for exposure concentrations with a different mean and variance can be determined. In addition, a cost can be estimated for each level of contaminant concentration reduction. Variability of safe exposure concentrations for a heterogeneous human population can also be estimated from laboratory animal studies or from epidemiological studies. By incorporating the variability of risk into technology performance evaluations, the least costly, most effective technology design and operating conditions can be identified in order to ensure negligible health risks for all members of a heterogeneous human population. Figure 2 represents an attempt to identify the variety of subsurface structural and compositional features of an aquifer that must be accounted for in developing models to describe contaminant transport. Not all aquifers have the same level of complexity, but it is important and necessary to begin with a comprehensive representation that can then be modified as information on the transport characteristics of an aquifer is acquired. Thus, a mathematical model will be selected that best represents the transport dynamics of a specific aquifer.

Contaminant Transport Models
In the process of model selection and development, the first considerations are the structural and compositional characteristics of the subsurface aquifer and the chemical/biochemical characteristics of the contaminant. As indicated in Figure 2, an aquifer can be composed of porous solids only or it can be composed of porous solids with large irregular fractures or cracks. Within the porous solids, there can be small channels or pores, and the composition of these porous solids can vary spatially. Contaminants that flow very rapidly through one type of solid material will not necessarily flow very rapidly through another. Contaminant flow rates can be retarded because of sorption to different porous media materials. For both polar and nonpolar organic chemical contaminants, the same chemical usually is sorbed more strongly to organic matter than to day minerals. Sorptive strengths are a function of the molecular structures of contaminants and the solid components of the subsurface. Metal ions may form coordination complexes with organic matter and clay minerals. Along with structural/compositional variabilities that can affect sorption-desorption, different types of organic chemicals are degraded at different rates and by different mechanisms. Consequently, the products of these transformations have different binding or sorption characteristics with subsurface materials and surfaces.
Because of their vapor pressures, some organic chemicals are volatile and can migrate in the vapor phase through the Figure 3. Interactions of exposure, environmental factors, and internal physiologic and pharmacologic processes that affect disease-adapted from Hulka (1). Characterizing uncertainty in all of these processes is important for understanding individual susceptibility. subsurface and enter buildings through their substructures. Part of the subsurface can be partially saturated or unsaturated by water or another liquid phase and another part of the subsurface aquifer can be completely saturated by water or another fluid. The dynamics of contaminant transport depend on all these features and in the following sections, important components of the models for each of these situations are presented.
The traditional approach in contaminant transport model development is to designate a volumetric element in the subsurface as representative of the entire aquifer. In this conservation of mass approach, contaminants are transported as a result of a set of interacting transport phenomena that produce a deterministic model. Describing contaminant transport by a set of physical laws produces a set of partial differential equations that are solved with appropriate initial and boundary conditions. The major problem with this approach is that the volumetric element is usually made very small compared to the dimensions and scale of the aquifer. The heterogeneity of the subsurface makes it difficult to know exact values for model transport parameters at all locations and difficult to use this microscopic model for a small element of volume to approximate macroscopic contaminant transport behavior for the entire aquifer.
A stochastic contaminant transport model avoids a complete description of the physical-chemical laws governing contaminant transport at the microscopic level by using probability theory. In many instances, stochastic components and deterministic components (such as in Monte Carlo simulation) are combined in the contaminant transport model. If X(t) is a function that represents the position of a single contaminant particle at time t(x(0) = 0), then its trajectory or path through the subsurface is given as {X(t):t.0}. The trajectory is random and the probability law governing the random path of the contaminant is determined by the average (macroscopic) conditions governing flow. The major problem in this approach is the development and justification of the appropriate probability law for {X(t):t20}.
In most instances, the average macroscopic conditions governing transport in the stochastic model are the same set of physical-chemical laws describing transport at the microscopic level. For this reason, an appropriate approach for modeling contaminant transport in the subsurface is to start with the deterministic partial differential equation model. To incorportate subsurface heterogeneity, each transport process parameter is defined as a random variable. Some variables may be correlated; others may not. This is the approach that has been followed by a large number of researchers during the last 20 years. These same physical-chemical mass conserving laws appear to describe a wide variety of contaminant transport phenomena situations in the subsurface that includes single contaminants, multiple contaminants and multiple phases.
In this section, the description of models for different transport situations is followed by a description of techniques that have been developed to reduce parameter uncertainty. Even after the selection of the modeling approach has been made, it should be remembered that the predictive capability of the model always is limited and uncertain. These limitations are the result of availability or nonavailability of model transport parameters and the ability to model intra-and interphase contaminant transport processes or account for their effects on transport.

Contaminant Models for the Saturated Zone
Starting with the saturated zones of both confined and unconfined aquifers, describing contaminant transport requires the solution of the equation for saturated water transport and the convection-dispersionreaction equation for contaminant transport through porous media (2)(3)(4). For the saturated zone of the subsurface, contaminant transport is assumed to be isothermal and a heat balance equation is not required. Microscopic contaminant transport model equations are developed using a small representative element of volume in the subsurface. Stochastic versions of the transport equations provide a method to project this microscopic model into a macroscopic model that can be used to descibe transport for an entire aquifer. Stochastic models account for subsurface heterogeneities for which microscopic deterministic models cannot easily account (5-28), but they also introduce greater uncertainty into the predictions of contaminant concentrations in time and space.
If only isothermal contaminant transport in the saturated zone of the subsurface is considered, the general three-dimensional (3-D) physical law equation for saturated water transport through a representative small volumetric element in the porous structure of the subsurface is given as St is the storage coefficient of the aquifer, dimensionless h is hydraulic head, cm W(x,y,z,t) is the volume flux per unit area source term (positive for outflow and negative for inflow) cm/hr V is the del operator defined as (a/ax + a/ay + a/az), cm-l * denotes the dot product of vectors and tensors Ksat is the hydraulic conductivity tensor, related to fluid velocity, V; by Darcy's law, cm/hr In Equation 1, the notation V * (Ksate Vh) is referred to as the divergence of the vector field formed by the dot product (e) of hydraulic conductivity, Ksat, and the gradient of hydraulic head, V h. The divergence of the vector formed by Ksat V h has a simple meaning: it is the net rate of groundwater efflux per unit volume. Translating this notation into partial differential equation form gives the net rate of groundwater efflux per unit volume as The general 3-D physical-chemical law, convective-dispersive-reaction contaminant transport equation for this same element of volume for a single chemical is given as [3] where S is contaminant concentration sorbed to soil or other solid surfaces, cm3/cm3 Environmental Health Perspectives * Vol 105, Supplement I * February 1997 Cis contaminant concentration, g/cm3 PB is bulk density, g/cm3 1is porosity, dimensionless Dh is the hydrodynamic dispersion tensor, cm2/hr V is the fluid velocity vector in a porous media calculated from Darcy's law, cm/hr A is the reaction rate coefficient of order w for transformation of the contaminant by either chemical or biological processes, (cm3/g)'1/hr S, is contaminant concentration in the source or sink fluid, g/cm3 Q is volumetric flow rate, cm3/hr Vol is the volume of the volumetric element used in developing this model, cm3 For equilibrium sorption approximated by a linear Freundlich isotherm, the amount sorbed to solid components of the subsurface is approximated as S= KWC, where KW is the equilibrium sorption coefficient, cm3/g. As used here, KW is an overall or composite equilibrium sorption coefficient and is a function of sorption to all solid components, e.g., silica, clay minerals, and organic matter. Substituting for S, and assuming that all biological and chemical transformations are first order reaction processes, gives Rf dC = V. (D5 .VC) -Ve VC _AC SI Q [4] where the retardation factor, Rf (dimensionless), for equilibrium sorption is defined as If sorption is not an equilbrium process, it can be described by an interphase mass transfer process (12,19,22) as S = iCc, [6] where K is an overall mass transfer coefficient with units of (cm3/g)'>1 and o is the order of the mass transfer reaction of the chemical with the solid surfaces of the porous media. With this modification, the retardation factor is given as Rf =1 + [PB iCc-l] [7] Equation 3 is valid for a single chemical contaminant that is present in the aqueous phase in dilute concentrations. For most chemical contaminants, single chemical models are appropriate because the solubility of the chemical in the aqueous phase is very low (29). Thus, changes in time and space of contaminant concentrations are substantial. The first two terms on the righthand side of Equation 3 are associated with the transport dynamics of flow in porous media, namely hydrodynamic dispersion and convective transport. The last term accounts for first-order chemical and biological reaction rate processes that transform or convert contaminants into other chemicals in the subsurface. For dilute solutions, the assumption that reaction rate processes can be represented as first-order processes is reasonable because only the concentration of the contaminant changes substantially during degradation or transformation. For example, with hydrolysis reactions the water concentration remains virtually unchanged, a pseudo first-order degradation mechanism. For the element of volume in the saturated zone of the subsurface, a general form of the concentration profile for contaminant transport in the subsurface is given as a function of these transport parameters and subsurface characteristics as C(x,y,z,t) = F(Dh,Ksat,AIKC,E,PB) [8] However, some subsurface contaminant transport problems occur in conduit-type flow aquifers (30) that do not behave like flow-through porous media. For these aquifers and flow regimes, equivalent hydraulic characteristics must then be defined and used in Equation 8 to describe the contaminant concentration profile. Contaminant Transport Models for the Unsaturte Zone Contaminant transport in the unsaturated zone of the subsurface is important for describing transport of agricultural chemicals, spilled chemicals on the surface, chemicals leaking from underground storage tanks, chemicals placed in lagoons, and waste chemicals leaching from buried waste sites. For unsaturated transport, equations describing moisture, heat, and contaminant transport in the gaseous liquid, and solid phases of the subsurface must be included and solved simultaneously to predict spatial-time changes in contaminant concentrations (29,(31)(32)(33)(34)(35). If the contaminant is an unreactive organic solute or liquid with a low vapor pressure, vapor phase transport of the contaminant is usually negligible and is not induded in model development. However, if the contaminant is volatile or can be transformed into a gas, vapor phase transport of the contaminant must be included in model development along with the water vapor transport equation (36). Such behavior and changes in phase have been observed during the biologically mediated dechlorination of tetrachloroethylene and trichloroethylene to dichloroethylene isomers and then to vinyl chloride and finally to ethylene (37,38).
For contaminant transport in the unsaturated zone of an unconfined aquifer or water table aquifer, the element of volume required for description of all transport phenomena must include the atmosphere above the subsurface. In addition to atmospheric transport phenomena, transport of moisture, heat, and contaminant in both the water and the water vapor phases must also be included in the model. In Figure 4, the characteristics of this expanded element of volume and the processes and process laws that are included in a deterministic model are shown. Transport of heat, moisture, and contaminant are occurring simultaneously, but on the right-hand side of this diagram, the three transport phenomena have been separated to show the physical-chemical processes that are applicable. Using the subscripts w for water phase and v for vapor phase and letting 9 be the moisture content or degree of saturation, the moisture transport equation is written as the heat transport equation as at{ (1 ) Figure 4. Moisture, heat, and contaminant transport in the unsaturated zone. In the unsaturated zone, the transport of moisture, heat, and chemical contaminant occurs simultaneously. The three panels are used to illustrate the operable process laws for each type of transport phenomenon in the subsurface and to illustrate the importance of the interaction with atmospheric transport processes in driving transport in the unsaturated zone. and the contaminant transport equation with equilibrium sorption of the contaminant to soil surfaces as Pw [11] where V. is the water velocity vector, cm/hr V, is the water vapor velocity vector, cm/hr e is moisture content, < e, dimensionless psat (T) is saturated density of water vapor at temperature T, g/cm3 hRH is relative humidity, dimensionless percent cai,, c,, and c501jd are heat capacities for air, water, and solids, cal/g°K Pair, pu. and psolid are densities for air, water, and solids, g/cm3 H3, is the heat conduction through solid particles vector = -kLsolid * V T, cal/cm2/hr Hsw is the heat transfer through water phase vector = -* V T + cPV T, cal/cm2/hr Hs, is the latent heat + heat conducted by air vector = -4 Datmatort * VPwv -Air * V T, cal/cm2/hr AsoRid Aw, and 2aijr are heat conductivity vectors for solids, water, and air, cal/cm°K/hr 4 is latent heat of vaporization, cal/g Da,. is the molecular diffusion coefficient of water vapor in air, cm2/hr atort is the tortuosity of the porous media, dimensionless Q,s are source and sink terms, g/cm3/hr, e.g., Qij,x Qe.tp Q so H, is Henry's law constant, dimensionless KW is the equilibrium sorption coefficient for sorption from the water phase, cm3/g KV is the equilibrium sorption coefficient for sorption from the water vapor phase, cm3/g qf, is the contaminant mass flux vector in the water phase = -D,w .VCW + VwC, g/cm2/hr qcv is the contaminant mass flux vector in the water vapor phase = -D eV C,v + VVC,, g/cm2/hr Cw, is contaminant concentration in the water phase, g/cm3 C,, is contaminant concentration in the vapor phase and is related to water phase concentration by Henry's law, Cv = HIICW, g/cm3 Dcw is the dispersion coefficient for the water phase, cm2/hr Environmental Health Perspectives * Vol 105, Supplement * February 1997 132 D,, is the dispersion coefficient for the vapor phase, cm2/hr AW is the overall first order rate coefficient for contaminant degradation in the water phase, hr-1 AV is the overall first order rate coefficient for contaminant degradation in the vapor phase, hr-1 Using Henry's law to relate concentrations in the liquid phase to concentrations in the vapor phase, the complexity of the reaction rate term is given as AwoCW + (E -6)AvCv -AWOCW + (E -6)HCAvCw =((A+rW)0 + (awpwo + a, (Ee)Hepwv(T)))Cw [12] where J3w is the rate coefficient for first order biodegradation in the water phase, hr1 rw is the rate coefficient for first order chemical degradation in the water phase, hrl aW is the rate coefficient for hydrolysis in the water phase, hr1 av is the rate coefficient for hydrolysis in the vapor phase, hr' pWV( T) is the density of water vapor at temperature T, g/cm3 Of particular concern in this model for unsaturated flow are physical-chemical laws for describing the dynamics of competitive sorption of volatile contaminants and water vapor to soil surfaces. As moisture content increases toward saturation, concentrations of volatile organic chemicals such as p-xylene sorbed to clay minerals are greatly reduced (39). The sorption model used in this study underestimates the competitive sorption behavior of water and p-xylene to clay minerals, indicating that new models describing competitive sorptive behavior of contaminants in the saturated zone are required.

Multicomponent, Multiphase Contminant Transport Models
For multiple organic contaminant transport, transport equations are required for each constituent in each phase. Multiple phase contaminant transport is also referred to as transport of nonaqueous phase liquids (NAPLs). NAPLs that have densities greater than water tend to sink to the bottom of the saturated zone of the subsurface, near the bedrock or confining layer of an aquifer, and are referred to as dense nonaqueous phase liquids (DNAPLs). Chlorinated solvents such as trichloroethylene and tetrachloroethylene are examples of DNAPLs. Nonaqueous phase liquids that have densities that are less than water tend to float on top of groundwater in the saturated zone of the subsurface and are referred to as light nonaqueous phase liquids (LNAPLs). Lower molecular weight aromatic and aliphatic components of gasoline are examples of LNAPLs. The numerical simulation of downconing for NAPL free-phase recovery is an interesting problem (40).
In the unsaturated zone, the forms of the liquid transport equations for each liquid phase will be similar to the water-water vapor equations describing moisture transport in the unsaturated zone, Equation 9. In the saturated zone, liquid transport equations for each phase will have a form similar to that of Equation 1. In addition, equations will be required to describe interphase liquid transport or solubilization rates of one phase in another. The form of the set of equations needed for multiple contaminant, multiple phase transport is similar to that of the single contaminant transport equation given in Equation 3. The phases are gas (g), solid (s), water (w), and organic liquid (o). A general form of this model is given as [Abriola (36), Kirkner and Reeves (41), and Reeves and Kirkner (42) K-afla) where K,ap is the partition function of the it constituent between the a and P phases.
Important insight concerning competitive sorption processes can be obtained from the chromatography literature (43) and from studies of multiple phase, multiple contaminant flow in packed beds (44)(45)(46)(47)(48)(49). Clearly, describing multiple contaminant, multiple phase transport in the saturated and unsaturated zones of the subsurface requires more complicated models. A more complicated solution algorithm is required because more equations must be solved simultaneously. In addition, these more complicated models contain transport processes and parameters for which an understanding of the process laws and data on transport parameters is usually not available or is incomplete. As an example of how complicated the description of multiple chemical, multiple phase transport can be, continuous flow column experiments with mixtures of benzene and toluene suggest that model predictions for degradation of these two chemicals are very sensitive to initial active biomass concentrations, maximum specific substrate utilization rate, and metabolic rate parameters for bacteria (50).

Solving Conuminant Transport Models
For single contaminants in the saturated zone, for multiple phase single contaminants for the unsaturated zone, and for multiple phase multiple contaminants in the saturated and unsaturated zones, the water transport equation, the heat transport equation, and the contaminant transport equations represent a coupled transport system of partial differential equations. For isothermal transport in the saturated zone, the general method for obtaining changes in contaminant concentrations is to solve the water transport equation for spatial-time changes in water velocity with the appropriate initial and boundary conditions. For nonisothermal situations observed in the unsaturated zone, the heat and moisture transport equations must be solved simultaneously to obtain velocity fields as a function of temperature and hydraulic head (29). These velocities are then used to obtain spatial-time changes for contaminant concentrations defined by Equations 3, 11, and 13. Equations 1, 3,4, 9, 10, 11, Environmental Health Perspectives * Vol 105, Supplement 1 * February 1997 and 13 are nonlinear partial differential equations, and numerical computational methods will be required to solve them. The nonlinearity occurs because some transport parameters are functions of other transport parameters. An example is the hydrodynamic dispersion coefficient, which is a function of fluid velocity. In this situation, one approach has been to use a Picard iteration method coupled with a Newton-Raphson convergence algorithm. For nonlinear models, both finite difference and finite element approximations have been used to solve the governing partial differential equations (51).

Stochastic Models and the Inverse Problem
In the previous sections, transport models for a wide variety of contaminant transport situations in the subsurface have been introduced. In formulating these models, it is essential to recognize that there are significant sources of uncertainty that impose limits on predicted spatial-temporal changes in contaminant concentrations. In turn, these model limitations impose additional limitations on the reliability of calculated contaminant exposure concentrations used in health risk assessments. The first limitation and source of uncertainty is the scale of heterogeneity in the subsurface. Is it unbounded so that any value for a contaminant transport property is reasonable or is it bounded so that only specified sets of values for these parameters are reasonable? The second is the availability of adequate data on many of the parameters and transport variables in these models. The third source of uncertainty is the status of understanding of the physical-chemical-biological laws and processes occurring in a contaminated aquifer. Misunderstanding important contaminant transport processes can produce a model that fails to represent all the contaminant transport processes for a specific aquifer. In this situation erroneous estimates of contaminant exposure concentrations are produced.
One approach to confronting all these sources of uncertainty is to represent the transport parameters in groundwater and contaminant transport models as random space variables described by distribution functions and use Monte Carlo simulation to estimate spatial-temporal changes in concentrations (52)(53)(54)(55)(56). In Monte Carlo simulation, initial estimates of distribution functions for transport parameters such as hydraulic water conductivity, dispersivity, soil porosity, and the equilibrium sorption coefficient are either obtained directly or are estimated. Initial distribution functions are obtained from pumping test experiments, laboratory column or lysimeter test studies, or from sourcebooks (3). In many instances, lognormal or normal distributions are reasonable first approximations for these transport parameters (9). By using sensitivity analysis, transport parameters with the greatest effects on calculated contaminant concentration means and variances at specified locations are estimated.
In situations in which one or more transport parameters have negligible effects on calculated concentration means and variances, these parameters may be approximated by their mean values, thus reducing the complexity of the model and reducing computational time required to make model calculations.
In Monte Carlo simulation, values for each transport parameter are randomly selected from their initial distributions. These values are then used to calculate changes in groundwater flow velocities and contaminant concentrations for a spacetime grid that has been superimposed on the aquifer. Then a new set of values for model transport parameters are selected randomly from their distribution functions and a new set of contaminant concentrations is calculated for the same space-time grid. This process is repeated until a set of concentrations is generated at each grid point and at each time step. The sets of concentrations calculated in this manner are analyzed statistically to determine means and variances. Concentration means and variances are calculated at specified grid points such as monitoring or pumping wells and specified time intervals in the contaminated aquifer. To have enough calculated concentrations at each grid point for statistical analysis, a large number of iterations is selected in a Monte Carlo simulation. A more appropriate method for selecting the optimal number of iterations for Monte Carlo simulation is initially to select a very high iteration number. Then, after several hundred iterations, the contaminant concentration means and variances at several selected locations can be calculated. When changes in the mean and variance at each location are within prescribed tolerances, the simulation is terminated.
In most instances, however, the initial distribution function estimates for groundwater flow velocities and dispersion coefficients are uncertain because they are difficult to measure. As a result, predictions of contaminant concentrations at specified locations may not be very accurate. Finding optimal distribution functions for model transport parameters or solving the inverse problem for the contaminant transport equation has two major problems . The first is the nonlinearity of the dispersion coefficient as a function of groundwater flow velocity. The second problem relates to the first-distinguishing between convective and dispersive transport. For any of the methods for solving the inverse problem to be effective, groundwater quality data must be available and reliable. Calculated mean contaminant concentrations are compared with observed mean contaminant concentrations at monitoring or observation wells. Distribution functions for model transport variables are altered to produce a more reliable calibrated model for predicting spatial-temporal changes in contaminant concentrations in a contaminated aquifer. Commonly used methods to solve the inverse problem are kriging, cokriging, maximum likelihood, simulated annealing, and Bayesian updating. To illustrate Bayesian updating, let the uncertainty of a random transport variable A be defined by the initial or prior distribution P(A). By calculating contaminant concentrations with P(A) and comparing them with observed contaminant concentrations that have been calculated with a conditional probability distribution, P(IA IA), from Bayes' theorem, we have -P(IAI A)P(A)dA [14] where P(AIIA) is the updated or posterior distribution of the random transport variable A that predicts contaminant concentrations at a pumping well that are very close to actual measured concentrations. Figure 5 represents the steps in the development of transport models for a contaminated aquifer and how distribution functions for model transport variables are adjusted by Bayesian updating to provide better estimates of these variables for a particular contaminated aquifer (53)(54)(55)(56)(74)(75)(76)(77) [15] for noncarcinogens R (D) = 0, for D< Dt [16] where D£ is the threshold dose below which there is no observed adverse effect. In actuality, zero exposure doses cannot be achieved, and risk assessment is concerned with the determination of safe exposure doses or the distribution of safe exposure doses that are greater than zero but produce a negligible increase in cancer incidence over background. Presently, many of the same procedures for carcinogenic risk assessment are being adapted to and modified for immunotoxicity (91)(92)(93), reproductive and developmental toxicity (94)(95)(96)(97)(98), and neurotoxicity (99).

Models ofCarcinogenesis
For most cancer risk assessments, doseresponse relationships that correlate exposure dose to cancer incidence are determined from laboratory animal bioassay studies for individual chemicals. The standard design of these experiments uses three dose levels-a high or maximum tolerated dose, an intermediate dose, and a low dose, and 50 male and female test animals, usually rats and mice, for each exposure dose (90). The duration of the study is 2 years, which is about the average lifetime of these two rodent species. The different doses are administered daily to the dose groups, which makes it possible to calculate an average daily lifetime exposure dose for each exposure group. A fourth dose is the zero dose, which is maintained as a control group of animals to determine the background frequency of tumor formation in the animal species being used in the bioassay. At the end of the 2-year bioassay study, the number and types of tumors at each exposure dose are determined in both male and female animals. If animals die before the end of the 2-year test period, the time of death is recorded along with the types and number of tumors. Early death results provide important information on the kinetics of time to tumor formation, which is used in the development and calibration of models to predict carcinogenesis. Bioassay studies are conducted at concentrations that are high enough to produce tumors. For the much lower concentrations observed in the environment, models describing the carcinogenic process are required in order to extrapolate doseresponse relationships to the low-dose region. Current theories of carcinogenesis are based on the assumption that normal stem cells are transformed into intermediate cells and finally into malignant cells in a series of multiple steps or stages. An excellent review of the development of Environmental Health Perspectives * Vol 105, Supplement I * February 1997 multistage models of carcinogenesis is given by Whittemore and Keller (100). Many of the current models (100)(101)(102)(103)(104)(105)(106)(107)(108)(109)(110) are modifications and expansions of earlier work by Armitage and Doll (101).
In the multistage model, exposure to carcinogenic chemicals is assumed to enhance the transition rates of normal stem cells from one stage to the next. The simplest method of representing the transition rate from one stage to the next as a result of carcinogen exposure is a linear function, given as = a+fPiD [17] where Xi is the transition rate from one stage to the next, ai is a rate constant for spontaneous transitions, e.g., transitions that occur in the absence of exposure, and .Bi is the transition rate per unit dose as a result of exposure. The probability of a malignant cell being developed as a result of unrepaired mutational events occurring in normal cells is rare and can be represented by a Poisson process as R(D,t) = 1-exp[-ctkf(ai + PID)] [18] where c is a constant to be determined. For a fixed time or age, t, usually 2 years for a bioassay, the exponent becomes a polynominal in the dose rate D and the probability of cancer is given as R(D) = 1exp [-ko-klD-k2D2 _k3D3-. ..kD [19] For carcinogenic chemicals, values for the parameters in these multistage models of carcinogenesis are determined by maximum likehood methods applied to doseresponse data from the bioassay (88,90).
In determining carcinogenic risk in humans from animal experiments, a lifetime equivalent dose (ED) is determined in target tissues and organs in humans. Historically (90,111), a simple empirical scaling function has been used to relate bioassay exposure dose (D) to an ED in target tissues and organs, e.g., D/body weight (bw) = F(ED, surface area) as Lifetime average equivalent dose =ED= D [20] 2 (BWA)3 where D is lifetime daily average dose given to test animals, mg/kg/day BWA is the average test animal body weight, kg (about 0.03 kg for mice and 0.3 kg for rats) To determine safe exposure doses in humans, the safe ED in animals is converted to a safe human dose by inverting Equation 20 and substituting average human body weight (about 70 kg) for animal body weight. The cancer risk for a particular ED is determined by substituting ED for D in the multistage model of carcinogenesis (Equation 19).
There are several significant concerns with using models of carcinogenesis and empirical formulas to determine equivalent doses for estimating cancer risk and safe exposure concentrations for the very low exposure concentrations encountered in contaminated groundwater. As with contaminant transport models, there are uncertainties about the processes that transform normal cells into malignant cells and uncertainties in carcinogenesis model parameters. The parameters in a multistage model of carcinogenesis are determined in the high and medium dose ranges of the doseresponse curve. The parameterized model is then used to determine cancer risks for low exposure doses. In the low-dose region, it is not at all certain how cancer risk varies with dose. It makes a big difference in the determination of safe dose concentration for a given cancer risk if the shape of the curve is linear or nonlinear in the low-dose region of the dose-response curve, as shown in Figure  6. For the same cancer risk, the safe dose for a linear dose-response curve can be substantially less than that of a nonlinear dose-response curve that has a sublinear or concave shape in the low-dose region, as shown in Figure 6. Another way to interpret these two curves is that the linear model curve will provide a more conservative cancer risk estimate than the concave model curve. To overcome this problem, DNA adducts formed with carcinogenic chemicals or their metabolites may represent an important method to quantify the transition of normal cells into malignant cells and quantify cancer effects at low doses and provide insight about how this occurs. In addition, DNA adduct formation may represent a better method of correlating actual contaminant exposure concentrations observed in pumping wells with the processes that transform normal cells into malignant cells (112)(113)(114)(115)(116)(117)(118)(119)(120)(121)(122).

Physiologicaily Based Pharmacokinetic Models
The method for determining the equivalent or biologically effective dose in target organs and tissues does not directly include physiological and pharmacological processes in the whole animal. Empirical scale factors are used to relate exposure doses in laboratory animal studies to predict exposure doses in humans. In addition, risk is computed by using dose-response data for an average experimental animal to predict the cancer risk for an average human. By using the risk of cancer for an average human to predict risk of cancer for a heterogeneous human population, individual variability in EDs is de-emphasized and not adequately accounted for.
To incorporate biological information on carcinogenesis into the risk assessment process, it has been suggested [Andersen et al. (123)] that per-day average concentrations of contaminants or their metabolites are better choices for EDs. Equivalent or biologically effective doses of contaminants in drinking water and/or their metabolites in target tissues and organs can be determined using physiologically based pharmacokinetic (PBPK) models (123)(124)(125)(126)(127)(128). The key assumptions of the PBPK models are that organs and tissues can be represented as compartments or stirred tank reactors in which everything is perfectly mixed, and that a contaminant is transported through each compartment via the circulatory system. For each compartment, mass balances . Shapes of dose-response curves in the lowdose region. The shape of the dose-response curve in the low-dose region can greatly affect the estimate of safe doses for a negligible cancer incidence. If the curve is quadratic or concave downward in the lowdose region, safe doses for a specific cancer incidence will be much greater than if the curve is linear. describe contaminant transport, partitioning between blood, and compartmental tissues and transformation or metabolism of the contaminant. These phenomenological models are used to predict contaminant concentrations as functions of time in various organ systems and tissues. Figure 7 shows the compartmental characteristics of a generalized PBPK model that includes exposure via ingestion and inhalation. The general form of the mass balances for a contaminant in each compartment that includes not only flow processes but also transformation of the contaminant by biochemical/metabolic processes is given as In Figure 7, rapidly perfused tissue includes the spleen, brain, and kidney. Slowly perfused tissue is muscle; liver is the richly perfused tissue. Taking Equation 21 and applying it to rapidly perfused tissues, there are no biochemical transformation terms. In this case, the mass balance for contaminant in this compartment, CR, is given as dt VR C _ CR PR [22] [23] where CVR is venous blood concentration in rapidly perfused tissue, ng Ca,7 is arterial blood concentration, ng PR is the partition coefficient between blood and rapidly perfused tissue For the liver, the mass balance for contaminant concentration in this compartment, CL, for a saturable Michaelis-Menton metabolic mechanism is given as Figure 7. Flow diagram for PBPK model. The compartments represent regions of the body and can be discrete organs, such as kidney or liver, or widely distributed anatomical regions, such as muscle and fat. Mass balances describe transport by blood flow of contaminants into and out of each compartment. Rapidly perfused tissues include the spleen, brain, liver, and kidney. Slowly perfused tissues include muscle. Metabolism occurs predominantly in the liver but can also occur in the lung. Parameters in these models can be grouped as anatomicalblood flow rates, organ volumes; thermodynamic-equilibrium partition coefficients; physiologic-metabolic and clearance rates for contaminants; and transport-diffusion and sorption coefficients for contaminants.  [24] C = CL VZ -L [25] CVL is venous blood concentration in liver tissue, ng PL is the partition coefficient between blood and tissue in liver tissue V,"a, and Km are the Michaelis-Menton variables Equations similar to Equations 22 through 25 can be written for each compartment in order to form a set of differential equations that can be approximated by linear equations and solved by matrix algebra along with iterative convergence methods because of nonlinearities in the metabolic transformation parameters. If metabolism occurs in more than one organ system and by more than one mechanism, multiple organ sites and multiple risks must be considered and evaluated. PBPK models provide a biological basis for estimating EDs in target organs and tissues, but as they are structured by Equation 21, variability in the rates of transport and transformation of contaminants and the partition coefficients for different compartments is not included. In an alternative risk assessment method, Portier and Kaplan (129) used random variables for the transport variables of the PBPK model. With this randomized PBPK model, the distribution of safe doses was estimated for a minimal cancer risk of one in a million from exposure to methylene chloride, thus expanding on work carried out earlier by Andersen et al. (123). Of the variables in the PBPK models, sensitivity analysis indicated that variability in metabolic rate parameters appeared to affect changes in biologically effective doses of the metabolites of methylene chloride to the greatest extent. Metabolism of methylene chloride can occur by two pathways-a saturable Michaelis-Menton P450 mediated mixed function oxidase (MFO) mechanism and a first-order kinetic glutathione S-transferase (GST) mechanism. Animal experiments indicate that the GST pathway appears to be the carcinogenic mechanism (130,131). Variability in metabolic rate parameters for both pathways was introduced by increasing the variance of the distribution function for metabolic rate parameters. For a constant exposure concentration, safe dose estimates using a onestage model of carcinogenesis exhibited an increasing range of variability when there was increasing variability in the variance of the distribution function for metabolic rate parameters. Figure 8 is a frequency distribution plot of calculated safe exposure doses for a cancer risk of one in a million. As the variability of metabolic rate parameters increases, the uncertainty in estimating a safe dose for this negligible cancer probability increases, as exhibited by a flattening of the frequency versus dose curve. For a heterogeneous human population, the increasing uncertainty in safe doses may be more representative of the range of responses than would be observed, and may be a significant method of incorporating and representing dose-response variability in health risk assessment (108).

Evaluating Remediation Technology Performance
The performance of a remediation technology in reducing health risks for the least cost can be determined using the second moment method (132). Distribution functions for health risks in heterogeneous human populations and distribution  (129). In heterogeneous human populations, variability in metabolic rate parameters can affect the range of safe doses for a negligible cancer risk from exposure to methylene chloride. As variability in metabolic rate parameters increases, the range of safe exposure doses increases. Heterogeneous human populations are composed of males and females of different ages and with differing states of health and genetic characteristics.
functions for contaminant concentrations at pumping and/or monitoring wells are required for this method of reliability analysis. With the second moment method, the probability that contaminant concentrations will exceed a prescribed safe exposure dose or distribution of safe exposure doses is determined (132). The concentration of a contaminant (pg/liter) in the subsurface or at a pumping well is represented as a random variable (C) with a prescribed mean and variance. Distribution functions for contaminant concentrations can be either measured or calculated by a stochastic contaminant transport model that has been calibrated for a particular aquifer. As described previously, stochastic contaminant transport models start with Equations 1, 3, 4, 9, 10, 11, and 13 and randomize model transport parameters. Monte Carlo simulation is used to calculate contaminant concentration means and variances at all locations in the contaminated aquifer at different times or at specific locations such as pumping wells at different times. Safe exposure doses for a low cancer risk (pg/liter) in a heterogeneous human population are also represented as a random variable (R) with a prescribed mean and variance. For a negligible cancer risk, safe exposure doses should be greater than the concentration of contaminant within the contaminant plume or as observed in a pumping well. In this situation, R> C. Similarly, if R < C, the cancer risk will be greater than the prescribed low cancer risk. Mathematically, a remediation technology fails to reduce exposure concentrations when Pf = P(R < C) [26] where Pf is defined as the probability of failure of the technology to reduce contaminant concentrations below safe dose concentrations in a heterogeneous human population. Letting fr(r) and fc(c) define the distribution of safe doses and the distribution of contaminant concentrations, respectively, the probability of the remediation method not achieving reductions in contaminant concentrations that are below the prescribed low cancer risk is defined by the convolution integral as (132) 00 Pf = JF7(C)f,(C)dc 0 [27] where Fr(C) is the cumulative distribution function of fr(R) evaluated at C. Equation 27 can also be viewed as a method to Environmental Health Perspectives * Vol 105, Supplement * February 1997 incorporate exposure information into health risk assessment. The second moment method was used to calculate cancer risk contours for a contaminated aquifer at Hill Air Force Base, Utah (133). In this study, contaminant plume concentrations for trichloroethylene were incorporated into the risk assessment process to produce cancer risk contours for the contaminated aquifer.
Determining FT(C) can be accomplished by using dose-response relationships from laboratory animal studies or by using regression formulas from epidemiological studies that relate mortality or disease incidence to contaminant concentrations. There is limited epidemiological experience in developing regression formulas relating disease incidence to contaminant concentrations in groundwater because of lack of extensive or accurate monitoring data on contaminant exposure concentrations in pumping wells. However, there is extensive experience in developing regression formulas for relating cardiovascular and respiratory mortality to air pollutant concentrations that would suggest that the use of regression formulas is feasible for this performance evaluation method. Recent reviews on the development of regression formulas relating mortality to particulate air pollutant concentrations are given by Schwartz (134), Pope et al. (135,136), and Shumway et al. (137).

Conclusions
There are many similarities in the sources of uncertainty in estimating exposure concentrations with contaminant transport models and in estimating health risks with biologically effective exposure concentrations in target organs and tissues using PBPK models. Progress has been made on methods to reduce uncertainty in single contaminant concentration predictions for heterogeneous aquifers that are fairly uniform, e.g., without large fractures that drastically alter transport dynamics. For more complex aquifers and for unsaturated contaminant transport, multiple contaminant, multiple phase transport, more research is needed on understanding physical-chemical-biological processes and how they interact and affect contaminant transport. In addition, another important continuing research need is the characterization of the variability in response in heterogeneous human populations as well as the use of this information in risk assessment. Based on this review, several specific research needs are offered.
Reliable groundwater quality data are essential for the development of models describing contaminant transport in the subsurface and for risk assessment in heterogeneous human populations. Reliable and inexpensive methods are required to monitor more extensively contaminant concentrations in pumping wells and biological tissues. In this regard, immunoassays have been developed for monitoring contaminants in water, soil, food, and biological tissues (138)(139)(140)(141)(142)(143)(144)(145)(146)(147)(148)(149)(150)(151). Data from these assay systems, along with standardized chromatography/mass spectroscopy chemical analysis methods, make it possible to randomize the transport parameters in contaminant transport and PBPK models in a more reliable manner. Using water quality monitoring data from observation wells, initial distribution functions for transport variables can be adjusted to develop transport models that predict more accurately the transport behavior of a specific contaminated aquifer.
Many models of varying degrees of complexity for incorporating subsurface heterogeneity have been developed. Criteria are needed to identify conditions for the application of a particular model for a specific aquifer so that the uncertainty in estimating contaminant concentrations at pumping and monitoring wells can be minimized and to ensure that the most appropriate transport model is selected for that aquifer.
PBPK models that incorporate variability in physiologic and pharmacokinetic parameters for determining exposure or biologically effective doses in target organs and tissues may make it possible to represent the cancer risks for heterogeneous human populations in a more realistic manner. For a specified low or negligible cancer risk, a distribution of safe doses for a specific contaminant can be calculated rather than a single or very narrow range of safe doses. This makes the process of risk assessment much more complicated and uncertain but may represent an approach that incorporates individual susceptibility and diversity and may reduce the importance of the shape of the dose-response curve in the low-dose region. In addition to cancer risk assessment, risk assessment methods need to be developed for reproductive toxicity, neurological diseases, and disease resulting from immunological dysfunction.
The performance and cost of aquifer remediation technologies for reducing health risks can be evaluated more effectively if the variance and mean of the distribution function for human health risks as a result of exposure to aquifer contaminants can be determined with a greater degree of certainty than is currently available. The characterization of health risk variability for homogeneous and heterogeneous populations is very important information that is needed to develop and evaluate fully the effectiveness and performance of remediation technologies to reduce health risks for the least cost. In this regard, greater use of biological indicators of exposure and disease in epidemiological studies may be an important research activity to continue and to expand (152)(153)(154)(155)(156).