published_or_final_version / Civil Engineering / Doctoral / Doctor of Philosophy

published_or_final_version / Applied Geosciences / Master / Master of Science

A method was developed for assessing the vulnerability of groundwater to contamination from contaminant sources at the soil surface, using a numerical groundwater flow model linked to a digital map database. The method was applied using the pcARC/INFO Geographic Information System (GIS) to input, store, and manipulate base maps, resulting in a database of digital maps for the alluvial aquifer system in the Willamette Valley of western Oregon. Digital elevation maps were created by digitizing topographic maps of land surface (1:250,000 scale), water surface, and the base of the Tertiary-Quaternary sedimentary deposits (1:500,000 scales). Soil association and aquifer unit maps digitized from 1:500,000 scale map sheets were also used. Data were extracted from ARC/INFO to the SURFER software package to create a 3-D surface model for each of the digital elevation maps. An ARC/INFO point coverage was then used to store and overlay these surfaces, allowing the creation of maps of depth to water, saturated thickness, and water table gradient. These data became the input to a numerical finite element groundwater flow model. The model solves a dual formulation problem for the potential function and the stream function to calculate the time-of-travel for water to flow from the surface to the water table and laterally for 100 meters as an index of groundwater vulnerability. A cluster analysis is used to condense the data and form a training data set for a multiple regression model. The regression model is fit to the results of the finite element model with an R-squared of greater than 0.96. The simpler regression model is then used for mapping travel times for the entire study area. When properly calibrated against the finite element model and when combined with the digital map database and Geographic Information System (GIS) procedures described, the regression model can be conveniently used to assess the vulnerability of groundwater to contamination over large areas. / Graduation date: 1989

This dissertation develops an analytical-numerical approach to deterministically predict the space-time evolution of concentrations in heterogeneous geologic media conditioned on measurements of hydraulic conductivities (transmissivities) and/or hydraulic heads. Based on the new conditional Eulerian-Lagrangian transport theory by Neuman, we solve the conditional transport problem analytically at early time, and express it in pseudo-Fickian form at late time. The stochastically derived deterministic pseudo-Fickian mean concentration equation involves a conditional, space-time dependent dispersion tensor. The latter not only depends on properties of the medium and the velocity but also on the available information, and can be evaluated numerically along mean "particle" trajectories. The transport equation lends itself to accurate solution by standard Galerkin finite elements on a relatively coarse grid. This approach allows computing without using Monte Carlo simulation and explicitly the following: Concentration variance/covariance (uncertainty), origin of detected contaminant and associated uncertainty, mass flow rate across a "compliance surface", cumulative mass release and travel time probability distribution across this surface, uncertainty associated with the latter, second spatial moment of conditional mean plume about its center of mass, conditional mean second spatial moment of actual plume about its center of mass, conditional co-variance of plume center of mass, and effect of non-Gaussian velocity distribution. This approach can also account for uncertainty in initial mass and/or concentration when predicting the future evolution of a plume, whereas almost all existing stochastic models of solute transport assume the initial state to be known with certainty. We illustrate this approach by considering deterministic and uncertain instantaneous point and nonpoint sources in a two-dimensional domain with a mildly fluctuating, statistically homogeneous, lognormal transmissivity field. We take the unconditional mean velocity to be uniform, but allow conditioning on log transmissivity and hydraulic head data. Conditioning renders the velocity field statistically nonhomogeneous with reduced variances and correlation scales, renders the predicted plume irregular and non-Gaussian, and generally reduces both predictive dispersion and uncertainty.

The development of an open pit phosphate mine in 1965 near Aurora, North Carolina, required groundwater withdrawals in excess of 50 million gallons per day for pressure relief from the underlying confined Castle Hayne aquifer. The effects of pumping from this limestone aquifer were widespread, extending over an area of more than 2,000 square miles in the first year before the pressure cone reached a stabilized condition. Salt water encroachment by lateral movement from downdip in the aquifer and by downward leakage from the Pamlico Estuary and Sound was possible if prolonged pumping were permitted. A joint project between state and industry in 1970-71 focused on analyzing the effects of five years of pumping for the purpose of making projections of future conditions of groundwater quality for continued and expanded groundwater development. The present study describes the method of flow net analysis used in the joint project to determine quantitatively the values of aquitard vertical permeability and aquifer transmissivity. These values provided the rational basis for making projections of groundwater movement and quality changes associated with this movement, A hydrologic projection model, developed in the present study, integrates the quality and volumetric flow of vertical leakage with that of lateral flow. Projection analyses are presented for two hypothetical cases of chloride distribution changed by continued pumping and for chloride changes associated with three different pumping regimes in the subject area.

The main aim of this study is to develop a suitable method for the calibration and validation of mathematical models of large and complex aquifer systems. Since the calibration procedure depends on the nature of the model to be calibrated and since many kinds of models are used for groundwater, the question of model choice is broached first. Various aquifer models are critically reviewed and a table to compare them as to their capabilities and limitations is set up. The need for a general calibration method for models in which the flow is represented by partial differential equations is identified from this table. The calibration problem is formulated in the general mathematical framework as the inverse problem. Five types of inverse problems that exist in modeling aquifers by partial differential equations are identified. These are, to determine (1) parameters, (2) initial conditions, (3) boundary conditions, (4) inputs, and (5) a mixture of the above. Various methods to solve these inverse problems are reviewed, including those from fields other than hydrology. A new direct method to solve the inverse problem (DIMSIP) is then developed. Basically, this method consists of transforming the partial differential equations of flow to algebraic equations by substituting in them the values of the various derivatives of the dependent variable (which may be hydraulic pressure, chemical concentration or temperature). The parameters are then obtained by formulating the problem in a nonlinear optimization framework. The method of sequential unconstrained minimization is used. Spline functions are used to evaluate the derivatives of the dependent variable. Splines are functions defined by piecewise polynomial arcs in such a way that derivatives up to and including the order one less than the degree of polynomials used are continuous everywhere. The natural cubic splines used in this study have the additional property of minimum curvature which is analogous to minimum energy surface. These and the derivative preserving properties of splines make them an excellent tool for approximating the dependent variable surfaces in groundwater flow problems. Applications of the method to both a test situation as well as to real-world data are given. It is shown that the method evaluates the parameters, boundary conditions and inputs; that is, solves inverse problem type V. General conditions of heterogeneity and anisotropy can be evaluated. However, the method is not applicable to steady flows and has the limitation that flow models in which the parameters are functions of the dependent variable cannot be calibrated. In addition, at least one of the parameters has to be preassigned a value. A discussion of uncertainties in calibration procedures is given. The related problems of model validation and sampling of aquifers are also discussed.

Theoretical and experimental studies of the chemical and physical factors which affect molecular diffusion of dissolved substances from fractures into a tuffaceous rock matrix have been made on rocks from G Tunnel and Yucca Mountain at the Nevada Test Site (NT8). Although a number of physical/chemical processes may cause nonadvective transport of dissolved species from fractures into the tuff matrix, diffusion in these rocks is controlled by the composition of the groundwater through multicomponent effects and several rock properties. The effective molecular diffusion coefficient of a particular species in the tuff can be related to its free aqueous diffusion coefficient by Dₑ = θ(m)(α/τ²)D₀ where bm is matrix porosity, α is the constrictivity, and τ is the tortuosity. The porosities of the samples studied ranged from 0.1 to 0.4. The parameter (α/τ²) ranged from 0.1 to 0.3, and effective matrix dif— fusion coefficients were measured to be between 2 to 17. x 10⁻⁷ cm²/s for sodium halides and sodium pentafluorobenzoate. Total porosity was found to be the principle factor accounting for the variation in effective diffusion coefficients. The constrictivity— tortuosity factor was found to have a fair correlation with the median pore diameters measured by mercury intrusion. Measurements of bulk rock electrical impedance changes with frequency indicate that the constrictivity factor, a, has a maximum value of 0.8 to 1, but may be smaller. If the larger values are correct, then the diffusion paths in tuff are more tortuous than in granular media. The diffusion coefficient matrix computed for various tracers in J-13 well water from the NTS indicates coupling of the diffusion fluxes of all ionic species. Multicomponent diffusion is a second order effect, however, which does not significantly affect experimental results. The results of a bench—scale fracture flow experiment revealed that the transport of ionic tracers (SCN ⁻ and pentafluorobenzoate) was affected by diffusion into the tuff matrix. The transport of a particulate tracer did not appear to be affected by diffusion.

Multiple aquifer systems similar to those that exist at coal gasification sites are complicated groundwater situations. In these types of systems, the aquifers are separated by aquitards through which interaction between aquifers can occur. The movement of the products of combustion into the coal seam and adjacent aquifers is a serious problem of interest. This dissertation presents two-dimensional finite element models for water flow and water quality in multiple aquifer systems. These models can be applied for general problems as well as the problems associated with the burned cavities in coal gasification sites. The Galerkin weightedresidual method is used in both models. Eight-noded isoparametric elements are used. Spatial numerical integration is performed using Gaussian quadrature. A weighted finite difference scheme is used, in both of them, for time integration. The two models are written in FORTRAN V for the CDC CYBER 175. They are applicable to one- or two-dimensional problems involving steady-state or transient flow. Each aquifer can have different initial conditions and boundary conditions. Boundary conditions, pumping rates, and the recharge can be specified as a function of time. The output of the flow program-nodal heads and velocity components is used as an input to the quality program. The numerical models were validated for simple problems that have available analytical solutions.

Subsurface fluids may travel along paths having a minimum permeabilility greater than the effective permeability of the rock. This may have an important impact on contaminant migration. A stochastic approach related to percolation theory is advanced to address the question of what is the probability that a high permeability path extends across a given volume of the subsurface. The answer is sought numerically through subdividing the volume of interest into a three-dimensional grid of elements and assigning a random permeability to each element. Four permeability processes are considered: 1) Stationary with independence between grid elements; 2) Stationary and autocorrelated; 3) Nonstationary due to conditioning on measured values; and 4) Random rock volume included in grid. The results utilizing data from fractured granites suggest that in large grids, at least one path having a minimum permeability in excess of the "effective" rock permeability will cross the grid. Inclusion of autocorrelation causes an increase in the expected value of the minimum permeability of such a path. It also results in a significantly increased variance of this permeability. Conditioning on field permeabilities reduces the variance of this value over that obtained by unconditional, correlated simulation, but still produces a variance greater than that obtained when independence was assumed. When conditioning is performed, the mean of the minimum permeabilities along these paths is dependent on the principal axis of the path. Finally, including a random rock volume by allowing the length of the grid to be random increases the variance of the minimum permeability.

Methodologies for estimating hydraulic and solute transport properties of unsaturated, fractured rock are developed. The methodologies are applied to networks of discrete fractures for the purpose of estimating steady fluid flow rates and breakthrough curves of entrained solutes. The formulations employ the boundary integral method to discretize the outer rim of each fracture and to solve a two dimensional flow equation within fracture planes. A three dimensional variant of the two dimensional boundary integral method is used to calculate flow through a permeable matrix with embedded permeable fractures. Exterior and interior surfaces are discretized using boundary elements to account for flow between fractures and the matrix, and between the matrix and fractures and the exterior boundaries. Synthetic fracture networks are created using planar fractures of finite areal extent embedded within a three dimensional rock matrix for the purpose of performing sensitivity studies of network hydraulic conductivity with respect to geometric parameters, such as fracture orientation and density. Results of the sensitivity studies show that: (1) The global hydraulic conductivity is linearly dependent on the product of fracture transmissivity and density for fractures of which fully penetrate the rock volume; (2) The effect of correlation between fracture length and transmissivity is to increase the global hydraulic conductivity; and (3) Results using a three dimensional coupled fracture— matrix flow regime compare favorably with analytic results. Flow through variably saturated fracture networks is modeled by assuming a constant capillary head within individual fractures. A free surface is found using an iterative procedure which locates nodal points at the intersection of constant total head and pressure head contours. The simulated free surface compares favorably with an approximate analytic solution and with laboratory results. Simulations indicate the presence of zones of water under both positive and negative pressure, as well as regions of air—filled voids. Travel times and breakthrough curves are determined by integrating the inverse velocity over a streamline, and then summing over all streamlines. For the fracture network examined, travel times decrease with decreasing fracture saturation. The effects of retardation and matrix diffusion are also examined.