Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.
2 Modeling of Processes INTRODUCTION This chapter describes what models are and how they work. It begins by explaining the processes that control ground water flow and contaminant transport. To understand models, it is neces- sary to describe these processes by using certain mathematical equa- tions that quantitatively describe flow and transport. The mathe- matical aspects of modeling are critical. The precise language of mathematics provides one of the best ways to integrate and express knowledge about natural processes. By developing an awareness of the natural processes, the mathematics should be understandable. Also, where the process is not well understood, this awareness pro- vides an appreciation of the limits of the mathematics. Methods of solving the mathematical expressions are presented at the end of the chapter. Subsurface movement- whether of water, contaminants, or heat is affected by various processes. These processes can be related to three different modeling problems: ground water flow, multiphase flow (e.g., soil, water, and air; water and gasoline; or water and a dense nonaqueous-phase liquid (MAPLE, and the flow of contami- nants dissolved in ground water. 28
MODELING OF PROCESSES 29 Ground Water Flow Of these three problems, ground water flow is the simplest to characterize and understand In most cases, models need to consider only two ground water flow processes: flow in response to hydraulic potential gradients, and the loss or gain of water from sinks or sources, recharge, or pumping from wells. Hydraulic potential gradi- ents simply represent the difference in energy levels of water and are generated because precipitation that is added to a ground water sys- tem at high elevations has more potential energy or hydraulic head than water added at a lower elevation (Figure 2.~. The result of these potential differences is that water moves from areas of high po- tential to areas of lower potential. As rainfall or other recharge keeps supplying water to the flow system, ground water continues to flow. On a cross section, it is possible to represent the spatial variability in hydraulic potential existing along a flow system by using what are called equipotential lines (see Figure 2.1~. The equipotential lines are contours of hydraulic potential within some area of interest. In some simple situations, the direction of ground water flow is perpendicular to these equipotential lines, as shown in Figure 2.~. The actual distribution of hydraulic head observed for an area depends mainly on two factors, how much and where water is added and removed, and the hydraulic conductivity distribution that exists in the subsurface. Consider a few examples. Figure 2.2 illustrates the hydraulic head distribution for two different water table config- urations. The water table effectively represents the top boundary of the saturated ground water system, and its configuration reflects different recharge conditions. In both cases, the bottom and sides of the section are considered to be impermeable (no flow). With a linear water table and recharge mainly at the right end of the sys- tem, a relatively smooth regional flow system develops (see Figure 2.2a). The second water table, representing significant local areas of recharge and discharge at three locations, shows a much different flow pattern (see Figure 2.2b). Instead of a broad regional trend, several small, local flow systems have developed. Ground water flow patterns also depend on the hydraulic con- ductivity distribution. Figure 2.3 compares the pattern of ground water flow along a cross section where all properties except the hy- draulic conductivity for each layers are kept constant. Each of the two layers shown is defined in terms of a hydraulic conductivity in the horizontal direction (Kh) and in the vertical direction (Kit), with the ratio Kh/KV describing the degree of directional dependence in
30 o g In c~i ~6 !r ~ - U] ~88 _ __ C`l N~J I 11 r , aL~ 1 \~' Y 11 ~n I,i.~ ~ ~ ' 1 1 \ Ott'. 1 ~ ("J' 11, 1 $~. j In\ ~ tIt; \ ~ \\\ \ 'trim :~t i '~;- .... l U] I i o U.l I - ° C) ~] ° 01U ::~ ._ 4. o ._ _ ~ ~ r td ~ _ ~ · - 41) d 4o ud P4 "0 · - td ·_ 3 ~ o o d Le ~ oo 4= ~ ~ d ._ _ 3 ~ o ~ o bO 4. ~3 ~ o _ ~ ~Po ~ I L. L4 "._ _= 8 ~d ~ U] ·~= ~ ~,0 o O ~ C.) ~, · _ ~d to rd ~ 3 d o Le · e~ ~ d _` R r~ 4= I_ p, ~ ~q
MODELING OF PROCESSES 0.2S co m > 0.2S up a: 31 i ' ' 1 ---Equipotential Lines a ~ Ground Surtace , Water Table O _ OK = 1i ~ 1 ~'\ \~ W~/ _ - ~I _ _ 1 ~1, I I K = 10 1 1 ~ , , 1, , 1, , 1, , 1, , 1, , 1, , 1, , 1, , 1, , 1, , 1, . . \. S ~I ~ l ~I l I I I I I I I I ~i ~I b o O S RELATIVE BASIN LENGTH FIGURE 2.2 Dependence of the pattern of ground water flow on the recharge rate, as reflected by the configuration of the water table. All other parameters are the same in the two sections (from Freeze, 1969b). 0.2S ~ O id <i' > 0.2S UJ o a o ~I ~ ~ ~ ---Equipotential Lines C~rolint] Pilirf~r!- WAter Ts~hl- I ~11 | I K = lo ~I j 11 ~71 S 1 ~: o RELATIVE BASIN LENGTH S FIGURE 2.3 Dependence of the pattern of ground water flow on the hydraulic conductivity distribution. The only difference in the two diagrams is the pattern of geologic layering defined in terms of the relative hydraulic conductivities shown (from Freeze, 1969b).
32 GROUND WATER MODELS hydraulic conductivity. Examination of the two flow patterns shows how changes in the hydraulic conductivity distribution can change the character of ground water flow. Adding or removing water also can have a significant impact on the pattern of flow. The most important sources and sinks in a ground water flow system are pumping or injection wells (i.e., point sources/~nks). These are considered internal flows of water (fluxes). Other possibilities such as recharge or evaporation are most often considered as boundary fluxes. Pumping lowers the hydraulic potential at the well and in its immediate vicinity, creating what is known as a cone of depression. The result of decreasing hydraulic potential toward the well is the flow of water to the well. Injection does the opposite and results in flow away from a well. So far, only steady-state flow, or flow that does not change as a function of tune, has been discussed. Often, however, flow systems are transient, which means that hydraulic heads change with time, leading to variations in flow rates. For example, water leveh decline when a pumping well is first turned on, providing an early transient response. In many instances when sources of recharge are available, water levels will eventually stabilize, providing a new equilibrium or steady-state flow system. The most important feature of a transient flow system is the ability of water to be removed from or added to storage in individual layers. The parameter describing the water storage capabilities of a geologic unit is called the Unspecific storage." For transient flow problems, its value contributes to determining the distribution of hydraulic head at a given tune. Note that the smaller the specific storage, the faster the ground water system will seek a new equilibrium. Readers wishing a more detailed explanation of this parameter and aspects of ground water flow should see Freeze and Cherry (1979~. M~tiphase Flow Multiphase flow occurs when fluids other than water are mov- ing ~ the subsurface. These other fluids can include gases found in the soil zone or certain organic solvents that do not appreciably dissolve In water (i.e., immiscible liquids). Examples of fluids that are immiscible with water include many different manufactured or- ganic chemical such as the cleaning solvent trichioroethylene and preservatives such as creosote. Petroleum products such as crude oil, heating oil, gasoline, or jet fuel are also examples.
MODELING OF PROCESSES 33 The process causing all of these phases to flow ~ again movement in response to a potential gradient. Now, however, the situation is more complicated because the potential causing each fluid to move is not necessarily the same as that for water. Thus each fluid can be moving in a different direction and at a different rate. Another complexity Is that many characteristic parameters are no longer constant when several fluids are present together and competing for the same pore space. For example, the relative permeability of a geologic unit to a particular fluid like water will be small if the proportion of water present in a given volume of porous medium is small and will tend to increase as the amount of water increases. As discussed previously for water, a fluid's potential also depends on any sources or sinks that add or remove fluid. The same idea ap- plies to multifluid systems, except that now the number of processes increases because the effects of pumping/injection and evaporation (volatilization) affect each of the fluids present and, in addition, there can be transfers of mass between fluids. An example of this latter mechanism is that some portion of a gasoline spin might dissolve in water. To illustrate these concepts about the theory of multiphase sys- tems, consider two problems of particular interest to this report-the flow of water in the unsaturated zone and the migration of organic contaminants that are either more or less dense than water. When studying the problem of water movement in the presence of soil gas in the unsaturated zone, it is sometimes assumed that only the water moves. The only effect of the gas on water movement is the variabil- ity in the parameters caused by the presence of several fluids. For example, hydraulic conductivity varies as a function of the quantity of water in the pores. Figure 2.4 shows a relationship between hydraulic conductivity (K) and pressure head (fib). According to Freeze and Cherry (1979), pressure head is one component of the total energy water possesses at a point. Several features should be noted. As the pressure head becomes smaller (more negative), the soil becomes drier and the hy- draulic conductivity decreases. Much less water wall move through a dry soil than through a wet soil. Another feature ~ that if the soil is drying out there is one ¢-K relationship and if it is wet- ting there is another. Further, repeated wetting and drying cause the relationship to be defined by the scanning curves that join the wetting and drying curves at intermediate points (Figure 2.43. In
34 GROUND WATER MODELS 0.03 c) I > 0.02 it 8 0.01 Cot o I -300 -200 -100 0 100 200 300 Unsaturated Saturated Drying ~;- . I I ~ K o -Scanning Curves Wetting PRESSURE HEAD, ~ (cm H2O) FIGURE 2.4 Example of the relationships between pressure head and hydraulic conductivity for an unsaturated soil (modified from Freeze, 1971a). most multifluid systems, hydraulic conductivity and other parame- ters commonly exhibit this kind of "hysteretic" behavior, and yet for many applications, these types of site-specific data are not available. The progress of a wetting front moving into a dry soil can be described in terms of either potentials or volumetric water content, 0.14 - t _ ~ MOISTURE CONTENT FIGURE 2.5 The distribution of water in the unsaturated zone can be de- scribed in terms of pressure head and moisture content. Results presented for a combined saturated and unsaturated flow system illustrate how pressure head in particular is continuous across the water table (from Freeze, 1971b).
MODELING OF PROCESSES 35 defined as the ratio of the volume of water in the voids to the total volume of voids. Figure 2.5 illustrates how both are used to define a wetting zone near the top of the ground. The water table is clearly illustrated by the zero pressure contour and the total porosity contour (complete saturation). Given that moisture contents are easier to measure than potentials, the former are used more frequently to describe real systems. A more complicated case to consider is a flow involving an im- miscible fluid and water in the subsurface. Eventually, a distinction has to be made between a fluid that is less dense than water and one that is more dense. However, where an organic liquid is spilled on the ground surface, both fluids will move much the same way through the unsaturated zone (Figure 2.6a and b). The free organic liquid in a homogeneous medium moves vertically downward, leaving a residual trail of organic contaminants. Each pore through which the free organic liquid moves retains some of the contaminant (resicl- ual saturation) in a relatively immobile state. Thus, if the volume of spilled liquid ~ small and the unsaturated zone is relatively thick, no free liquid may reach the water table. Of course, free liquid may reach the water table over extended periods of time, and dissolved organic liquid may be conducted by water flow. It is when the free liquids begin to approach the top of the capillary fringe above the water table that the differences in density begin to affect transport. The capillary fringe is a zone above the water table where the pores are completely saturated with water but the pressure heads are less than atmospheric (Freeze and Cherry, 1979~. A contaminant that is lighter than water will mound and spread, following the dip of the water table (Figure 2.6a). A fluid that is heavier than water will spread slightly and keep moving downward. This fluid will ultimately mound on the bottom of the aquifer or on a low-permeability bed within the aquifer and move in whatever direction the unit is dipping (Figure 2.6b). Thus water and the organic liquid need not move in the same direction. To understand the details of multicomponent flow, it is essential to study the concepts of wettability, imbibition and drainage, and relative permeability. A discussion of these topics is, however, beyond the scope of this report. Readers can refer to Bear (1972) and Greenkorn (1983) for en overview of the basic theory. The key point to remember is that, as in the case of water, the permeability of the material through which these fluids are moving plays a major role in controlling the direction and rate of flow. In the case of multiphase
36 GROUND WATER MODELS a Ground Surface . _ _ Oil Phase -Oil Components Unsaturated Zone Capillary Fringe Dissolved in Water~ --__ Saturated Zone ~/////'~0~: b Ground Surface _v _ ~- ~ , ~ ~ ~ 1: .~#~. 1 t..~ :! ~ CHC Phase Unsaturated Zone Capillary Fringe K2 K1 CHC Dissolved in Water Saturated _ _ _ _ _ _ _ _ ' Zone l _ I_ ~////////////////////////////////////////i/////~, FIGURE 2.6 The flow of a nonaqueous-phase liquid that is (a) less dense than water (oil) and (b) more dense than water (chlorohydrocarbon, CHC) in the unsaturated and saturated zones. In both cases the contaminants are also transported as dissolved compounds in the ground water (from Schwille, 1984~. systems, a relative permeability is defined for each fluid with values ranging between zero and one as the relative abundance of each fluid (i.e., saturation) changes. The key point here is that as more fluids are introduced to the pore space, more of the pore space is devoted to the relatively immobile state of each fluid and therefore less pore space is devoted to liquid flow.
MODELING OF PROCESSES 37 The distribution of an NAP L in the subsurface is described quan- titatively in terms of the relative saturation of the NAPL, which is given by the ratio of the volume of the NAPL to the total pore vol- ume. In other words, it describes what proportion of the pore volume is filled with the NAPL. A relative saturation can be defined for each one of the organic liquids and water. This kind of description is generally not used in field settings because of the detailed study that is necessary. Instead, presence/absence indications are used, as illus- trated in Figure 2.6. The results of computer simulations normally characterize relative fluid saturations. Expressing the distribution of fluids in terms of relative saturation is analogous to expressing the moisture content in terms of unsaturated flow. Dissolved Contanunant Passport One of the reasons why problems involving dissolved contam- inants are so difficult to model is the number and complexity of controlling processes. The processes can be divided into two groups: (1) those responsible for material fluxes and (2) sources or sinks for the material. For the problem of contaminant migration these are the mass transport and mass transfer processes, respectively (Table 2.1~. A brief discussion of each of the processes listed in Table 2.1 follows, with a general assessment of its impact on contaminant transport. Advection Advection is the primary process responsible for contaminant migration in the subsurface. Mass is transported simply because the ground water in which it is dissolved is moving in a flow system. In most cases, it can be assumed that dissolved mass is transported in the same direction and with the same velocity as the ground water itself. For example, given the conditions of flow described by the equipotential lines and flowlines of Figure 2.7a, it is a simple matter to define the plume of dissolved contaminants in terms of the streamtubes that pass through the source. A streamtube is defined an the area between two adjacent Bowlines. When Bowlines are equally spaced, the discharge of water through each is the same (Freeze and Cherry, 1979~. This simple approach assumes that the density of the contaminated fluid is about the same as that of the ground water. The mean velocity of contaminant migration can also be assumed to be the same as the mean ground water velocity (or seepage velocity).
38 t-TROUND WATER MODELS TABLE 2.1 A Summary of the Processes Important in Dissolved Contaminant Transport and Their Impact on Contaminant Spreading Process Definition Impact on Transport Mass transport 1. Advection Movement of mass as a Most important way of consequence of ground transporting mass away water flow. from source. 2. Diffusion Mass spreading due to An attenuation mechanism molecular diffusion in of second order in most response to concentration flow systems where gradients. advection and dispersion dominate. 3. Dispersion Fluid mixing due to effects An attenuation mechanism of unresolved hetero- that reduces contaminant geneities in the per- concentration in the meability distribution. plume. However, it spreads to a greater extent than predicted by advection alone. Chemical mass transfer 4. Radioactive decay 5. Sorption 6. Dissolution/ . . . precipitation Irreversible decline in the activity of a radionuclide through a nuclear reaction. Partitioning of a contaminant between the ground water and mineral or organic solids in the aquifer. The process of adding contaminants to, or removing them from, solution by reactions dissolving or creating various solids. 7. Acid/base Reactions involving a reactions transfer of protons (H+). An important mechanism for contaminant attenuation when the half-life for decay is comparable to or less than the residence time of the flow system. Also adds complexity in production of daughter products. An important mechanism that reduces the rate at which the contaminants are apparently moving. Makes it more difficult to remove contamination at a site. Contaminant precipitation is an important attenuation mechanism that can control the concentration of contaminant in solution. Solution concentration is mainly controlled either at the source or at a reaction front. Mainly an indirect control on contaminant transport by controlling the pH of ground water.
MODELING OF PROCESSES TABLE 2. 1 Continued 39 Process 8. Complexation Definition Impact on Transport An important mechanism resulting in increased solubility of metals in ground water, if . . ac sorption IS not enhanced. Major ion complexation will increase the quantity of a solid dissolved in solution. Often hydrolysis/substitution reactions make an organic compound more suscep- tible to biodegradation and more soluble. Combination of cations and anions to form a more complex ion. 9. Hydrolysis/ substitution 10. Redox reactions (biodegradation) Reaction of a halogenated organic compound with water or a component ion of water (hydrolysis) or with another anion (substitution) . Reactions that involve a transfer of electrons and include elements with more than one oxidation state. Biologically mediated mass transfer 11. Biological Reactions involving the transformations degradation of organic compounds, whose rate is controlled by the abun- dance of the microorgan- isms and redox conditions. An extremely important family of reactions in retarding contaminant spread through the precipitation of metals. Important mechanism for contaminant reduction, but can lead to undesirable daughter products. The close relationship between advective transport and ground water flow means that the factors considered for flow, the location and quantity of the inflow and outflow to the flow system, the hydraulic conductivity distribution, and the presence of pumping/injection wells also play a major role in determining where contarn~nants mi- grate. Indeed, the process of advection is often so dominant that the mean velocity predicted by flow models can be used to estimate patterns of contaminant transport with surprising accuracy. Diffusion Diffusion is an important process that results in mass mixing. Diffusion is mass transport in response to a concentration gradi- ent. Thus contaminants present in a plume will diffuse away from
40 25 10 o GROUND WATER MODELS a >, b 25 20 15 10 5 o o 40 80 120 DISTANCE (m) FIGURE 2.7 Plume produced (a) by advection alone and (b) by advection and dispersion (adapted from Frind, 1987~. the plume in all directions in response to concentration gradients. Although this process occurs in most contaminant problems, its overall contribution to the spreading of contaminants is usually neg- ligible. There are situations, mainly in fractured rock settings and low-permeability units, where diffusive mass transport is of primary importance and needs to be considered. Dispersion Dispersion refers generally to phenomena that cause fluid mix- ing. Dispersion is more accurately described as the apparent mixing due to unresolved advective movement at scales finer than captured by the mean advection model. Essentially, dispersion produces a mixing zone between the contaminated water and the native ground water. This effect can be illustrated by considering a plume devel- oped due to advection alone (Figure 2.9a) and modifying it to also include dispersion (Figure 2.7b). A comparison of parts a and b in Figure 2.7 shows that dispersion has expanded the plume size be- yond that expected due to advection alone. Contaminants spread into adjacent streamtubes and farther down the streamtube where the contaminants are migrating. The overall plume becomes larger,
MODELING OF PROCESSES 41 and, in general, the concentration is less than was the case with advection alone. Dispersion in a direction perpendicular to the mean direction of ground water flow is termed transverse dispersion, while disper- sion parallel to the mean direction of flow is termed longitudinal dispersion. There are actually two directions of transverse spreading (Figure 2.~. These different components of dispersion usually need to be considered separately in models because spreading upward or downward is often considerably less than spreading in horizontal or subhorizontal planes. Hydrodynamic mixing occurs as a consequence of nonidealities at various scales that result in local variability in velocity around some mean Precocity. For example, at the scale of pores (Figure 2.9a) this variability may be caused by velocity variations within a pore or by subtle changes in the flow network that cause the mass to spread out or finger into adjacent pores of the pore network. At the macroscopic scale, the variability can be due to heterogeneities in the hydraulic conductivity distribution (Smith and Schwartz, 1980) of the kind shown in Figure 2.9. In terms of the relative magnitude or mixing due to hydrodynamic dispersion and diffusion, the former is by far the more significant. Readers interested in learning more about dispersion should refer to Schwartz (1975, 1977), Anderson (1979, 1984), Tennessee Valley Authority (1985), Mackay et al. (1986), Freyberg (1986), and Sudicky (1986~. Transverse Dispersion Transverse Dispersion T . Longitudinal Dispersion FIGURE 2.8 Idealized pattern of plume spreading in three dimensions is characterized by a longitudinal and two transverse dispersion components.
42 GROUND WATER MODELS a Microscale Dispersion _ .~- ~ Am. ::. D..-0 . . . do, O. Mixing in Individual Pores (from Freeze and Cherry, 1979) b Macroscale Dispersion Mixing in a Pore Network (from Cherry et al., 1975) 5 m .~` 5 m Mixing Caused by Variability in Hydraulic Conductivity (from Schwartz, 1984) [> [> FIGURE 2.9 (a) Microscale and (b) macroscale variability contributing to the development of dispersion. Radioactive Decay Radioactive decay, the transformation of one element into an- other through the loss of atomic particles from the element's nucleus, is a process that has been thoroughly characterized and is well under- stood. Radioactive decay leads to the loss of the original radioactive isotope from ground water over a period of time, but daughter prod- ucts are produced that may also be of environmental concern. A
MODELING OF PROCESSES 43 simple rate law can be applied to the decay of any radioactive iso- tope, and it has been included in various transport models for many years. It describes an exponential decrease with time in the concen- tration of the dissolved radioactive component. The best-known example of radioactive decay in ground water is that of tritium (3H). Tritium is produced naturally by interaction of cosmic rays (various nuclear particles coming in from outer space) with gases in the upper atmosphere. Consequently, trace amounts of tritium are found in all natural waters. During the 1950s, large amounts of new tritium were injected into the atmosphere as a result of the testing of fusion bombs. The anomalously high concentrations of bomb-produced tritium led to much interesting work in dating and tracing the patterns of flow of ground waters. Tritium, with a half- life of 12.5 yr, decays to stable helium (3 He) by emission of a beta particle. Because of its short half-life, tritium produced by the testing of atomic weapons in the atmosphere is gradually disappearing from natural waters. Of environmental importance are radioactive species that may be inadvertently released into ground water from such activities as mining, milling, and storage of wastes. In particular, concern ex- ists about the escape and potential hazards of radium, uranium, and lead in ground water adjacent to uranium mills and processing plants, and about the leaching and movement of radioactive isotopes (including isotopes of uranium, plutonium, cesium, neptunium, eu- ropium, iodine, selenium, and others) away from geologic repositories for high-level radioactive wastes from commercial power plants and defense installations (Bates and Seefel~t, 1987; Fried, 1975~. Because the radioactive isotopes of concern in the environment undergo other chemical reactions in addition to radioactive decay, many years of research will be required before their behavior can be modeled with confidence, even though radioactive decay is well understood. When the half-life for radioactive decay is of the same magni- tude or smaller than the residence time of the contaminant in the subsurface, decay significantly affects contaminant migration. This is illustrated in Figure 2.10. Figure 2.10a illustrates what a hypo- thetica] plume might look like if advection and dispersion were the only controlling processes. The plume is much larger than the one in Figure 2.10b, where it is assumed that radioactive decay is also operative. As chemical and biological processes are discussed, it will become apparent that in general their effect is to attenuate the spread of
44 15. _ 10. _ 5. _ LL.' O _ O. Hi: G Ul > 15. ~ GROUND WATER MODELS a b 10. ~ 40 80 120 5 ~(~-~.1 1 1 1 1 o 40 80 120 HORIZONTAL DISTANCE (m) FIGURE 2.10 Many of the geochemical processes like radioactive decay and sorption attenuate the spread of contaminants. Compared on the figure are map views of half-plumes (a) without attenuation and (b) with attenuation (based on Grind, 1987~. contaminants relative to that caused by advection and dispersion alone. It is for this reason that so much emphasis in modeling has been placed on accounting for these processes to the fullest extent possible. Sorption Adsorption reactions remove contaminants from ground water and add them to the surfaces of minerals or the solid organic carbon of the unit through which the contaminants are moving. The term sorption is a general one that includes adsorption (attraction to a surface), absorption (incorporation into the interior of a solid), ion exchange (adsorption, with a charge-for-charge replacement of the ionic species on a surface by other ionic species in solution), and Resorption (the opposite of each of the above adsorption reactions). Sorption will affect virtually all dissolved species in ground water to some degree. Sorption is such a complex process that it is not really possible within the scope of this review to provide a complete appreciation of what causes contaminants to move from solution onto solids. Metal ions are sorbed primarily because of the positive charges they carry or chemical reactions that bind them to the surface. Clay minerals in particular have large surface areas carrying an overall negative
MODELING OF PROCESSES 45 charge. This surface charge is balanced by positively charged ions that are attracted to the surface. Contaminant ions in ground water can in many cases preferentially replace these positively charged species. Other surfaces (e.g., metal oxides or metal oxyhydroxides) are reactive in the sense that a metal or certain metal-conta~ning compounds can be chemically bound to the surface. These surfaces are particularly interesting because the reactivity of the surface, or the ability to sorb contaminants, is controlled in part by the pH of the ground water. Another major class of contaminants, noncharged organic mole- cules, sorb mainly onto solid organic material. The force that drives the exchange in this case is the hydrophobic (water hating) character of some organic compounds. For such chemicals, sorption onto solids increases as the solubility in water decreases. Thus the more water soluble the contaminant, the less likely it is to be sorbed. When an organic phase is present, either solid or liquid, organic contaminants prefer to reside in that phase. Many sorption reactions discussed are completely or partially reversible. In other words, if the concentration of the contaminant in the ground water decreases, Resorption will occur to maintain an equilibrium between the contaminant in solution and that sorbed on the solids. Thus sorption does not permanently remove a contami- nant from solution but instead only stores it. A number of empirical or semiempirical methods have been de- veloped for describing sorption equilibrium. The so-called distribu- tion coefficient (K&) is the most simple of these; Ka is defined as the concentration of a given contaminant sorbed on the solid phase (commonly in micrograms per gram) divided by the concentration of the same contaminant in solution (in micrograms per milliliter), with the resulting units being milliliters per gram. A large Ka value indicates strong sorption, or that the compound distributes itself primarily onto or into the solid phase. A small Ka indicates that the compound stays mainly in the water phase. Therefore Ka serves as a qualitative guide to the relative tendency toward sorption of various dissolved species in a given solution. The usefulness of Ka is diminished somewhat by the fact that it may vary as a function of concentration, ionic strength, competing ions, and other factors. In soils, sediments, and some aquifers, solid organic matter is the primary solid material onto which organic compounds sorb. For many organic compounds, empirical relationships can be derived to
46 GROUND WATER MODELS predict the Ka value as a function of the amount of organic ma- terial in the solid phase and a measure of the organic compound's hydrophobic nature, most usually its octanol/water partition coef- ficient. A large octanol/water partition coefficient signifies a highly hydrophobic compound, which will have a large Ka. The overall effect of sorption is to retard or delay the spread of contaminants. This effect is not unlike that shown in Figure 2.10. When sorption occurs, the rate at which the contaminant appears to move is lower than would be the case for an unretarded or neutral tracer. This behavior helps to reduce the spread of contaminants but also makes it more difficult to remove contamination from the ground; that is, it tends to increase the time required to remediate to a cleanup level. Precipitation and Dissolution Dissolved contaminants can be either lost from solution or brought into solution by the processes of precipitation and disso- {ution of a solid phase. RunnelIs (1976) gives exa~nples of precipita- tion of dissolved contaminants in ground water caused by reactions with other dissolved species, hydrolysis, and reduction or oxidation. Examples of contaminants that could be reduced to lower concentra- tions in ground water through the formation of precipitates include arsenic (by reaction with iron, aluminum, or calcium), lead (by reac- tion with sulfide or carbonate), and silver (by reaction with sulfide or chIoride). Hydrolysis can lead to the precipitation of iron, man- ganese, copper, chromium, and zinc contaminants. Oxidation or reduction could favor the precipitation of chromium, arsenic, and selenium. The precipitation of dissolved contamination plays an important role in contaminant attenuation. Although this process is not well described in case studies, theoretical work shows that it will atten- uate the spread of contaminants by removing mass from solution as saturation is exceeded. Unlike the sorption processes that also partition mass between the solid and solution, these reactions are less reversible. For example, metals precipitating as metal-sulfides are virtually immobilized for as long a time as the general chemical environment remains constant. Contaminant dissolution is an important reaction that can occur at a source to initially bring contaminants into solution. However, for some minerals that dissolve relatively rapidly, it is possible for
MODELING OF PROCESSES 47 contaminants initially immobilized by precipitation to be remobilized as the plume moves further down the flow system. A natural analog of the repeated precipitation and remobilization of metals is a uranium roll-front deposit (Galloway and Hobby, 1983~. Acid/Base Reactions Reactions involving the gain or loss of the hydrogen ion (H+) are called acid/base reactions. Acids are chemical species that give up or donate a H+ ion, while bases are species that accept a H+ ion. Many potential contaminants are susceptible to change in speciation because of changes in pH. For example, under oxidizing conditions, dissolved arsenic should be present in normal ground water (with a pH of 7 to 8) in the form of HAsO42-. However, if the pH of the water is lower than about 6, the dominant form of oxidized arsenic is H2AsO4- or, at very low pH, H3AsO4-. Depending on the number of protons attached to the arsenate, the chemical behavior of the arsenic in solution may be quite different. For example, the sorption behavior of H3AsO4 is quite different from the sorption behavior of H2AsO4-. There are cases where the chemical reaction being considered does not include the contaminant. In terms of understanding trans- port, this means that in some cases the reactions that control im- portant geochemical parameters of a system (e.g., pH) must be con- sidered even though the compounds or ions involved in the specific reaction are not contaminants. For this reason, sophisticated mass transport models often need to include reactions related to the CO2- water system, one of the dominant controls of the pH of ground water. Comp le x at ion The process of complexation is the combination of simple cations and anions into more complex aqueous species. According to More! (1983), complexes can be classified as ion pairs of major constituents, inorganic complexes of rare metals, and organic complexes. Following are examples of reactions forming each of these complexes: Ca2+ + SO42- = CaSO4, Cu2+ + H2O = CuOH+ + H+ Cu2+ + Y~ = CuY+
48 GROUND WATER MODELS where Y~ is an organic species such as glycine. Conventionally, complexation reactions are modeled as equilibrium processes where equations like those above are characterized in terms of mass law ex- pressions and equilibrium constants. This relatively straightforward treatment enables the concentration of the individual complexes in water to be easily calculated. In terms of mass transport, complexation reactions are impor- tant mainly because of the role they play in increasing the mobility of metals. For example, over the range of pH common to most natural ground waters, metals (present as ions) will occur at relatively low concentration because of the solubility constraints provided by solid phases including metal-hydroxides, carbonates, or sulfides. However, when metals complex to a significant extent, the total quantity of a particular metal dissolved in water can be much larger than simply the concentration of the metal ion itself. Overall then, complexation can enhance the quantity of metals being moved in a contaminant plume. It is for this reason that most models involving metal trans- port need to consider complexation reactions. Another instance where complexation needs to be considered is in the sorption of metals on surfaces whose charge changes as a function of pH. Examples of such surfaces include kaolinite, metal oxides, and metal oxyhydroxides. For such solids, the sorptive behavior changes depending upon which metal species (ion and complexes) are present in the solution. Additional complexity arises from the fact that changing the composition of the water (e.g., pH) also changes the concentration of various metal species in solution. Thus, in situations where this type of sorption can occur, the metal complexation must be included to fully characterize the sorption reactions. Hydrolysis/Substitution Hydrolysis and substitution are abiotic transformation reactions that affect organic contaminants in ground water. The term hydrol- ysis refers specifically to substitution reactions involving water or a component of water, for example: RX+H2O ROH+HX, where R refers to the main part of the organic molecule and X is a halogen (e.g., CI~, Br~) (Jackson et al., 1985~. However, not all substitution reactions involve water. For example HS- can react and
MODELING OF PROCESSES 49 substitute for a halogen (e.g., Br~) in a reaction of the following kind (Jackson et al., ~ 985~: RCHX + HS- _ RCHSH + X-. Kinetic rate laws for hydrolysis/substitution reactions can be com- plex. However, in some cases, they can be approximated as first-order reactions or, in other words, reactions like radioactive decay that can be described simply in terms of a half-life. The reason transformation reactions are important is that products are often more susceptible to biodegradation and more soluble. Redox Reactions Any element that can have different valences is potentially sub- ject to transformation via oxidation and reduction reactions, known in short as "redox" reactions. Included are abundant species, such as nitrogen, sulfur, carbon, and phosphorus, and minor trace species, such as iron, manganese, uranium, selenium, copper, and arsenic. Redox reactions involve the movement of electrons from one species to another. Redox reactions change the speciation of the dissolved elements, and they can result in the removal of a dissolved element when the product species is involved in a phase transfer reaction. As an example of a redox reaction, consider the oxidation of Fez+, described by the following reaction: O2 + 4Fe2+ + 4H+ = 2H2O + 4Fe3+ In this reaction, the exchange of electrons changes the oxidation number of oxygen from (O) to (-II) and iron from (+IT) to (+III). For such reactions, there are reductants (electron donors, e.g., Fez+) and oxidants (electron acceptors, e.g., 02~. Fundamental problems exist in the conceptualization of redox re- actions. Most explanations make the assumption that the reactions among aqueous species are reversible and at equilibrium. However, it is abundantly clear from various lines of evidence (see, for exam- ple, Lindberg and Runnelis, 1984) that many redox reactions are essentially irreversible. That is, the reactions may go in one direc- tion easily but cannot be reversed (without biological intervention) to go in the other direction. Well-known examples include the re- action at low temperatures between such reduced-sulfur species as HS- or S2- and oxidized sulfur in the form of SO42-; the oxidation
50 GROUND WATER MODELS to SO42- readily occurs in the presence of oxygen, but the reduction of SO42- cannot take place unless sulfate-reducing bacteria are ac- tively involved. Similar examples of irreversible reactions are known for selenium, arsenic, nitrogen, and many other elements. Much re- search is needed to identify which redox reactions can be modeled as reversible reactions and which are irreversible and should not be included in an equilibrium model. Biological Transformations When organic and some inorganic compounds are present as contaminants, biological transformation can be an important pro- cess, because the original contaminant is destroyed. The advantages of biological transformation are two: (~) the contarn~nant can be completely mineralized to innocuous products, and (2) the process is not saturated, an with sorption, exchange, or filtration. Microor- ganisms can be involved with redox reactions, as described in the previous section, or with substitution and hydrolysis reactions. The microorganisms produce enzymes that allow the reactions to proceed much more rapidly. Therefore the kinetics of m~crobially mediated reactions are faster than those of the same reactions in the absence of the microorganisms. The metabolic capabilities of the microorganisms that can exist in soils and aquifers are quite diverse and can allow biodegradation of almost all types of organic contaminants. For instance, the fungi are known for their ability to degrade complex polysaccharides and other polymers of natural origin. Their capability to degrade xenobiotic (i.e., man-made) chemicals is, however, thought to be small. Bacteria have wide-ranging capability to degrade natural and xenobiotic organic compounds. Significant advances have been made over the past 10 years in elucidating the broad capabilities of bacte- ria toward xenobiotic chemicals. Table 2.2 lists classes of xenobiotic organic compounds that are known to be degraded by aerobic bac- teria, while Table 2.3 lists compounds that are degraded by strictly anaerobic bacteria. For further information, several more thorough reviews can be consulted (Alexander, 1985; Atlas, 1981; Rittmann et al., 1988~. Current research is being directed toward further defining the capabilities of bacteria for degradation of xenobiotics. Biological reactions are driven by the ultimate goal of producing new cell mass. In order to accomplish this goal, rn~croorganisms must transform environmentally available nutrients to forms that
MODELING OF PROCESSES TABLE 2.2 Classes of Xenobiotic Organic Compounds Known to Be Biodegraded by Aerobic Bacteria Unsubstituted aromatics (phenols, benzenes, benzoates) Nitro-substituted aromatics Halogen-substituted aromatics Polycyclic aromatics PCBs Most pesticides Phthalate esters TABLE 2.3 Classes of Xenobiotic Organic Compounds Known to Be Biodegraded Under Strictly Anaerobic Conditions Halogenated aliphatic solvents Most unsubstituted aromatics Some PCBs Some pesticides 51 are useful for incorporation into cells. Then they must synthesize the useful components into the polymers that make up the cell mass. The environmentally available nutrients often are not in the form needed by the cells. In general, ceils utilize reduced forms: e.g., NH4~-N, HS--S, and CH2O-C. However, the commonly available forms often are oxidized, e.g., NO3--N, SO42--S, and CO2-C. In order that the needed forms can be made, a source of electrons is necessary. Hence an essential feature for growing cells is having an electron donor. Reducing nutrients takes energy as well as electrons. Synthe- sizing the polymers that make up ceils, repairing or replacing cell constituents, transporting nutrients across the cell membrane, and motility are also significant energy sinks. Therefore cells must have an energy source if they are to grow and sustain themselves. In most cases the energy source is the electron donor. To generate energy the electron donor must donate its electrons to an electron acceptor, making available the energy for cell synthesis. Following is an example of a biologically mediated redox reaction in which an organic compound typified as CH2O is oxidized to simpler compounds: (1/4)CH2O + (1/4~02 = (1/4)CO2 ~ (1/4)H2O.
52 GROUND WATER MODELS In this reaction, the transfer of electrons between the organic com- pound CH2O (the electron donor) and O2 (the electron acceptor) provides the energy required for cell growth. In addition to oxy- gen, other potential electron acceptors could include NO3-, NO2-, SO42-, CO2, and certain organic compounds. To summarize, growth of cells requires environmentally available nutrients, an electron donor, an electron acceptor, and an energy source. If all these factors are present and if the environment is not toxic, microbial growth is possible. Whether or not a contaminant of interest fulfills any of these needs, some materials must fill the need if microorganisms are to accumulate in the environment. WHAT IS A MODEL? A mathematical mode! is a replica of some real-worId object or system. It is an attempt to take our understanding of the pro- cess (conceptual model) and translate it into mathematical terms. Therefore the mathematical mode} is only as good as our concep- tual understanding of the process. A mathematical mode! differs from other models (e.g., physical, analog) in its attempt to simulate the actual behavior of a system through the solution of mathemat- ical equations. In this sense, a mathematical mode! is much more abstract than physical models. The three main components of a mode! are the specific infor- mation describing the system of interest (e.g., what processes are important), the equations that are solved in the model, and the mode} output. A requirement for solving the equations embedded in a mode} is data about the user's particular problem. These data include specified numerical values for parameters describing the pro- cesses, for simulation parameters that are part of the procedure used to solve the equations, and for parameters describing the shape or geometry of the region of interest. This information provided by the user customizes the mode} to the particular problem. These input data in conjunction with the governing equations determine system behavior under the specified conditions. As an example of the steps in modeling, consider a ground water flow problem (Figure 2.~. Information that needs to be provided to describe a real system (Figure 2.12a) could include the following: the shape of the modeled region, the hydraulic conductivity and specific storage distributions in that region,
MODELING OF PROCESSES ~ - - ~ ~1 1 ~ K = sol ~ , ,) , , , , ~ l 53 Input: Region Shape 2. Hydraulic Conductivity 3. Boundary Conditions a Model Control Parameters ~Or | Governing | I Equation I Output: Predicted Hydraulic Head K = ~ it_ FIGURE 2.11 The components of a model: input data, a governing equation solved in the code, and the predicted distribution, which for this example is hydraulic head. the boundary and initial conditions, and ~ mode} control parameters. The mode! control parameters represent information like grid size and time step sizes that is necessary for the numerical solution of the differential equations. These details are used, given the capability of the model, to solve the partial differential equationts) describing ground water flow. Output from the mode! is a predicted hydraulic head distribution for the specified region (Figure 2.12b). Each of the different types of flow problems of interest (e.g., ground water flow, multiphase flow, and contaminant transport) will have a different governing equation, reflecting the fact that different processes are involved in ground water flow in comparison to other processes such as dissolved contaminant transport. Not unexpect- edly, then, the information that a user supplies will also be different because the various processes have unique parameters. The model control parameters will also vary from one mode! to another because,
54 E ~n - o 7:, a, 3 o (n m cl5 3 ~ O a, m ~o ~ ~ ~q C~ C _ .o ~c C E _ ,= W C C ~ ~ W ~ il., o' D .m E ~ :, _ I,o ,C Cd w ~ ~4,, _ n W cC {D ~ C W ~Q C W ~a ~ o' m ~5 3 O cn m 1 0 0 ° O 0 ~_ - CM (I S w ~ I) NOI1VA313 o 0 ~n u~ C C -~o 00 _ _ ~y> q, 1 ; ~ - iL o z ~a ~V w V . ._ _ _ ~ - p W J'-' W ~ oo^~ ~4 O R ·.,,, O ._ L' od $ ~ dR :~. O d .° P o _ ~ ~, w d R W X d · _ ~ f" ~ E E L) c, o o 11 11 c y> o · _ W w · W ·_l _,= W `_ o R- 4. O _ O X R ~ W d _ C~ ~ C~ ~ d
MODELING OF PROCESSES 55 in general, many different mathematical techniques are available to solve a given equation. The following sections examine these issues, beginning with the governing equations and how the processes are actually incorporated. Next, some of the different methods available to solve flow and transport equations are discussed. Governing Equations The development of the flow and contaminant transport equa- tions is relatively straightforward because all of the flow problems of interest here ground water flow, multicomponent flow, and dis- solved contaminant transport are developed from the same funda- mental principle, namely, the conservation of fluid or dissolved mass. Given a block of porous medium, the general conservation equation for the volume can be expressed as rate of mass input-rate of mass output + rate of mass production/consumption = rate of mass accumulation. (2.1) The differential equations for flow or mass transport are simply mathematical expressions of this conservation statement incorporat- ing the relevant processes. For fluid flow, the process responsible for moving fluid mass into and out of the volume element is simply flow in response to a potential gradient, while the rates of mass production or consumption include processes such as injection or pumping that add or remove fluid mass directly to or from the volume element. The same ideas hold for mass transport, with advection and dispersion responsible for moving mass into or out of the volume element, and the chemical and biological processes acting as source-sink terms. The details of how the governing equations are developed are relatively complicated and are best left to textbooks (Bear, 1972, 1979; Freeze and Cherry, 1979; Greenkorn, 1983~. A brief overview of how the governing equations are formulated is included here. However, those not interested in the mathematics or the details of how the equations are developed should simply be aware that equations can be developed that describe flow problems of interest; they can continue the general overview of models with the section "Solving Flow and Transport Equations." As a simple example in the development of the differential equa- tions that might be incorporated in a model, consider the simple case of ground water flow. The development begins with a mathematical statement of fluid mass conservation having the following form:
56 GROUND WATER MODELS · (~(Pq~) + (squid) + sa ] ~ pW §t ' (2.2) where p is the fluid density, q is the specific discharge or Darcy velocity,; ~ is porosity, W is a source-sink rate term, t is time, and x, y, z represents the system of Cartesian coordinates. This equation is developed by taking a small cube of porous medium and accounting for inflows and outflows, fluid mass storage, and sources or sinks (Freeze and Cherry, 19793. What is more important than the mathematical complexity of this equation is what it means. The first three terms on the left-hand side of (2.2) incorporate all the processes contributing to fluid mass or dissolved mass movement, and the fourth accounts for sources and sinks. Every problem of flow that can be modeled is described by one or more equations of this form. The continuity equation (2.3) is a general mathematical frame- work that needs to be refined by providing a more precise description of each of the processes involved. More specifically, this means ex- pressing mass fluxes in terms of the driving force (i.e., head gradients) and source-sink terms in terms of specific rate equations for the par- ticular processes. To illustrate these ideas, consider the development of the basic mass flux equations beginning with the simplest case, that of ground water flow. Freeze and Cherry (1979) show how to simplify (2.2) by inter- preting the right-hand side in terms of water released from storage because of a decline in head, and assuming fluid density is constant they divide all terms by p. With these modifications the statement of continuity can be rewritten an [. (q=) ~ (squid) + (sq8) ~ ~ W = S8§h/§t, (2.3) L ox oy on -- T AX where S., is specific storage and h is hydraulic head. The continuity equation (2.3) needs to be refined by providing a more precise description of the flow process. More specifically, this means expressing mass fluxes in terms of the driving force (i.e., head gradients). For ground water flow, the step involves replacing specific dis- charge by using the well-known Darcy equation. In its simplest form, the Darcy equation states that the flow of water through a porous medium with a unit cross-sectional area is related to the product of the hydraulic gradient and a constant of proportionality termed
MODELING OF PROCESSES 57 hydraulic conductivity. This latter parameter is related to the per- meability of the medium. Mathematically, this important physical law can be written as go =-K~§h/§x, qy =-Ky~h/bY, qx =-Kz th/§z, (2.4) where qua is the specific discharge or Darcy velocity and K is the hydraulic conductivity with components in the x, y, and z directions that are aligned with the principal material property axes. The negative signs in (2.4) mean that water is flowing in the direction opposite to increasing hydraulic potentials. Substitution of (2.4) in (2.3) provides one form of an equation suitable for modeling ground water flow, or ~ (KX ~ ~ + ~ (Ky ~ ~ + ~ (Kz ~ ~ ~ W = Seth/§t. (2.5) In this equation, the source-sink term does not need further elaboration because for the case of ground water flow, it is a simple constant related to the pumping or injection rate per unit volume. Thus (2.5) in one form or another is the differential equation used to mode! ground water flow in response to a potential gradient and subject to the effects of pumping/injection. In some cases, fluid properties such as density or viscosity vary significantly in time or space because of changes in temperature or chemical composition. When the system is nonhomogeneous, the relations among water levels, heads, pressures, and velocities are not straightforward. Calculations of flow rates and directions then require information on intrinsic permeability, density, and viscosity (rather than hydraulic conductivity) and fluid pressures and eleva- tions (rather than hydraulic heads). The same set of steps (1) development of the continuity equa- tion and (2) substitution of some form of the Darcy equation can be used with every fluid flow problem. However, for multicomponent flow problems, one continuity equation is required for each flowing fluid. Tables 2.4 and 2.5 summarize the development of basic equa- tions for unsaturated ground water flow and two-component liquid flow (organic liquid and water), respectively. The general steps are the same in both cases. The solution to the unsaturated flow equation (2.6) (see Table
58 GROUND WATER MODELS TABLE 2.4 Summary of the Steps Involved in Developing the Unsaturated Flow Equation Write the continuity equation for water an equation for the air phase is not required because the air phase is assumed to be immobile. Also no source term is assumed. ~ is the volumetric moisture content; O' is the degree of saturation. Make assumptions that let p be removed, i.e., constant density. Right-hand side has been expanded and small terms removed, with (A = em'. Substitute the Darcy equation for unsaturated flow in which hydraulic conductivity K is a function of pressure head ¢. Substitute qx, qy, qz into continuity equation. Rewrite in terms of pressure head `4, where h = ~ + z and dry =-·-or C(O at ' where C(~) is the specific moisture capacity. r qX = d(pqx) d(Pqy) . + AX BY day ~qZbo + + =- by fezat -K(O- ax qy = -K(O- dy ah qz = -K(O LIZ + 6(Pqz)1 = 5(P60') Liz I fit d dh d ah _ K - d [K(O d ] + d [K(O d ] + d [ (A d ] = ax [K(~) dx at _ ] by [ dy] + dz [K(O ( dz + 1 - C(O0¢ (2.6) NOTE: Solution of the working equation requires knowledge of the characteristic curves K(O and C(O. SOURCE: From Freeze and Cherry, 1979.
59 to s Q + ~ ~o - ~ -- - 51 -~ Cal o .~ Ct :: CQ o 3 3 ~ lo: o Em .= C C ~ O ~ O . a CUB c ,~ ° , c ~ ~ C c a 4,, <,, 8 . In : ~ ,, o y ~ O ~ E 3 + 3 it :r Q Q US Q - 1 S. Q - + UJ X ~ `7, Q 1 11 .£ S o o Q - 1 Ct ._ 4 - - ._ 1 11 ~x,1 ~ ~1 ^ 31 ~ 3 Cal UJ X Q 11 1 ._ ~ - o ~ . cd s:~ ~ c~ ~ . - ~ - ~4 C-) · cd . · ~ o ~ ~ . o . ~ ~ - >' ;^ ·~ ~ ~ ~ - o g cd . u>, s 3 ~ C ~ ~ ° ;;~. ~ 3 ~ ·~- o ,, ~ Is" ~ 1 1 C~ ~1 ~ '1 ~ 1 1 11
60 'a 11 i - oc 3 Q I ant 1 + i1 ~- 3 3 Q l ~3 =1 X it so VO in, ~ O ·- Cd ~ O O ~ Cal a' 3 o V) o 4- ~ Cd ~ o cram at, o o o Cd ;> ~ Cd ~ ,3 .b ~ U) a' V: Ct oo it ~ Q 11 1 S. Q + I1 ~ ~ 1 ~ - 1 i_ m1 ~ =1 S ~ ~ in; 1 1 ·£ ~ 1 X o
MODELING OF PROaESSES 61 2.4) provides the pressure head (`b) at selected points in the unsat- urated zone as a function of time. In Table 2.5, the two equations describing the problem of multicomponent flow, (2.7) and (2.8), can be used together with other relationships to predict pressure and saturation distributions for the two fluids as a function of time and space. Clearly, these equations do not describe all of the complexity of multiple-fluid systems. For example, Abriola and Pinder (1985a,b) describe a comprehensive approach to modeling the migration of a chemical contaminant in the nonaqueous phase, in the gas phase, and in the water phase as a soluble component. The differential equations describing the transport of mass dis- solved in ground water also are developed from conservation state- ments. However, the processes involved are quite different. Let the flux of a particular dissolved constituent into and out of a volume element of porous medium be represented by ]. The change in nota- tion to HI simply reflects that dissolved mass rather than fluid mass is being transported. The continuity equation has the following form: tbtJ) [(J) [(J) ] ~ _ [m ~ fix + TV sz J (2.9) where ~ is a source-sink term accounting for mass lost or gained within the volume element and m is the mass per unit volume. The mass of a particular contaminant dissolved in a unit volume is the concentration (C; mass per unit volume of solution) multiplied by the porosity Gil, and so (2.9) becomes ~(Ji) STAN + §(Ji)] Or FOCI [x by [x (2.10) As before, the mass flux and the source terms on the left side of (2.10) are replaced by more detailed expressions describing the processes. The mass transport of a dissolved species is controlled by three processes: advection, diffusion, and dispersion. Mathemati- cally, the advective flux of a contaminant is described by the following equation: ]$ = VxC6, (2.11) where V$ is the mean ground water velocity or seepage velocity in the x direction. Dispersive mass fluxes are commonly assumed to be driven by concentration gradients, the so-called diffusive mode} of dispersion, and as such are described by Fick's law:
62 GROUND WATER MODELS Jx=-~DX§C/6X, (2.12) where Do is the dispersion coefficient in the x direction. The d~per- sion coefficient includes both the contributions from true molecular diffusion and hydrodynamic mixing (Bear, 1972, 1979; Freeze and Cherry, 19793. A detailed discussion of the mathematical formula- tion of dispersion concepts is presented by Bear (1972, 1979) and Anderson (1979, 1984~. The combined effects of advection and dispersion can be ac- counted for simply by adding (2.~) and (2.12) to give Jx =-eD2§C/§x + VxCe, Jy =-eDy[C/by + vyCe, Jx =-eDx§C/[x + VxCe. (2.13) Substituting (2.13) into (2.10) provides a useful form of the advection-~spersion equation. Because this equation in three di- mensions is so unwieldy, the following one-dimensional version, -it; ~-EDX I; + vxCe}] ~ r = ~ I. fit (2.14) is used In the committees discussion of chemical and microbiological attenuation mechanisms. Assuming porosity to be constant in space and time and Do constant in space, (2.14) can be rewritten as 62c [C r [C Dab 2-V3 ~ +-= fit. (2.15) Equation (2.15) ~ a common form of the advection-dispersion reaction equation that is just about ready to use except for the source-sink term. It is at this point that any of the reactions of interest (e.g., radioactive decay or ion exchange) have to be specified in detail. For example, examine the case of a first-order reaction describing either radioactive decay or hydrolysis: r = d(eC)/dt = - )~eC, (2.16) where ~ is the decay constant, related to the half-life for decay. All that is required to come up with one form of a simplified contaminant transport equation ~ to substitute {2.161 into {2.151 to Rive s2c Do [$2 Vat [X ~_ (2.17)
MODELING OF PROCESSES 63 The case of mass transport accompanied by sorption that is described in terms of a simple linear isotherm can be developed in the same way. Again, an expression for the source term has to be developed by starting with the isotherm S= K&C, (2.18) where S is the quantity of mass sorbed on the surface and Ka is the distribution coefficient. The appropriate rate expression is the product of the bulk density of the medium (Paq) and the time derivative of (2.18), or [S [C -tam Fit = PaqKd fit (2.19) Substituting (2.19) into (2.15) provides the governing mass transport equation, or 82C [C Ka [C [C D2--v ~--Paq--=- Rearranging terms yields Do [s C _ in [sC = [C t1 + Palm (2.20) where the quantity in parentheses on the right side of (2.20) is a constant known as the retardation factor (Rf ). It is beyond the scope of this overview to examine how the many different chemical and biological processes are specifically developed into transport equations like (2.17~. The process is essentially the same as just described. What makes formulating equations involving complex reactions somewhat difficult is that it is necessary to write an equation like (2.17) for each contain ant that is being transported and participating in the reaction. For example, a relatively complete description of the aerobic biodegradation of an organic contaminant requires an equation describing the transport of the organic contam- inant and oxygen as well as a growth mode! describing how the mass of the microbial population changes with time (Mole et al., 1986~. Each additional equation increases the data requirements. An addi- tional problem is that some of these more complex reactions are not well understood, which adds uncertainty to the mathematical models that represent them.
64 GROUND WATER MODELS Boundary and Initial Conditions and Parameter Values The equations of flow and transport in themselves are general statements of how fluids or dissolved mass in a system should be- have as a consequence of controlling processes. However, before one proceeds to actually solve an equation, information is needed about the system. There are essentially three features that need to be described: (~) the size and shape of the region of interest, (2) the boundary and initial conditions for that region, and (3) the physical and chemical properties that describe and control the processes in the system. To illustrate these ideas, consider a problem involving the appli- cation of a steady-state ground water flow mode} to the field (Figure 2.12~. At the particular site shown in Figure 2.12, glacial till overlies sandstone and shale bedrock. The first problem to address is defini- tion of the region of interest. In other words, the lateral and vertical dimensions of the area to be modeled must be determined. The thick lines on Figure 2.12a define the area selected. The bottom boundary is assumed to coincide with the top of the Bearpaw shale, the upper boundary is the water table, and the lateral boundaries are vertical lines drawn at a major topographic high and a major topographic low. When a region has been defined, it is implicitly assumed that the rest of the geologic system can be ignored (Figure 2.12b). However, the simulation has to account for the effects of conditions outside of the region being modeled. This job is handled by the boundary conditions applied on four sides. The boundary conditions are what make it possible to isolate a specific region of interest for detailed study. There are three commonly used boundary conditions: (1) spec- ified value, (2) specified flux, and (3) value-depenclent flux (Mercer and Faust, 1981~. These are briefly described in Table 2.6. It is important to realize that every differential equation included in a mode! requires a unique set of boundary conditions. For the problem in Figure 2.12, the bottom and lateral side boundaries are assumed to be no-flow boundaries. The choice of a no- flow boundary on the bottom can be justified by geologic arguments; that is, the hydraulic conductivity of the shale is several orders of magnitude smaller than overlying units. The side boundaries are no-flow by virtue of the assumed symmetry of flow on either side (the boundaries represent flowlines). By intentionally placing these boundaries at a topographic high and a topographic low, the
MODELING OF PROCESSES TABLE 2.6 Typical Boundary Conditions for Ground Water Flow and Transport Problems r ~ . ype Specified value 65 Description Values of head, concentration, or temperature are specified along the boundary. (In mathematical terms, this is known as the Dirichlet condition.) Flow rate of water, contaminant mass, or energy is specified along the boundary and equated to the normal derivative. For example, the volumetric flow rate per unit area for water in an isotropic medium . . . Is given by ah q,l = -K an where the subscript n refers to the direction normal (perpendicular) to the boundary. (A medium that is isotropic with respect to hydraulic conductivity is equally permeable in all directions.) A no-flow (impermeable) boundary is a special case of this type in which q,, = 0. (When the derivative is specified on the boundary, it is called a Neuman condition. ) Specified flux Value-dependent flux The flow rate is related to both the normal derivative and the value. For example, the volumetric flow rate per unit area of water is related to the normal derivative of head and to head itself by -K-= q,'(hb) an where q,, is some function that describes the boundary flow rate given the head at the boundary (hb) SOURCE: Mercer and Faust, 1981. Bowlines will generally parallel these boundaries and make them no- flow boundaries. The actual decision of what the modeled region would be like in this example was made in part to provide a simple set of boundary conditions. Moving the lateral boundaries could create a much more complex set of boundary conditions. However, the modeler does have the option of placing the boundaries anywhere. Note that similar boundary conditions (e.g., specified concentrations or mass fluxes) are required to solve each mass transport equation applied to the domain. The water table boundary is an example of a specified-value
66 GROUND WATER MODELS condition. Values of hydraulic head are assumed to be known at all points along the water table. In some problems, flow across the top boundary is represented by recharge and discharge fluxes, with the configuration of the water table actually determined as part of the simulation. The last information needed about the system Is the value of parameters controlling the various flow and transport processes. As- suming the example system to be at steady state, the only parameters necessary are values of the hydraulic conductivity SKY and Ky) for each geologic unit (Figure 2.12b) and injection or withdrawal rates for sources or sinks. At this point, the flow equation is ready to be solved to provide the unknown hydraulic head at points within the region. If the flow problem in Figure 2.12 was transient, it would be necessary to provide the initial conditions, or, in other words, the distribution of hydraulic head in the region, at the start of the simulation. In addition, values for specific storage must also be specified. With this information, it would be possible to simulate the changing conditions of hydraulic head not only as a function of space but also as a function of time. The discussion of the three preparatory steps to Cudgeling is re- lated particularly to a simple problem of ground water flow. However, the same steps are followed for multicomponent flow and dissolved contaminant transport. All that change are the type and number of parameters because of the type and number of processes involved. For example, to mode} the transport of a single dissolved contami- nant that may degrade in a first-order kinetic reaction requires values of velocity Eve, vy) everywhere in the domain of interest, longitudinal and transverse dispersivities, and the decay rate constant for the re- action. In many instances, it will be necessary to run a flow mode! to provide the necessary description of the velocity field. As indicated previously, several transport equations may be necessary, depending on the complexity of the degradation reaction, e.g., its dependence upon hydrogen, oxygen, or other substrates. Sol~rmg Flow and Transport Equations There are two basic ways to solve the flow and transport equa- tions. The analytical methods embody classical mathematical ap- proaches that have been used for more than 100 years to deal with differential equations. The numerical approaches have also existed
MODELING OF PROCESSES 67 for many years but were not fully exploited until the development of computers to solve approximate forms of the governing equations. The greatest strength of the analytical methods lies in their capabil- ity in many cases to produce exact solutions to a flow or transport problem in terms of the controlling parameters. Being able to estate fish the functional form of the solution, and the interrelationships among parameters, provides a great deal of physical insight into how the processes control flow and transport. Another useful way in which the analytical solutions are used Is to provide a check on the accuracy of numerical models, which can be subject to a variety of different errors. In terms of their usefulness in solving practical problems, the numerical approaches are superior to the analytical methods because the user can let the controlling parameters vary in space and time. This feature enables detailed replications of the complex geologic and hydrologic conditions that exist in nature. Analytical methods have a role to play in field applications (e.g., theory of well hydraulics), but, in general, they are appropriate only for a narrow range of simple problems. Practical problems involving the flow of more than one fluid or contaminant are sufficiently complex that only numerical approaches are suitable. Nearly all the numerical procedures involve replacing the contin- uous form of the governing differential equation by a finite number of algebraic equations. To develop these equations, it is necessary to subdivide the region into pieces. For the flow example discussed previously (Figure 2.12), the region can be subdivided by using rect- angles (Figure 2.13~. Other geometric shapes (e.g., triangles and quadrilaterals) are also used, depending on the solution technique. For transient problems, it is also necessary to subdivide the total simulation period into a number of smaller time steps. For the example problem, values of hyciraulic head are calculated at the nodes, located at the center of ceils, with one algebraic equation written for each node. Hydraulic conductivity values are supplied for each rectangular cell. This flexibility in assigning parameter values helps create in the mode! a Attribution of geohydrologic properties that closely approximates that observed in the field. A variety of analytical and numerical solutions have been de- veloped for use in ground water applications, and a comprehensive discussion of each would require a modeling textbook. However, a summary of techniques that are commonly applied to the solution of various flow and mass transport problems is provided in Table
68 G IIJ m z ~ 10 o a: 1d GROUND WATER MODELS , Constant Head Nodes it. . . ~Kh= ~=00 Kh = v = 10-5 cm/s Kh = 10-7 cm/s; Kv = 10~ cm/s Kh = It = 10 3 cm/s COLUMN NUMBERS FIGURE 2.13 Model system from the previous figure subdivided by a rect- angular grid system. Nodes are defined in the center of each grid cell. By assigning Kh = Kv = 0.0 in the area at the top left, it is effectively excluded from the calculation. 2.7. Included are a brief description of the methods and a few key references that can be used to obtain more detailed information. The last topic that needs to be addressed in this section involves the mathematical techniques for solving matrix equations. In most models, the major computational effort comes in solving the system of mode! equations. In general, there are two basic methods. In one approach the entire system of equations is solved simultaneously with direct methods, providing a solution that is exact, except for machine round-off error. In the second approach, iterative methods obtain a solution by a process of successive approximation, which involves making an initial guess at the matrix solution and then improving this guess by some iterative process until an error criterion is satisfied. Direct methods have two main disadvantages. The first is that a computer may not be able to store the large matrices or solve the system in a reasonable time when the number of nodes is large. Sometimes this problem can be dealt with to some extent by us- ing sparse matrix solvers and various node-numbering schemes. The second problem with direct methods is the round-off errors. Be- cause many arithmetic operations are performed, round-off errors can accumulate and significantly influence results for certain types of matrices. Iterative schemes avoid the need for storing large matrices. This
69 o EM cat cat cd a' c°~ cat ._ cd 3 o .= o cat - o o cat . _ au o o v, o o an o m 2 co Do car ~I Go ~^ ~_ '_ ~ C)d _ ~ c ~ ~ ~x 2 ~ ~ = ~ ~ ~_ ~ jet ~· _ o |, ~ e ~ · 3 ~ 3 ~ 3- ~ v: cat c or .o ~,3 ~ s ~ c u: Em ¢ 3 o sit a' cd 3 0
70 s is in He o au o 4- o . . ~ a) ~ .= ·4- ~ o=o - o ~ - ~ ^ ^ - ~ oo E Go, ~ ~ ox X on X w ~ =,^ X ~ ~ . id ~ 0N 0) ·0 ~ ~ ~ ~ ~ O ~ ~ X ~0 E o ~ E ~ ~ E E ~ ~ ~ ~ ~ it_ ~ ~ ~ it. ~ ~ 3 o ~ v o ~ · ~ ILL C:: (O ~ ~ E~ ~ C ~, E ~ ~ ~ " 3 ~ o con ~ ~ ~<u ~ ~ ~a' ~ ~Z v, ~cn ~v' ~ ~C) O ~al ~ ~ ~ad ;^ ~ <( C ~ ~<t o 3 3 he o on ~ o 3 ~ == ·= ~3 o ~ ~ ~.C~ ~
71 Go - ~ oo To 'x ~ - ~ ~ -= .= - ~ ~ c == ^ 3 E hi, ~ ~ 0 ~ ~ O ~ ~ ~ ~ m ~ ~ of; ._ ._ ct v, so . ~ a' > Q Ct Go Do cry - . ^ ~ 00 3 car $, _ ~ Is: 0 4- ~ so ~ ct i, cat c~ it,) ~0 In ~ ~ ~ ~ , c, ~ O &.5 9 Y . ~ ~ y _ m ~0 ~ se c,0 ~ · s,.( cat ~0 ~ 3 ,~ if: (=,, ~ c ~ .= ~ ~ cat ~ ·= S ~ S O O - C) ~ O > O Cal ~ _
72 GROUND WATER MODELS feature makes them attractive for solving problems with many un- knowns. Numerous schemes have been developed; a few of the more commonly used ones include successive overrelaxation methods (Varga, 1962), the alternating-direction implicit procedure (Douglas and Rachford, 1956), the iterative alternating-direction implicit pro- cedure (Wachpress and Habetler, 1960), and the strongly implicit procedure (Stone, 1968~. Because operations are performed many times, iterative methods also suffer from potential round-off errors. The efficiency of iterative methods depends on an initial estimate of the solution. This makes the iterative approach less desirable for solving steady-state problems (Narasimhan et al., 1978~. To speed up the iterative process, relaxation and acceleration factors are used. Unfortunately, the definition of best values for these factors commonly is problem dependent. In addition, iterative approaches require that an error tolerance or convergence criterion be specified to stop the iterative process. This, too, may be problem dependent. All of these parameters must be specified by the mode} user. According to Narasimhan et al. (1977) and Neuman and Nara- simhan (1977), perhaps the greatest limitation of the iterative schemes is the requirement that the matrix be well conditioned. An ill-conditioned matrix can drastically affect the rate of conver- gence or even prevent convergence. An example of an ill-conditioned matrix is one in which the main diagonal terms are much smaller than other terms in the matrix. More recently, a semi-iterative method has gained popularity (Gresho, 1986~. This method, or class of methods known as conju- gate gradient methods, was first described by Hestenes and Stiefel (1952~. It is widely used to solve linear algebraic equations where the coefficient matrix is sparse and square (Concus et al., 1976~. One advantage of the conjugate gradient method is that it does not require the use or specification of iteration parameters, thereby elim- inating this partly subjective procedure (Manteuffe} et al., 1983~. Kuiper (1987) compared the efficiency of 17 different iterative meth- ods for the solution of the nonlinear three-dimensional ground water flow equation. He concluded that, in general, the conjugate gradient methods did the best. Numerical methods by their very nature, yield approximate so- lutions to the governing partial differential equations. The accuracy of the solution can be significantly affected by the choice of numerical parameters, such as the size of the spatial discretization grid and the length of time steps. Those using ground water models and those
MODELING OF PROCESSES 73 making management decisions based on mode! results should always be aware that trade-offs between accuracy and cost will always have to be made. If the grid size or time steps are too coarse for a given problem, it is possible to generate a numerical solution that converges on an answer that has an excellent mass balance but is still inaccu- rate. furthermore, if iteration parameters are not properly specified, the solution may not converge. It is hoped that this will show up as a mass balance error, which will be noted by the user. This indicates, however, the importance of a mass balance in numerical models. REFERENCES Abriola, L. M., and G. F. Pinder. 1985a. A multiphase approach to the mod- eling of porous media contamination by organic compounds, 1. Equation development. Water Resources Research 21~1), 11-18. Abriola, L. M., and G. F. Pinder. 1985b. A multiphase approach to the modeling of porous media contamination by organic compounds, 2. Numerical simulation. Water Resources Research 21~1), 19-26. Ahlstrom, S. W., H. P. Foote, R. C. Arnett, C. R. Cole, and R. J. Serne. 1977. Multicomponent Mass Transport Model: Theory and Numerical Im- plementat ion (discrete-parcel-random-walk version). BNWL-2 127. Battelle Northwest Laboratories, Richland, Wash. Alexander, M. 1985. Biodegradation of organic chemicals. Environmental Science and Technology 19, 106. Anderson, M. P. 1979. Using models to simulate the movement of contaminants through ground water flow systems. Critical Reviews in Environmental Control 9~2), 97-156. Anderson, M. P. 1984. Movement of contaminants in groundwater: Groundwa- ter transport-advection and dispersion. Pp. 37-45 in Groundwater Con- tamination. Studies in Geophysics. National Academy Press, Washington, D.C. Atlas, R. M. 1981. Microbial degradation of petroleum hydrocarbons: An environmental perspective. Microbiological Reviews 45, 180. Baehr, A. L., and M. Y. Corapcioglu. 1987. A compositional multiphase model for groundwater contamination by petroleum products, 2. Numerical solution. Water Resources Research 23~1), 201-214. Bates, J. K., and W. B. Seefeldt. 1987. Scientific Basis for Nuclear Waste Management X. Materials Research Society, Symposium Proceedings, Vol. 84, 829 pp. Bear, J. 1972. Dynamics of Fluids in Porous Media. Elsevier, New York, 764 PPe Bear, J. 1979. Hydraulics of Groundwater. McGraw-Hill, New York, 569 pp. Bredehoeft, J. D., and G. F. Pinder. 1973. Mass transport in Bowing ground- water. Water Resources Research 9, 194-210. Cherry, J. A., R. W. Gillham, and J. F. Pickens. 1975. Contaminant hydroge- ology: Part 1, Physical processes. Geoscience Canada 2~2), 76-84. Cleary, R. W., and M. J. Ungs. 1978. Groundwater Pollution and Hydrology, Mathematical Models and Computer Programs. Rep. 78-WR-15, Water Resources Program, Princeton University, Princeton, N.J.
74 GROUND WATER MODELS Concus, P., G. Golub, and D. O'Leary. 1976. Sparse Matrix Computations. Academic Press, New York, pp. 300332. Corapcioglu, M. Y., and A. L. Baehr. 1987. A compositional multiphase model for groundwater contamination by petroleum products, 1. Theoretical considerations. Water Resources Research 23~1), 191-200. DeMarsily, G. 1986. Quantitative Hydrogeology. Academic Press, Orlando, Fla., 440 pp. Domenico, P. A., and G. A. Robbins. 1985. A new method of contaminant plume analysis. Ground Water 23~4), 476-485. Douglas, J., Jr., and H. H. Rachford, Jr. 1956. On the numerical solution of heat conduction problems in two and three space variables. Transactions of the American Mathematics Society 82, 421-439. Faust, C. R. 1985. Transport of immiscible druids within and below the un- saturated zone: A numerical model. Water Resources Research 21~4), 587-596. Freeze, R. A. 1969a. Regional ground water Bow Old Wives Lake drainage basin, Saskatchewan. Inland Waters Branch, Department of Energy, Mines, and Resources, Canada, Scientific Series, No. 5, 245 pp. Freeze, R. A. 1969b. Theoretical analysis of regional ground water Bow. Inland Waters Branch, Department of Energy, Mines, and Resources, Canada, Scientific Series, No. 3. Freeze, R. A. 1971a. Influence of the unsaturated flow domain on seepage through earth dams. Water Resources Research 7~4), 929-941. Freeze, R. A. 1971b. Three-dimensional, transient, saturated-unsaturated Bow in a groundwater basin. Water Resources~Research 7~2), 347-366. Freeze, R. A., and J. A. Cherry. 1979. Ground Water. Prentice-Hall, Englewood Cliffs, N.J., 604 pp. E`reyberg, D. L. 1986. A natural gradient experiment on solute transport in a sand aquifer: II. Spatial moments and the advection and dispersion of non-reactive tracers. Water Resources Research 22~13), 2031-2046. Fried, J. J. 1975. Ground Water Pollution: Theory, Methodology, Modeling, and Practical Rules. Elsevier, Amsterdam, 330 pp. Frind, E. O. 1987. Simulation of ground water contamination in three dimen- sions. Pp. 749-763 in Proceedings of Solving Ground Water Problems with Models. National Water Well Association, Denver, Colo. Galloway, W. E., and D. K. Hobday. 1983. Terrigenous Clastic Depositional Systems. Springer-Verlag, New York, 423 pp. Greenkorn, R. A. 1983. Flow Phenomena in Porous Media: Fundamentals and Applications in Petroleum, Water and Food Production. Marcel Dekker, New York, 550 pp. Gresho, P. M. 1986. Time integration and conjugate gradient methodsfor the incompressible Navier-Stokes equations. Pp. 3-27 in Finite Elements in Water Resources, Proceedings of the 6th International Conference, Lisbon, Portugal. Springer, Berlin. Hanks, R. J., A. Klute, and E. Bresler. 1969. A numeric method for estimating infiltration, redistribution, drainage, and evaporation of water from soil. Water Resources Research 5, 1064-1069. Hantush, M. S. 1964. Hydraulics of wells. Advances in Hydroscience 1, 281-432. Hestenes, M., and E. Stiefel. 1952. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards 49~6), 409-436.
MODELING OF PROOESSES 75 Huyakorn, P. S., and G. F. Pinder. 1983. Computational methods in subsurface flow. Academic Press, New York, 473 pp. Jackson, R. E., R. J. Patterson, B. W. Graham, J. Bahr, D. Belanger, J. Lockwood, and M. Priddle. 1985. Contaminant Hydrogeology of Toxic Organic Chemicals at a Disposal Site, Gloucester, Ontario. 1. Chemical Concepts and Site Assessment. NHRI Paper No. 23, Iwo Scientific Series No. 141, National Hydrology Research Institute, Inland Waters Directorate, Environment Canada, Ottawa, Canada, 114 pp. Jacob, C. E. 1940. On the flow of water in an elastic artesian aquifer. Traneac- tions, American Geophysical Union 2, 574-586. Javandel, I., C. Doughty, and C.-F. Twang. 1984. Groundwater Transport: Handbook of Mathematical Models. Water Resources Monograph 10, Amer- ican Geophysical Union, Washington, D.C., 228 pp. Jeppson, R. W. 1974. A'cisymmetric infiltration in soils Numerical techniques of solution. Journal of Hydrology 23, 111-130. Konikow, L. F., and J. D. Bredehoeft. 1978. Computer model of two-dimensional solute transport and dispersion in ground water. Technical Water Re- sources Inventory, Book 7, Chap. C2. U.S. Geological Survey, Reston, Va., 90 pp. Kruseman, G. P., and N. A. de Ridder. 1983. Analysis and Evaluation of Pumping Test Data. Bulletin 11. International Institute for Land Recla- mation and Improvement/ILRI, P.O. Box 45, G700 AA Wageningen, The Netherlands. Kuiper, L. K. 1987. A comparison of iterative methods as applied to the solution of the nonlinear three-dimensional groundwater flow equation. SIAM Journal of Scientific and Statistical Computing 8~4), 521-528. Lappala, E. G. 1980. Modeling of water and solute transport under variably saturated conditions: State of the art. Prepared for Proceedings of the Interagency Workshop on Radioactive Waste Modeling, December 2-4, 1980, Denver, Cola. Lenhard, R. J., and J. C. Parker. 1987. A model for hysteretic constitutive relations governing multiphase ~ow, 2. Permeability-saturation relations. Water Resources Research 23~12), 2197-2206. Liggett, J. A., and P. ~. Liu. 1983. The Boundary Integral Equation Method for Porous Media Flow. George Allen and Unevin, London. Lindberg, R. D., and D. D. Runnelle. 1984. Ground water redox reactions: An analysis of equilibrium state applied to Eh measurements and geochemical modeling. Science 225, 925-927. Lohman, S. W. 1979. Ground-water hydraulics. U.S. Geological Survey Profes- sional Paper 708. U.S. Government Printing Office, Washington, D.C. Mackay, D. M., D. L. E`reyberg, P. V. Roberts, and J. A. Cherry. 1986. A natural gradient experiment on solute transport in a sand aquifer, 1. Approach and overview of plume movement. Water Resources Research 22~13), 2017, 2029. ManteuRel, T. A., D. B. Grove, and L. F. Konikow. 1983. Application of the conjugate-gradient method to ground-water models. R~p. 83-4009, Water Re~ources Investigations, U.S. Geological Survey, 24 pp. Mercer, J. W., and C. R. Faust. 1981. Ground-Water Modeling. National Water Well Association, Worthington, Ohio, 60 pp.
76 GROUND WATER MODELS Molz, F. J., M. A. Widdowson, and L. D. Benefield. 1986. Simulation of microbial growth dynamics coupled to nutrient and oxygen transport in porous media. Water Resources Research 22, 1207-1216. Morel, F. M. M. 1983. Principles of Aquatic Chemistry. Wiley, New York, 446 PP Narasimhan, T. N., S. P. Neuman, and A. L. Edwards. 1977. Mixed explicit- implicit iterative finite element scheme for diffusion-type problems: II. Solution in strategy and examples. International Journal for Numerical Methods in Engineering 11, 325-344. Narasimhan, T. N., P. A. Witherspoon, and A. L. Edwards. 1978. Numerical model for saturated-unsaturated flow in deformable porous media, 2. The algorithm. Water Resources Research 14~2), 255-261. Neuman, S. P. 1972. Finite element computer programs for flow in saturated- unsaturated porous media. Second Annual Report. Project No. A10-SWC- 77, Hydraulic Engineering Laboratory, Technion, Haifa, Israel, p. 87. Neuman, S. P., and T. N. Narasimhan. 1977. Mixed explicit-implicit iterative finite element scheme for diffusion-type problems: I. Theory. International Journal for Numerical Methods in Engineering 11, 309-323. Nielsen, D. R., M. Th. Van Genuchten, and J. W. Biggar. 1986. Water flow and solute transport processes in the unsaturated zone. Water Resources Research 22~9), 89S-108S. Ogata, A. 1970. Theory of dispersion in a granular medium. U.S. Geological Survey Progress Paper 411-I. Osborne, M., and J. Sykes. 1986. Numerical modeling of immiscible organic transport at the Hyde Park landfill. Water Resources Research 22~1), 25-33. Parker, J. C., and R. J. Lenhard. 1987. A model for hysteretic constitutive relations governing multiphase flow, 1. Saturation-pressure relations. Water Resources Research 23~12), 2187-2196. Philip J. R. 1955. Numerical solution of equations of the diffusion type with diffusivity concentration dependent. Transactions of the Faraday Society 51,885-892. Phillip, J. R. 1957. The theory of infiltration: 1. The infiltration equation and its solution. Soil Science 83,345-357. Pinder, G. F., and W. G. Gray. 1977. Finite Element Simulation in Surface and Subsurface Hydrology. Academic Press, New York, 295 pp. Prickett, T. A., T. G. Naymik, and C. G. Lounquist. 1981. A Random-Walk Solute Transport Model for Selected Groundwater Quality Evaluations. Bulletin 65, Illinois State Water Survey, Champaign, Ill., 103 pp. Pruess, K., and R. C. Schroeder. 1980. SHAFT79, User's Manual. LBL-10861, Lawrence Berkeley Laboratory, University of California, Berkeley. Reddell, D. L., and D. K. Sunada. 1970. Numerical simulation of dispersion in groundwater aquifer. Hydrol. Paper 41, Colorado State University, Fort Collins, 79 pp. Reeves, M., D. S. Ward, P. A. Davis, and E. J. Bonano. 1986a. SWIFT II Self-Teaching Curriculum: Illustrative Problems for the Sandia Waste- Isolation Flow and Transport Model for Fractured Media. NUREG/CR- 3925, SAND84-1586, Sandia National Laboratory, Albuquerque, N. Mex. Reeves, M., D. S. Ward, N. D. Johns, and R. M. Cranwell. 1986b. Data In- put Guide for SWIFT II, The Sandia Waste-Isolation Flow and Transport
MODELING OF PROCESSES 77 Model for Fractured Media. NUREG/CR-3162, SAND83-0242, Sandia National Laboratories, Albuquerque, N. Mex. Reeves, M., D. S. Ward, N. D. Johns, and R. M. Cranwell. 1986c. Theory and Implementation for SWIFT II, The Sandia Waste-Isolation Flow and Transport Model for Fractured Media. NUREG/CR-3328, SAND83-1159, Sandia National Laboratories, Albuquerque, N. Mex. Reisenauer, A. E. 1963. Methods for solving problems of partially saturated steady flow in soils. Journal of Geophysical Research 68, 5725-5733. Remson, I., G. M. Hornberger, and F. J. Molz. 1971. Numerical Methods in Subsurface Hydrology. Wiley, New York, 389 pp. Rittmann, B. E., D. Jackson, and S. L. Storck. 1988. Potential for treatment of hazardous organic chemicals with biological processes. Pp. 15-94 in Biotreatment Systems, Vol. III, D. L. Wise, ed. CRC Press, Boca Raton, Fla. Runnells, D. D. 1976. Wastewaters in the Vadose zone of arid regions: Geo- chemical interactions. Ground Water 14~6), 374-385. Schwartz, F. W. 1975. On radioactive waste management: An analysis of the parameters controlling subsurface contaminant transport. Journal of Hydrology 27, 51-71. Schwartz, F. W. 1977. Macroscopic dispersion in porous media: The controlling factors. Water Resources Research 13~4), 743-752. Schwartz, F. W. 1984. Modeling of ground water {low and composition. Pp. 178-188 in Proceedings of First Canadian/American Conference in Hydro- geology, B. Hitchon and E. I. Wallich, eds. National Water Well Association, Banff, Alberta. Schwartz, F. W., and A. S. Crowe. 1980. A deterministic-probabilistic model for contaminant transport. NUREG/CR-1609, Nuclear Regulatory Com- mission, 158 pp. Schwille, F. 1984. Migration of organic fluids immiscible with water in the unsaturated zone. Pp. 27-48 in Pollutants in Porous Media: The Unsatu- rated Zone Between Soil Science and Groundwater, B. Yaron, G. Dagan, and J. Goldshmid, eds. Ecological Studies, Vol. 47. Springer-Verlag, Berlin. Smith, L., and F. W. Schwartz. 1980. Mass transport, 1. A stochastic analysis of macroscopic dispersion. Water Resources Research 16~2), 303-313. Stone, H. K. 1968. Iterative solution of implicit approximations of multidi- mensional partial differential equations. Society of Industrial and Applied Mathematics, Journal of Numerical Analysis, 5~3), 530-558. Sudicky, E. A. 1986. A natural gradient experiment on solute transport in a sand aquifer: Spatial variability of hydraulic conductivity and its role in the dispersion process. Water Resources Research 22~13), 2069-2082. Tennessee Valley Authority. 1985. A review of field-scale physical solute trans- port processes in saturated and unsaturated porous media. EPRI EA-4190, Project 2485-5, Electric Power Research Institute, Palo Alto, Calif. Theis, C. V. 1935. The relationship between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage. Transactions, American Geophysical Union 2,519-524. van Genuchten, M. T., and W. J. Alves. 1982. Analytical solutions of the one-dimensional convective-dispersive solute transport equation. Technical Bulletin 1661, U.S. Department of Agriculture, Washington, D.C., 149 pp. Varga, R. S. 1962. Matrix Iterative Analysis. Prentice-Hall, Englewood Cliffs, N.J., 322 pp.
78 GROUND WATER MODELS Voss, C. I. 1984. SUTRA-Saturated Unsaturated Transport A finite-element simulation model for saturated-unsaturated Iluid-density-dependent ground- water Bow with energy transport or chemically-reactive single-species solute transport. Water Resources Investigations Rep. 84-4369, U.S. Geological Survey, Reston, Va., 409 pp. Wachpress, E. L., and G. J. Habetler. 1960. An alternating-direction-implicit iteration technique. Journal of Society of Industrial and Applied Mathe- matics 8, 403-424. Walton, W. C. 1970. Groundwater Resource Evaluation. McGraw-Hill, New York, 664 pp. Wang, J. F., and M. P. Anderson. 1982. Introduction to Groundwater Modeling. Freeman, San Francisco, Calif., 237 pp. Welch, J. E., F. H. Harlow, J. P. Shannon, and B. J. Daly. 1966. The MAC Method, A Computing Technique for Solving Viscous, Incompressible, Transient Fluid-Flow Problems Involving Free Surfaces. LA-3425, Los Alamos Scientific Laboratory of the University of California, Los Alamos, N. Mex. Yeh, G. T., and D. S. Ward. 1980. FEMWATER: A Finite-Element Model of Water Flow Through Saturated-Unsaturated Porous Media. ORNL-5567, Oak Ridge National Laboratory, Oak Ridge, Tenn. Yeh, G. T., and D. S. Ward. 1981. FEMWASTE: A Finite-Element Model of Waste Transport Through Saturated-Unsaturated Porous Media. ORNL- 5601, Oak Ridge National Laboratory, Oak Ridge, Tenn. Zienkiewicz, O. C. 1977. The Finite Element Method, 3rd ed. McGraw-Hill, London.
