Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media

The presence of fractures in a low-permeability medium can greatly influence ground-water flow paths and contaminant migration rates. Fractures represent preferential pathways along which a solute can migrate rapidly; however, contaminant diffusion from the fractures to the porous matrix can significantly reduce migration rates along the fractures. Because it is critical to include such processes as matrix diffusion and advection along the fractures, and perhaps advection in the matrix if it is sufficiently permeable, efficient numerical techniques must be used to simulate large-scale flow and transport problems in such media. A discrete fracture, saturated-unsaturated numerical model is developed where the porous matrix is represented in three dimensions and fractures are represented by two dimensional planes. This allows a fully three-dimensional description of the fracture network connectivity. Solute advection and diffusion in the porous matrix are also directly accounted for. The variably-saturated flow equation is discretized in space using a control volume finite element technique. Because the relative permeability and saturation curves for fractures may be highly non-linear, and in strong contrast to those of the matrix, the robust Newton-Raphson iteration method has been selected to solve the variably-saturated flow equation. Upstream weighting of relative permeabilities and adaptive time stepping further enhance the efficiency of the solution process. The Laplace Transform Galerkin technique for cases where flow is at steady-state or a standard Galerkin finite element technique for transient flow situations are used to discretize the solute transport equation. Although the methodology is developed in a finite element framework, a finite difference discretization for both groundwater flow and solute transport can be mimicked through a manipulation of the influence coefficient technique. The use of an ILU-preconditioned ORTHOMIN solver permits the fast solution of matrix equations having tens to hundreds of thousands of unknowns. A fully three-dimensional simulation of groundwater flow and solute transport through a fractured aquitard and into an underlying aquifer reveals that the interpretation of field observations of the watertable position, hydraulic heads and solute concentration in the fractured aquitard-aquifer system can be exceedingly difficult without knowledge of the nature and location of the fractures. Unlike predictions based on simpler two-dimensional models, a fully three-dimensional treatment shows that even the plume that forms in the underlying aquifer will he highly irregular in space because of the complex pattern of fracture traces at the interface between the fractured aquitard and the aquifer. Under variably-saturated conditions, a further complexity arises because the fractures can act as either preferential flow paths or flow barriers depending on the contrasts between the hydraulic properties and constitutive relations of the fractures and the porous matrix.