A mathematical model for the transport and fate of organic chemicals in unsaturated/saturated soils.

A mathematical model, simulating the transport and fate of nonionizable organic compounds in unsaturated/saturated porous media (soils) in a terrestrial microcosm has been developed. Using the principles of water mass, momentum, heat energy and chemical mass balance, the three fields: moisture, temperature, and liquid phase chemical concentration are solved for simultaneously by coupling the soil slab to an environmentally realistic air-soil interface (a dynamic free boundary) conditions and a prescribed height water table. The environmental conditions at the soil surface-air chamber interface are easily changed, via geometric scaling factors, to simulate either an open agricultural field or a landfill type of situation. Illustrative simulation runs examine the effects of different soil-chemical characteristics on hydrological and chemical concentration profiles.


Introduction
Transport and fate models for organic chemicals in soils have been in the literature for years. They are almost always, however, specific models, which have been designed to simulate the transport and fate of organic compounds in isothermal, homogeneous, water saturated, sorbing porous media. Piver and Lindstrom (1) and Lindstrom and Piver (2) have surveyed the existing literature on these models. To date, there does not appear to be any published transport and fate model for the nonisothermal, nonhomogeneous, unsaturated/ saturated sorbing porous media case except that of Lindstrom and Piver (2).
It is no longer adequate to ignore the influences of moisture and temperature variations (spatially and temporally) on chemical reactions and/or biological processes in soils. It is well known that microbial degradation rates are very temperature sensitive for most organic compounds. It therefore is incumbent upon the transport and fate modeler to retain as much realism, with respect to known physical and biological processes, as is possible. Just such a model, called NEWTMC, has been evolving for several years now. Presenting  features of this model is the subject of this paper. Only the major assumptions and resulting transport equations are presented below. The reader, interested in the mathematical details of the integrated system of transport and fate equations, should consult Lindstrom and Piver (2). Also, the current listing of the Fortran V computer code together with an example run both on a CDC Cyber 172/720 and a DEC VAX can be found in Lindstrom and Piver (2).

System Assumptions and Major Features
NEWTMC is the name of the Fortran V coded algorithm for simulating the transport and fate of nonionizable organic compounds in a terrestrial microcosm (TMC). Although the model was originally designed to simulate the transport and fate of a chemical in a terrestrial microcosm whose scale is about 1 to 2 m3 in volume (Fig. 1), the version described here is general enough to span the range of physical situations from small TMCs to large fields or landfills. The large field simulation is obtained by geometric parameter specification. Figure 2 shows a diagrammatic sketch of the CERL TMC, while Figure 3 shows a thin slab of soil in the TMC within which certain physical and biological transport and chemical conversion processes are-occurring. A nomenclature list explaining all the mathematical symbols can be found in Tables 1 and 2. Basic Assumptions for NEWTMC For the air chamber, the following assumptions have been made: (1) the contents of the air chamber are instantaneously well mixed; (2) the air flow in the box is turbulent and not laminar; (3) the air inflow rate is constant in time; (4) there are no sources or sinks of air in the system anywhere; (5) the incoming air temperature and relative humidity can be characterized by Fourier series methods; (6) the light source can be readily characterized over time via off-on switch type step functions (indicator functions); (7) the chemical input and/or initial distribution is known; (8) the rain intensity and frequency are prescribed; (9) the air chamber is in intimate contact with the soil surface; For the soil chamber the following assumptions have been made. (1) The soil can be characterized in terms of known parameters (parameteric functions of space and time possibly) such as the local bulk density Pb, the local porosity E, the local percent, sand, clay and organic matter present, the local specific heats of the solids and water, the local heat conductivities of the soil solids and water, the local latent heat of vaporization, the local soil moisture tension function, the local soil hydraulic con-ductivity function and the density of water (assumed space-time and temperature independent) water vapor, and saturated water vapor in the soil voids. (2) The air chamber and soil surface are coupled together for all three fields-moisture, heat and chemical-via heat and mass transport rules operating in the top boundary layer. (3) The water table is spread uniformly over the bottom of the TMC and is assumed to be at external atmospheric conditions for water, heat, and dissolved chemical, and chemical vapor transmission at all times. (4) The initial chemical distribution in the soil is prescribed. (5) There is no flow of moisture and no flow of chemical through the lateral walls of the TMC. (6) The macroscopic transport equations are the outcomes of averaging the microscopic field equations (stochastic) over volumes of at least 1 cm3. heat and chemical fields in both the air chamber over the soil slab and the soil slab itself. (2) The model simulates an open field situation or those conditions over a landfill by setting the height of the air chamber to a very large value as well as setting qa, the volumetric air flow rate through the air chamber, simultaneously large so WS, the wind speed, remains fixed at the desired value. (3) The timing and intensity of both light and rain events can be very freely specified a priori. (4) A great many different atmospheric conditions (highlow humidity, high-low temperature, high-low wind speed) can be represented with ease. (5) Chemical, water and heat are all balanced and accounted for at all points in time and space with this balance information printed out at user specified times. (6) The free air-soil surface boundary is modeled dynamically so as to allow field variables to seek their appropriate values (they are not rigidly fixed or otherwise specified over time).
(7) Chemical may be administered via: incoming air (vapor phase), rain water, flow in from the water table or initially distributed in the soil slab. (8) The soil slab need not be spatially homogeneous but may in fact grade from one soil type to another with the latitude of soil properties being input information.

System Equations
In the air chamber the moisture balance equation is given by Eq. (1): dXmajdt = Qmwaain -Qmwaaot + Qevapi + Qvap (1) where dXma/dt is the instantaneous rate of change of water mass in air chamber per hour (g/hr).
Qmwaain, the mass of water input via incoming air per hour (g/hr) is given as: qPsatrI(Tai in(t)]RH ( (2) Qmwaain = WaPVw 'air in( air in(t) Qmwaaot, the mass of water lost via air per hour (g/hr) is given by Qevap, the evaporation of water from soil system (g/hr) Qevap, =AEconOo(l + 13vapWS) [pvvv(To)h -P" V(Ta)Oa] (4) Qvap, the diffusion of water vapor from the soil surface is given by and (6) Qevap = Qvapi + Qevap with Econ = 0.0Olp0,/pair (7) The relative humidity expression h=pWV(T)/pWV(T) = exp{tg/RT} (8) comes from Edlefson and Anderson (3). We represent the saturated water vapor density via a low-order spline function, with the parameters obtained, via least-square regression, over the range 0°C-50°C (4). The heat balance for the air chamber is given by Eq.
(9): dXhaldt =Qheatin -Qheatot + Qhtevp + Qhtlwrs -Qhtlwra + Qhtssl + Qhtswr -Qhtwalls (9) where dXhaldt is the instantaneous rate of change of heat in air chamber per hour (cal/hr), Qheatin is the heat input via incoming moist air (cal/hr), Qheatot is the heat output via exiting moist air (cal/hr), Qhtevp is the heat input via evaporation of water from the transfer of latent heat from soil surfaces (cal/hr), Qhtlwrs is the heat input via long wave radiation from the top of the soil slab (cal/hr), Qhtlwra iS the heat output via long wave radiation to the top of the soil slab (cal/hr), Qht58l iS the heat input via transfer of sensible heat to air from top of soil (a signed quantity) (cal/hr), Qhtswr iS the heat input via reflection of shortwave radiation off the top of soil slab (cal/hr) and Qhtwalls is the heat output or input via conduction through the air chamber walls (a signed quantity) (caL/hr). The various heat flow functions in Eq. (9) are defined analogously to those of the moisture field introduced in Eq. (1). The complete set of definitions of all heat flow functions can be found in Lindstrom and Piver (2). The chemical balance equation in the air chamber is given by: dXcaldt =Qmchain -Qmchaot + Qmcvls + Qmcv2s + Plach -Plap (10) where dXcaldt is the instantaneous rate of change of chemical mass (vapor phase) (,ug) per unit time (hour), Qmchain is the chemical vapor mass input via incoming air (,ug/hr), Qmchaot iS the chemical vapor output via exiting air (,ug/hr), Qm5vls is the chemical input via volatilization from liquid phase from the top of the soil slab (,ig/hr), Qmcv2s iS the chemical input via vapor-phase chemical transports out of the top of the soil slab (,ug/ hr) (signed quantity), Plach is the chemical lost via hydrolysis processes in air chamber (,ug/hr) and Plap is the chemical lost via photolysis processes in the air chamber (tig). Again, the chemical vapor mass flows, have been completely defined in Lindstrom and Piver (2).

Air Chamber Relative Fields
The air chamber moisture content (relative humidity) is defined via the equation: (11) while the temperature field (°K) in the air chamber is defined by the expression: and the chemical concentration in the air chamber (g/cm3 air) is defined via the rule:

Soil Slab Field Equations
In NEWTMC we assume there are no sources or sinks of water within the soil chamber and that elemental soil volumes and surfaces are neither deforming nor translating in space or time. All sources and sinks of moisture occur in the boundary conditions and the arbitrary elemental volume v is time-independent (Fig. 2).
Using the balances of water mass and momentum in v we are ultimately lead to the conservation equation:  (19) with the individual components defined as: Lindstrom and Piver (2) discuss the representations for the latent heat function, £(T), and the average soil solids heat conductivity function, Xofid,I(z), as well as the functions, Csolidr(Z), and PSOIidS(Z).

Chemical Field Equations
The total chemical flux (g chemical/cm2-hr) through an incremental portion of the surface E can be written as: (16) and: where the flow path averaged liquid phase chemical flux is: q=-QDC VC1 + V1C1 (24) and the flow path averaged vapor phase chemical flux is represented by the expression: The six indicated second-rank tensor field functions have been defined in Lindstrom and Piver (2).

Heat Field Equations
By defining the heat balance equation (integral form) in analogy with the moisture balance equation we arrive at the form: Assuming linear sorption and/or rapid equilibrium partitioning of the organic compound into the organic portions of the soil (e.g., dissolution into humus), a simple Henry's law partition between the chemical in the liquid phase and the vapor phase, and three types of first-order loss processes such as biodegradation, hydrolysis and irreversible sorption and/or first-order with q = q, + q, (15) (23) (17) (25) We remind the reader that the coefficients PI, Yl, al and av are all usually strongly temperature-dependent. It is very easy to incorporate this temperature dependence in NEWTMC for the three process laws for loss of chemical (structural decomposition or irreversible sorption) from each infinitesimal volume element because they are defined in subroutine FLOWS and Fortran code changes are easily made.

Boundary and Initial Conditions
Before we can solve the system of six coupled nonlinear equations [(1), (9), (10), (14), (18) and (26) (2) give a detailed discussion of the boundary conditions currently built into NEWTMC. It suffices to say that the "free boundary" which exists at the air chamber-soil slab interface is very complex and must be handled with care. This has been done in NEWTMC.

Discrete System Approximations
The usual ways of approximately solving these nonlinear parabolic equations are those of finite differences (5-7) or finite elements (8)(9)(10)(11)(12)(13). In either case we trade off the continuous system. This discrete approximation is formally analogous to a compartment or conceptualized system. For the remainder of this work we make the following additional assumptions: (1) the mass and heat transport phenomena in the soil are one-dimensional only, vertically, with the positive orientation being downward; and (2) the soil "slab" is partitioned into a finite number of adjoining soil layers (NSLYRS) with each layer being specified by its characteristic thickness Azi and its induced volume Vi. Note that Yi = AAzi. Figure 4 illustrates the partitioning of the soil "column" or "slab" together with the major mass and

Moisture Field in the Soil Slab
Beginning with Eq. (13) the assumption of an average mass of water all concentrated at the geometric center of the kth slab (layer) of soil bounded by planes Zk-l and Zk is now made so that we have the ordinary differential equation approximation for layer k given as:  (14), and p"t is three to four orders of magnitude less than pw for most field situations (clearly controllable in the TMC), the second term in Eq. (31) is neglected, leaving the reasonable approximation for total water in the kth layer as: XMk -AAZkPWO(Zk -1/2) (32) from which we readily obtain: where Zk l + AZk k = 1, 2,.., NSLYRS Z= 0 Interestingly enough, had we started with the conceptual approach we would have arrived at exactly this same ordinary differential equation. q,, is now substituted into Eq. (27) with due regard to the fact that only part of the surfaces at Zk-X and Zk contribute to the flow of moisture in each of the two respective phases. When this substitution is made, the equation is: Equation (27) is a mass balance equation for water. The instantaneous time rate of change of water in the kth soil layer is equal to the difference between the inflow (or outflow) of moisture from the k-1 layer and the outflow (or inflow) of moisture from the k layer to the k + 1 layer. Since the quantities pwV,z and PWVV are components of the total moisture flux vector they can represent either positive (downward transport) or negative (upward transport) flow of moisture across any and all of the artificially introduced boundaries (soil layer interfaces). Thus, for example, the flow of moisture from the bottom of soil layer k-1 becomes the source of moisture into the top of layer k, and so forth, for the other layers. At each soil layer interface the principles of continuity of flux in the normal component of moisture flow and continuity in moisture tension (pressure field) are applicable.

Heat Field in the Soil Slab
Proceeding exactly as before in the moisture field case, we consider the flow of heat into and out of the two bounding planes zk-1 and Zk (one-dimensional heat transport assumed), and the total heat content of soil layer k is placed at the geometric center of the layer. This heat content is given the interpretation of being an average heat value for the kth layer. The space integrals indicated in Eq. We give Eq. (35) the equivalent heat balance interpretation as we gave mass balance interpretation to Eq. (28). That is, the instantaneous time rate of change of heat in the kth soil slab is equal to the difference between the inflow (outflow) of total heat (conduction + convection + latent heat) from (to) soil layer Zk l and the outflow (inflow) of total heat to (from) soil layer Zk +1 A The heat field (average) and the temperature field in layer k are related via the expression:

Simulations of Moisture and Chemical Transport in Unsaturated Soils
To examine representative responses in an unsaturated soil, four simulation runs were made. Detailed documentations describing the computer code are available in Lindstrom and Piver (2). The program is of modular design using a main program to call different subroutines in a prescribed sequence with COMMON statements used to link the different subroutines. The sequence and order of call are shown in Figure 6 requirements are greatly reduced over those that would be required for a program using a finite element or finite difference approximation for the coupled set of nonlinear partial differential equations. Simulations of 40 days using complex soil structures and weather conditions were solved on a VAX/DEC computer.

System Variables
Values for parameters that describe and drive transport of moisture, heat and chemical in unsaturated soils are given in Tables 3-6. The variables that remained constant for the four simulations are listed in Tables 3-5, and those that were varied from run to run are listed in Table 6. Because transport in the soil column is sensitive to the external driving forces occurring in the air chamber, two different rain schedules were used for the 40-day simulations. For both runs, the duration of light and dark, however, were maintained at 14 hr of sunlight followed by 10 hr of darkness. The intensity of sunlight was held constant at 4.16 cal/cm2-hr, a value corresponding to bright sunlight at 10 A.M., at a latitude of 45°N on June 21 (the latitude of Corvallis, OR). In the first rain schedule, rain accumulation over a period of 10 hr was equal to an inch of rain. In the 40-day time interval, five rain events occurred at a frequency of one every eight days, a rain schedule similar to temperate climates (with this schedule, average annual rainfall is 42 in.). In the second rainfall schedule, rain fell for 2 hr at five different times during a time interval of 40 days, conditions comparable to an extended dry season in a temperate climate. The daylight and darkness schedule and the two rainfall schedules are given in Table 3.  Tables 4 and 5 list values for physical-chemical properties of the air-soil system that were held con-   Table 4, properties that are not varied with depth are listed. Included in this group are basic physical-chemical properties and parameters that describe the moisture tension function and the water conductivity function. In Table 5, properties that are functions of depth are listed. These include layer thickness, porosity, soil constituents, tortuosity, and rate constants for chemical and biological reactions. It should be noted that rate constants for all removal processes are between 10-5 and 10_6 hr-1, values indicative of nonreactive chemicals. Table 6 lists variables that are changed from run to run. Variables in this group include dispersivity and adsorptivity coefficients for different soil components. Precise data for these variables were not available, and several values were examined in lieu of either using stochastic methods to select values as a function of  4.90 x 10-9 cal/hr-cm2-°K 4 system variables or approximating these variables with power law functions of moisture content.
Hydraulic conductivity and moisture tension have been represented as power law functions of moisture content. The algorithms of NEWTMC have been constructed to allow the hydraulic conductivity and moisture tension functions to be functions of depth in the soil, making it possible to develop empirical representations for these variables that are sensitive to the characteristics of the soil at that position. The power law representation for moisture tension is given as: where a, and P,B are empirically determined constants; Osat(Z) is the moisture content at saturation for each soil layer; and 0(z) is the moisture content of the soil in layer k. The power law representation for hydraulic conductivity is given as: where K(Osat) is the saturated hydraulic conductivity in the kth soil layer and -y,c is an empirically determined constant. Other representations for hydraulic conductivity are given by Gardner (19). For the soils presented in the illustrative simulations, a clay silt loam was used 22 (45)    (15). The length of the soil column was 716 cm; the soil porosity was assumed to be uniform, and the power law constants for the moisture tension and hydraulic conductivity expressions were: Oxo = 2.5 cm; P36 = 6.0; -y,o = 6.5; and K(Osat) = 0.3 cm/hr (46) Discussion In Figures 7 and 8, the hydraulic and chemical responses during single rain events for the two rainfall schedules given in Table 3 and the soil characteristics of  Table 5 are shown. In both instances, initially the soils had been allowed to drain to their field capacities. Figure 7 shows the hydraulic response of the soil column. The responses are similar except in degree of moisture wetting and redistribution with time and position in the soil column. This is due to the amount of water applied to the surface. When the rain stops, the combined effects of evaporation at the surface and gravity drainage quickly drain the water from the top layers of the soil column.
The effects of the amount of rainfall during each event on chemical transport are shown in Figure 8. The mechanisms of many of the processes in the chemical mass balance [Eq. (38)] are not well understood, making it necessary to estimate values for many of the chemical transport variables. Variables that make up this group include reaction rate constants, dispersion coefficients, dispersivities, and adsorption coefficients. The transport of moisture in unsaturated soils, however, has been extensively examined and hydrological data exist for many different soils and for many conditions of saturation and partial saturation (15)(16)(17)(18). Because the transport of moisture, heat, and chemicals in soils are coupled phenomena and because the transport of moisture is behaving correctly as shown in Figure 7, there is a high probability that the chemical response is also correct.
The intensity and duration of a rainfall event has its most important effect in the upper soil layers. In Figure  8 for the 10-hr rainfall schedule, the concentration varies over a wide range during and after the rain event. There are changes in concentration in the upper soil layers during the 2-hr rainfall schedule, but they are not as extensive. Much of the variation is due to the manner in which concentration is being reported, as micrograms per cubic centimeter of soil water. The more water there is, the greater the dilution. However, at a depth of only several centimeters below the surface, these surface effects damp out quickly, and the responses for both rainfall schedules are similar except penetration is not as great for the soil subjected to the shorter rainfall schedule.
Variations in chemical transport properties and their influence on temporal variations in chemical concentra-   FIGURE 8. Changes in chemical concentration as a function of depth in the soil column for (A) the 10-hr rain event producing 1.0 in. of rain, and (B) for the 2-hr rain event producing 0.2 in. of rain. Soil properties are given in Table 5; dispersivity = 1.0 cm; adsorption coefficient = 400 cm soil water/cm soil.  Table 3; soil properties are given in Table 5; dispersivity = 1.0 cm; adsorption coefficient = 400 cm soil water/cm soil. tion proffles are shown in Figures 9-12. The values for these properties that were changed from run to run are given in Table 6. In each figure, both the heavy and light rainfall schedules are represented, and the chemical concentration profiles are for unreactive chemical substances. The chemical was added to the top soil layer at the start of the simulation. The effects of these physical-chemical properties on chemical transport are very apparent. When adsorptivities are high and dispersivities low, there is very little movement of the chemical past the top layers of the soil column. On the other hand, for low values of the adsorption coefficient and higher values of the dispersivity coefficient, penetration into the soil column is much greater; however, penetration is still strongly influenced by the amount and frequency of rainfall. Frequency and intensity of rainfall still appears to be the major process that drives the transport of chemical in the soil column.
Because the processes of adsorption and removal of chemicals by biological and chemical means are not well understood, methods of estimating appropriate values for parameters that describe these processes are important topics of investigation. With regard to chemical transport, variations in adsorptivity coefficients for a chemical onto different soil components significantly effect movement of that chemical in unsaturated soils. Other equally unknown mechanisms of actions are removal of chemicals by biological and chemical processes. In the examples presented here, the chemicals were assumed to have a low degree of reactivity both by biological and chemical processes. In addition, it was assumed that the rate of removal could be approximated by first-order reaction-rate kinetics. There is no guarantee that this is an accurate description of these processes. In such a complex structure such as soil, the mechanics of adsorption and reaction on surfaces with variable adsorptivities have not been investigated. To Changes in chemical concentration with depth and time for a high water conductivity, low adsorptivity soil and for (A) an average temperate climate rainfall schedule and (B) for an extended dry period rainfall schedule. Rain schedules are defined in Table 3; soil properties are given in Table 5; dispersivity = 10.0 cm; adsorption coefficient = 400 cm soil water/cm soil. . Changes in chemical concentration with depths and time for a low water conductivity, high adsorptivity soil and for (A) an average temperate climate rainfall schedule and (B) for an extended dry period rainfall schedule. Rain schedules are defined in Table 3; soil properties are given in Table 5; dispersivity = 1.0 cm; adsorption coefficient = 4000 cm soil water/cm soil. 10 10 FIGURE 12. Changes in chemical concentration with depth and time for a high water conductivity, high adsorptivity soil for (A) an average temperate climate rainfall schedule and (B) for an extended dry period rainfall schedule. Rain schedules are defined in Table 3; soil properties are given in Table 5; dispersivity = 10.0 cm; adsorption coefficient = 4000 cm soil water/cm soil.
initiate remedial action for purifying contaminated aquifers, such data are essential. Degradation of chemicals by microbial metabolism in soils is largely uninvestigated for industrial chemicals. Much attention has been given to the role of microorganisms in nutrient cycling but very little is understood about the ability of microorganisms to metabolize industrial chemicals either as primary or secondary carbon sources. Until recently the unsaturated zone of the soil below the root zone was considered to be biologically dead (20). Recent evidence for a limited number of compounds, however, suggests that the number of microorganisms at all depths in the soil does not decline but remains fairly uniform (21). The characteristics of these microbial populations, however, are not well understood.
One final note about this model relates to its ability to simulate multicomponent chemical transport. The transport of chemicals from landfills involves the migration of many chemicals. In this situation, terms must be added to account for interactions among the different chemicals and heat and mass transfer equations must be written for each component. In landfills organic solvents can be present resulting in partitioning of chemical components between different liquid phases. Because this situation exists, a single chemical model acting as a surrogate for a mixture of chemicals may only be appropriate when the chemicals are present in low concentrations. If concentrations are low, interactions of all types are usually negligible and the model is a very accurate method for simulating multi-chemical transport in unsaturated soils under a variety of soil conditions and external driving functions.
The initial modeling portions of this research were carried out by F. T. Lindstrom during a six-month IPA with the Corvallis Environmental Research Lab (CERL) division of the U.S. EPA. Initial software development for NEWTMC was supported in part by contract number 808864 for the Toxics and Hazardous Materials Branch, CERL. E I CL-LLI C: