This thesis developed a numerical model capable of simulating the processes involved in subsurface remediation. It is a general purpose model, which provides flexibility in input parameters and configurations. It includes the processes of two-dimensional saturated flow; Freundlich and Langmuir sorption isotherms; single, double, and competitive Monod kinetic models; and intermediate toxicity. In addition, it can simultaneously simulate the movement of multiple solutes.

A literature search of existing numerical models discovered that no single existing model provided all of the processes described above. Therefore, this thesis developed a new model. It chose construction of a new model over modifying an existing model for various reasons. The two primary reasons were to avoid programming errors due to unfamiliarity with existing code, and the choice of C++ as the programming language.

Chapter 3 discussed the construction of the model. A finite difference algorithm solved the fluid flow and was easy to integrate with the other numerical method. The Integrated Operator-Splitting method provided easy implementation and flexibility with respect to adding different types of biological and abiotic reactions. The model uses an explicit forward time backwards space finite difference scheme called the Mean Concentration Slope method used by Poulsen (1994). The Linear Integrated method simulated the reaction portion of the IOS method. In addition, the chapter developed maximum and minimum time step criteria for the model.

Chapter 4 evaluated model performance both quantitatively and qualitatively. Due to the lack of simple analytical solutions to advection and dispersion coupled with kinetic reactions, the chapter evaluated advection and dispersion processes independently from reaction processes. The Relative Sum of Squared Error quantified the differences between numerical and analytical solutions. The chapter reported it as a function of Peclet and Courant numbers. The error in advection and dispersion processes was highly dependent on the choice of Courant and Peclet numbers. This has implications on the choice of timestep and grid spacing.

Chapter 5 illustrated the use of the model as an aid in site remediation design. The model investigated different remedial strategies in a fictitious one-dimensional aquifer. The system was highly sensitive to TCE dissolution rate and relatively insensitive to grid mesh size. Model results suggested that pump-and-treat is the only feasible strategy for the spill. Biodegradation appeared to be unfeasible.

Some suggestions for model improvement and expansion surfaced during the use and evaluation of this model. First, execution speed should be increased. The numerical code needs to be optimized for speed. Secondly, the model is currently limited to two-dimensional simulations, and should be expanded to three. Another improvement is adding transient flow, and multiple fluid flow to the model. These improvements would help increase the utility of this program as a modeling tool.