Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous 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 which ensures fluid conservation both locally and globally. Because the relative permeability and saturation curves for fractures may be highly nonlinear, and in strong contrast to those of the matrix, the robust Newton-Raphson iteration method is implemented according to the efficient procedure of Kropinski (1990) and Forsyth and Simpson (1991) to solve the variably-saturated flow equation. Upstream weighting of the relative permeabilities is used to yield a monotone solution that lies in the physical range and adaptive time stepping further enhances the efficiency of the solution process. A time-marching Galerkin finite-element technique is 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. Verification examples are presented along with illustrative problems that demonstrate the complexity of variably-saturated flow and solute transport in fractured systems.