Software

MAFIC Software Information (continued)

Matrix Equation Solver

MAFIC uses an incomplete-Cholesky conjugate-gradient (ICCG) algorithm to solve the matrix equations of flow. This approach was first suggested by Meijerink and van der Vorst (1977), and is based on the conjugate-gradient method for efficient location of the minimum of a function (Hestetnes and Steifel, 1952). The incomplete Cholesky technique creates an approximate Cholesky decomposition of the flow matrix by decomposing only the non-zero terms in the matrix. For typical MAFIC meshes, the matrix is very sparse, and storing only the non-zero terms results in a substantial savings in memory requirements.

The approximate Cholesky decomposition is used to pre-condition the conjugate-gradient algorithm, and results in quite rapid convergence compared to PSOR or L-U solution algorithms. There is thus a substantial savings in memory, and a useful savings in CPU time using this algorithm. It has been shown that the conjugate gradient algorithm will always converge, and will require no more iterations than the number of degrees of freedom in the system.

The incomplete Cholesky factorization does not work unless the matrix is diagonally dominant. When finite elements containing obtuse angles or extreme aspect ratios are used, the matrix may not be diagonally dominant. In this case, a modification suggested by Kershaw (1978) is employed in MAFIC: if the diagonal term of the decomposed matrix becomes negative (so the necessary square root cannot be taken), the original diagonal is increased by a small arbitrary amount. (It is set to BOOST times the amount which is subsequently subtracted from it, where BOOST is a user-input parameter with a default value of 2.0.) This adjustment, made in subroutine ICDCMP, stabilizes the decomposition and does not impede the convergence of the ICCG algorithm.

The user specifies a parameter, TOL, which defines the conditions necessary for convergence of the algorithm. TOL, which typically has a value of 0.0001, is used in two tests of convergence. The largest absolute change in head at any node must be less than TOL times the largest nodal head minus the smallest nodal head. The largest unbalanced (residual) nodal flow rate must be less than TOL times the sum of the absolute values of all net nodal flows (which is approximately twice the total throughflow of the system).

In general, the algorithm appears to converge to a high level of accuracy using a number of iterations roughly equal to 0.1 times the number of degrees of freedom (the number of nodes which are not at a fixed head or on a group-flux boundary).

For steady-state problems, the flow matrix becomes indefinite (and hence unsolvable) if there are any finite elements which are not connected, directly or indirectly, to a fixed-head node. MAFIC uses a subroutine called XFLOAT which identifies groups of elements/nodes which "float" in this manner. As each "floater" is identified, its nodes are forced to have a fixed head equal to 999.99. Removal of any floaters is performed automatically by MAFIC, and sometimes results in a significantly smaller system of equations to solve. Removal of floaters may be disabled by setting input parameter RFLOAT false.

Modeling of Solute Transport

As discussed in the introduction, solute transport is currently simulated primarily within the fracture elements. Solute transport in the rock matrix is modelled only within the context of matrix diffusion between fractures and generic rock matrix blocks.

Solute transport in one dimension can be described by the equation (Bear, 1972):

(2-12)

where: V = the true groundwater velocity (L/T)
d_L = the longitudinal dispersivity (L)
D = the coefficient of molecular diffusion (L² /T)
x = the distance coordinate (L)
c = solute concentration (dimensionless)
t = time (T)
Q_s = external source (T^-1)
C_s = concentration of solute at external source (dimensionless)

The coefficient of molecular diffusion is ignored in the current formulation. This is a reasonable assumption for flow through fracture networks since fluid velocities are normally sufficiently large that V × d_L >> D.

The solution to Equation (2-12) subject to the boundary conditions:

c(0,0) = C_o
c(x,0) = 0.0
V(x,t) = a constant

is given by Bear (1972) as:

(2-13)

The form of Equation (2-13) is identical to that of a normal distribution with a mean of Vt (the convective travel distance) and a variance of 2d_LVt. Hence dispersion can be modeled as a random stochastic process about the mean convective travel distance. The dispersive travel distance in the direction of motion (longitudinal dispersion) is given by Prickett et al. (1981):

(2-14)

where: X_d = longitudinal dispersive travel distance (L)
X_c = convective travel distance: V×t (L)
d_L = lateral dispersivity (L)
N_L = a normal random variable with a mean of zero and a variance of one
(dimensionless)

A second dispersive component (transverse dispersion) orthogonal to the first dispersive component, is defined in a similar manner:

(2-15)

where: Y_d = transverse dispersive travel distance (L)
d_T = transverse dispersivity (L)
N_T = a normal random variable with a mean of zero and a variance of one (dimensionless)

Description of Particle Tracking Approach

Solute transport is simulated using a particle tracking approach. Similar versions of this approach have previously been applied to two-dimensional discrete fracture models (Smith and Schwartz, 1984; Hull et al., 1987). The particle tracking approach represents the concentration of solutes within the solvent (normally water) by a finite number of discrete particles of equal mass. Each particle represents a fraction of the total mass of solute present. Particles are introduced into the modeled region at specified solute sources and removed at encountered sink nodes. At each time step, particles are moved according to a deterministic convective component and a stochastic dispersive component. The convective component, Vt, is proportional and parallel to the velocity vector at the current particle location. The dispersive component, given by Equations (2-14) and (2-15), is proportional to the square root of the convective component.

Average concentrations are determined at the end of each timestep by summing the particle mass within each fracture element and dividing by the mass of fluid within the element. The equation for determining particle concentrations within a specific element is:

(2-16)

where: N_p = number of particles in element (dimensionless)
M_p = mass of a single particle (M)
P_w = density of solvent (M/L³)
A_e = area of fracture element (L²)
b_e = fracture aperture (L)

[more]