Processes of Interest

At each scale discussed so far, we may look, categorically, at three coupled, namely, microbiological, hydrogeological, and geochemical processes in the subsurface. Here, we will emphasize the representation of these processes in "averaged" forms rather than the details of hydrodynamics and geochemical reactions at microscopic level. So you are forewarned to bear in mind that, although the equations that we are going to see are frequently cited in the literatures, they at best represent averaged or upscaled processes in a strict sense. For our subsequent discussion of modeling the Melton Branch Watershed, we will use several finite element computer codes written by my colleagues and myself, and by Professor George Yeh at the Pennsylvania State University. The equations that you are going to see are implemented in these codes, which have been tested extensively in previous applications.

To lots of us, science is a detective work; yet to the others, it may look more like piecing puzzles together. So, bear with me when I’m trying to make sense to you all the components that we need in a watershed scale contaminant transport model.

Taking a reductionist viewpoint, one may want to use the Navier-Stokes equations to model the momentum and fluxes in the pore space. The capillary forces will need to be examined in each of the pores that are partially saturated. For large-scale modeling, it is very unlikely that we can afford this amount of calculation. The capacity of our computational facilities today is not going to accommodate this resolution of the hydrodynamic processes. And, more importantly, we are not going to be able to afford characterization of the pore geometry. For watershed scale modeling, it is more likely that we are going to use the Darcy’s equation and a mass conservation principle for the transport of water and solutes. But beside this necessary choice of an averaged process, what can we do with the additional information that the soils are structured and macroporous and the fact that we have some measurement of fracture spacing and aperture at hand? First of all, back in 1953, Sir George Taylor derived the famous "virtual coefficient of diffusion," later known as the hydrodynamic dispersion coefficient in the quantitative hydrogeology community, for a circular pipe:

, (1)

where a is the radius of the pipe, U the mean velocity of fluid, and D₀ the molecular diffusion coefficient of a solute in the fluid. The fractures and microfractures I showed you earlier, however, are in planar or plate form. In 1959, Professor Wooding at the Emmanuel College in Cambridge arrived at a similar expression for this kind of geometry:

, (2)

where l is the aperture of a fracture. In reality, we are looking at solute molecules with a certain thermodynamic nature of moving randomly across streamlines at space while at each location they are sampling the velocity of the carrying fluid molecules. This process along causes the hydrodynamic dispersion that the above two equations are trying to describe. More recently, Prof. Howard Brenner and his colleagues at MIT have extended these previous works to include periodically varying pore structures and has gained a lot of insight in the upscaling of the dispersive process. What I am trying to point out here is that the information of pore geometry and the upscaling relationship done in the past have enabled us to better integrate information at the microscopic scale to macroscopic scale applications. There is a catch here. The walls of the pipe in Tayler and Wooding’s relationships are impermeable. At mesoscopic scale there is at least a diffusive component that my result in diffusion of solutes into or out of the matrix material that is permeable. In fact, lots of literatures have shown that the matrix material is not only permeable but also consists of most of the pore volume and surface space that may store lots of contaminants in the soils and groundwater. The movement of contaminants in and out of the matrix pores are know as matrix diffusion among geologists and mass transfer among soil scientists. If one chooses to use an "equivalent" porous medium (EPM) approach to represent a soil, this mass transfer component is usually "averaged" out. I believe this is one of the reasons why an EPM model is not likely to perform as well as a mobile-immobile model or a multiple-pore-region model.

Now, if the matrix material is to be viewed as an EPM, I have to assert that, given the width of the matrix chunk or the distance between two adjacent fractures, we can use the following advective/diffusive or CD equation to calculate the mass flux across the fracture-matrix interface:

, (3)

where t is time, x is distance into the matrix, c_m is the concentration of a solute in the soil matrix, and q_m, v_m, D_m are the water content, fluid velocity, and dispersion coefficient in the matrix, respectively. The mass transport process can also be represented similarly:

, (4)

where c_f is solute concentration in the fracture, and q_f, v_f, D_f are the water content, fluid velocity, and dispersion coefficient in the fracture, respectively. The mass flux across the fracture-matrix interface will be taken as the boundary conditions of the matrix blocks. And, this is the catch here! Writing the equations this way and relying on the interface "internal" boundary condition to couple the two equations, one will need pretty detailed information of matrix block size or fracture spacing data. This approach is intuitive, easy to conceptualize, and has been taken into the implementation of FRAC3D and FRACTRAN, computer codes developed at the University of Waterloo by Professor Sudicky and his colleagues. FRAC3D and FRACTRAN have been very successful in modeling a wide variety of hydrogeological problems, but one would suggest that it is not always possible to implement detail fracture networks in a watershed scale model. So, it would be nice that eqs. (3) and (4) can be explicitly coupled through a mass transfer term as the following:

, (5a)

, (5b)

, (5c)

where the interface flux can be q_fm can be calculated, instead, using the following equation:

. (6)

where e^f_fm and e^t_fm are the advective and diffusive mass transfer coefficients, respectively. The advantage of such formulation includes the advection component in the matrix transport equation and the upscaling of mass transfer processes between pore domains, without the cost of characterizing pore geometry. The disadvantage is that we introduce two additional parameters, the advective mass transfer coefficient and the diffusive mass transfer coefficient, to the equations. The trick now becomes how do we incorporate the interface flux effected by a distribution of matrix blocks, as shown in this rendering of bedding plane fractures on one end of an undisturbed soil column, into eq.(6). The key is the interregion CD equation and the distribution of matrix block sizes. Because we can theoretically calculate the dispersivity in the soil matrix (Wooding’s equation) and the other parameters for the interregion CD equation can virtually be either estimated or calculated using experiment data, the mass transfer coefficients in eq.(6) can be determined without much uncertainty. Assuming that the fractures extend from the top to the bottom of a soil column, one can determine the matrix size distribution as shown here. From the sizes of matrix blocks (as plates) and using the CD equation, one can in turn calculate the mass transfer coefficients in eq.(6). The velocity terms in eq.(5), in the case of soil column tracer injection, can usually be measured and calculated using appropriate measurements of hydraulic conductivities. For real world, three-dimensional modeling, determining the hydraulic conductivities of a multiregion soil is more involved, which we would not have time to cover here today. This, perceivably, is the disadvantage of taking a multiple-pore-region approach. There are, however, at least two references here: my coauthor Glenn Wilson, Phil Jardine and I published a paper in the SSSAJ in 1992 on the subject of how to determine a three-pore-region hydraulic conductivity functions and water retention functions; and Dr. Smettem and Kirkby at CISRO, Australia, also published a paper in the Journal of Hydrology in 1990 that included a method used to determine a two-region hydraulic conductivity curve for their aggregated soils. More recently this year (1997), Dr. Mohanty and his colleagues have a paper in the WRR that deals also with this subject.

In addition to groundwater flow and solute transport, the biological and chemical processes in the groundwater are also very important to large-scale water resources and water quality management. The chemical speciation calculations could be based on the assumption of equilibrium or kinetic geochemical systems, depending on the flow rates and chemical residence times. Processes considered in the codes we are using, HYDROGEOCHEM developed by Professor George Yeh at the Pennsylvania State University, include aqueous complexation, adsorption/desorption, ion-exchange, precipitation/dissolution, oxidation/reduction, and acid/base reactions. HYDROGEOCHEM consists of a kinetic-equilibrium chemical speciation module, KEMOD, that calculates chemical concentrations, based on mole balance and mass action laws, in the fluid and solid phase of the soil. For example, the walls of the fractures in the saprolite of the C horizon soils I mentioned earlier are usually coated with Fe- and Mn-oxides. It is very likely that we have a DOC rich water in the upper A and B horizons that dissolves and mobilizes the irons on the surface of solids that have been weathered for years to exposed the mineral. The translocated irons are then moved down the soil horizon by infiltrating rainwater and absorbed onto the less weathered surface of fracture walls in the deeper C horizon. Mathematically, these reactions, for example, an aqueous phase complexation and an aqueous-solid phase adsorption/desorption reaction, can be represented as:

, (7a)

for kinetic species, and for equilibrium species

, (7b)

where T_j is concentration of species j, Q_j is external source, G_b,j is the production or reduction rate of the species through biogeochemical reactions, a_i,j is the stoichiometric coefficient, c_i the component species i, P_b,j is a derived species j, and N_a is the number of aqueous component species. Writing the equilibrium reaction in the form of (7b) is for the convenience of computational representation, where only the concentrations of the component species will be used as state variables. This is what we adopted in HYDROGEOCHEM, where the total analytical concentrations of component species are calculated in the transport step. Concentrations of derived species are calculated using the concentrations of relevant component species. So, we can imagine two small loops within a big loop in the modeling of hydrogeochemistry. The first small loop takes care of solute transport of every single component species and the second takes care of the speciation of the entire geochemical system. The outer, bigger loop couples the transport and geochemistry processes into a hydrogeochemistry system. From a computational viewpoint, one may think of each location in a model as consisting of a beaker that contains a local batch and the batch is interconnected by tubes with its neighbors with a certain flux rate in between. This conceptualization leads to the well-known, highly favorable understanding of nearly perfect parallel system in parallel computing. We will discuss this idea further in the application section.

In HYDROGEOCHEM, Microbial growth and metabolism are modeled using a modified Monod kinetics as the following:

, (8a)

, (8b)

where X_i is the concentration of microbe species i, m_i is the specific growth rate, l_i is the metabolic potential, n_i is the decay/endogenous rate, m₀ is the maximum growth rate, K_j is the substrate half-saturation constant, and N_b is the number of species participating in the specific microbiological reaction. In addition to microbial growth, we are also implementing a microbial deposition mechanism that uses an empirical filtration constant relationship.

The biogeochemistry itself defines a fairly complicated and computationally time consuming system even for today’s high speed workstations. Our capability to model field scale biohydrogeochemistry, currently, depends a lot on the capacity of our computers. We are trying to address this problem from the viewpoint of parallel computing, which I’ll be discussing in just a few minutes. Before that, let us return to the subject of watershed scale solute transport modeling and allow me to recount a research project initiated ten years ago by Bob Luxmoore at ORNL and my coauthors Glenn and Phil. The project has been completed in 1992, but we hope that we will be able to revamp the site, largely because we have had so many experiments conducted and so much data collected at the site, and because we have run numerous solute transport models, trying to understand the watershed-wide transport characteristics.