Introduction

PHREEQC, the U.S. Geological Survey’s aqueous geochemistry model, is capable of modeling dual-porosity reactive transport in one dimension. A primary row of connected cells is used to simulate multispecies advective-dispersive transport, while columns of adjacent, transverse cells can be used to model diffusion. This conceptualization can be used for the simulation of reactive transport through fractured porous media, for example.

As a demonstration project, I have been working on creating visualizations for supergene copper ore enrichment in a block of fractured rock. Supergene ores form when pyrite and other sulfide-bearing primary minerals above the water table are oxidized, creating an acidic percolate that carries copper and other dissolved minerals to the water table, where local reducing conditions lead to the formation of secondary sulfide mineral phases. In fractured rock setting, these processes will naturally occur within fracture zones, where percolation of recharge water is favored over the less-permeable surrounding matrix. Constructing such a simulation, as well as processing the results for visualization, can be greatly facilitated by means of python scripting.

Pre-processing

Setting up a PHREEQC simulation for a dual-porosity system involves employing the TRANSPORT and MIX keyword blocks (refer to the PHREEQC User’s Guide at ftp://brrftp.cr.usgs.gov/pub/charlton/phreeqc/Phreeqc_3_2013_manual.pdf for more details). The TRANSPORT block establishes the parameters for the advective-dispersive transport model through the advective/mobile-water cells, whereas the MIX keywords describe the extent of diffusive mixing between adjacent immobile-water cells, and at the mobile-immobile cell interfaces, per time step.

Manually computing mixing parameter values and setting up the MIX keyword blocks would be extremely tedious, particularly if parameter values are to be adjusted and a model re-run to refine the simulation. Fortunately, it is a straightforward task to automate writing keyword blocks for both TRANSPORT and MIX in response to more fundamental user input parameterization (e.g., cell-spacing, time step size) so that appropriate model values are assigned consistently between mobile/advective-diffusive and immobile/diffusive cells.

A second set of keyword block entails the initial condition definitions for the solid phase, such as EQUILIBRIUM_PHASES and KINETICS, which may be spatially variable. Because PHREEQC employs a specific numbering scheme to define cells in a dual-porosity configuration (mobile cells are listed first, then immobile cells in columns parallel to the mobile cells), keyword definitions must follow consistent numbering scheme. Again, writing these blocks can be easily scripted once basic model parameters are defined.
A python script for populating the each of the above keywords is provided in my GitHub repository. Any additional keyword blocks used to fully define a PHREEQC model run (e.g., RATES blocks, for kinetic rate laws) can be edited manually. It should be noted that PHREEQC can be called directly from within python (using, for example, the IPhreeqcCOM add-in for Windows-based computers), but a solution that involves assistance with writing the external input files read by PHREEQC is more general.

The code for accomplishing the above pre-processing steps can be found in my GitHub repository.

Post-processing

A second issue for my example problem is post-processing the results, particularly for a such a fractured system residing within a three-dimensional block of rock, with fractures characterized by arbitrary orientations. PHREEQC can write transport simulation results in a tabular text file via the SELECTED_OUTPUT keyword block. For a dual-porosity simulation, this output will usually include concentrations of various aqueous species, pH, pe, and masses of equilibrium phases and those species determined via kinetics expressions. Concentrations are specified by cell number, by distance along the mobile-cell column, and by time step. For a particular model run, this output corresponds to a set of two-dimensional model results that must be generalized to three dimensions via extrusion along a third axis, multiaxial rotation, and translation parallel to the background grid imposed by the full matrix rock block.

Projection of the PHREEQC results onto a three-dimensional lattice, with multiple fractures with arbitrary orientations, as well as a background matrix chemical composition, becomes an exercise in analytic geometry that is amenable to manipulation with python/numpy (and would be very tedious otherwise with spreadsheets, etc.). An approach I have taken is to first extrude the PHREEQC results along the y-axis (from the x-z plane referenced output) and to determine the normal vector to the user-specified strike and dip for each fracture. Each fracture can then be turned in three-space via a rotation matrix so that the fracture normals are aligned with the normal surface to the y-z plane, (1, 0 ,0):

class Fracture:

    def __init__(self, fileName, nz, nx, ny, dz, dx, dy, w, phi, theta, datum, y0, grid_x_intercept, nzDivide):
        # fracture properties (and simulation results, as data frames)
        self.name = fileName.split('.')[0]                          # fracture name
        self.nz = nz                                                # discretization along fracture
        self.nx = nx                                                # discretization transverse to fracture
        self.ny = ny                                                # extruded discertization, along y-direction
        self.dz = dz
        self.dx = dx
        self.dy = dy                                                # fracture depth (for gridding, independent of PHREEQC)
        self.w = w                                                  # width of advective zone
        self.phi = phi*pi/180.                                      # dip below horizontal (input as degrees)
        self.theta = theta*pi/180.                                  # strike, clockwise from N (input as degrees)
        self.datum = datum                                          # z value at coordinate system rotation point
        self.y0 = y0                                                # y-axis extrude starting location
        self.grid_x_intercept = grid_x_intercept                    # grid x-axis fracture crossing location
        self.nzDivide = nzDivide                                    # number of gridding subdivisions along advection axis
        self.phreeqc = read_csv(fileName, delim_whitespace=True)    # read associated PHREEQC transport output file
        self.fracVec = self.NormPlane()                             # normal vector to fracture plane
        self.rotMatrix = self.RotationMatrix()
        self.eqn = self.EqnPlane()                                  # equation of fracture plane

    def WithinFrac(self, x, y, z):
        # determine if point (x, y, z) is within the fracture zone
        interior = (self.DistPtPlane(x, y, z) <span 				data-mce-type="bookmark" 				id="mce_SELREST_start" 				data-mce-style="overflow:hidden;line-height:0" 				style="overflow:hidden;line-height:0" 			></span><= (self.w + self.nx*self.dx))
        return interior

    def DistPtPlane(self, x, y, z):
        # normal distance between point (x, y, z) and fracture plane
        # see College Calculus with Analytic Geometry, 3rd Edition, Protter and Morrey (1982); p. 410
        num = abs(self.eqn[0]*x + self.eqn[1]*y + self.eqn[2]*z + self.eqn[3])
        den = sqrt(self.eqn[0]**2 + self.eqn[1]**2 + self.eqn[2]**2)
        return num/den

    def EqnPlane(self):
        # determine equation of fracture plane in the form Ax + By + Cz + D = 0
        # see http://mathonline.wikidot.com/point-normal-form-of-a-plane
        A = self.fracVec[0]
        B = self.fracVec[1]
        C = self.fracVec[2]
        D = -A*self.grid_x_intercept
        return array([A, B, C, D])

    def NormPlane(self):
        # calculate the vector normal to a plane of strike theta and dip phi in 3-space
        # see http://eqseis.geosc.psu.edu/~cammon/HTML/UsingMATLAB/PDF/ML1%20FaultNormals.pdf
        nx = sin(self.phi) * cos(self.theta)
        ny = -sin(self.phi) * sin(self.theta)
        nz = -cos(self.phi)
        return array([nx, ny, nz])

    def RotationMatrix(self):
        # methodology: https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d
        v = cross(refVec, self.fracVec)
        c = dot(refVec, self.fracVec)
        vx = array([[0., -v[2], v[1]], [v[2], 0., -v[0]], [-v[1], v[0], 0.]])
        R = identity(3) + vx + linalg.matrix_power(vx, 2) * 1./(1. + c)
        return R

Next, for my initial approach to visualization, the objective was to populate a lattice, with each lattice point representing a chemical or mineral concentration that can be interpolated across three dimensions. Three-dimensional data visualization tools, such as Golden Software’s Voxler, operate using this approach. However, generation of the lattice to represent both the background matrix composition and the different conditions within the fracture zone(s) requires some algorithm design considerations. These include:

  1. A large number of points might be needed, and so vectorizing operations via numpy can be advantageous with respect to execution speed.
  2. The (rotated) fracture zone(s) are included in the lattice, making it irregular; any extension of the fractures beyond the block of rock represented by the matrix can be simply clipped off. However, background matrix composition lattice points cannot be permitted within the fracture zone(s) as this would adversely impact the resultant chemical/mineral distribution interpolation.

The solution to the first issue is simply to use numpy methods to quickly spawn a regular lattice (as opposed to for-loops and if-then statements to explicitly avoid placing background matrix points in the fracture zones):

def FillGrid(self, fracs):
     # gross fill-in of background grid points
     indxSets = mgrid[0:self.nx, 0:self.ny, 0:self.nz].T.reshape(-1, 3)
     indxCols = transpose(indxSets)
     x = self.x0 + (indxCols[0] + 0.5) * self.dx
     y = self.y0 + (indxCols[1] + 0.5) * self.dy
     z = self.z0 + (indxCols[2] + 0.5) * self.dz
     df = DataFrame({'x': x, 'y': y, 'z': z})
     # remove points within fracture zones
     df['fracZone'] = 0
     for frac in fracs: df['fracZone'] += frac.WithinFrac(x, y, z)
     df = df[df['fracZone']==0]
     df.drop(['fracZone'], axis=1, inplace=True)
     # add background quantities
     for i, item in enumerate(self.bkgName):
          df[item] = self.bkgNum[i]
     return df

The second issue can be addressed by converting the lattice into a pandas dataframe and then deleting those points that fall within fracture zone(s) according to a point-to-plane distance formula (see method, listed under Code Snippet #1).

The full code for processing the PHREEQC SELECTED_OUPUT results for each dual-porosity simulation can be found in my GitHub repository, along with some simple instructions for running the script.

 

lattice
Lattice generate by post-processing script. Dark points represent fracture zones, while light points represent the background rock matrix.

Results

In the supergene copper demonstration scenario, two fracture zones cut across a block of rock, 2.5 m x 2.5 m on a side and 15 m in height. Both fracture zones strike (clockwise from north) and 30 degrees and dip 70 degrees below horizontal. Advective-transport through the fracture on the left is approximately five times faster than the adjacent fracture to the right. Flow is from top to bottom; the water table (i.e., transition from oxidizing to reducing conditions) occurs approximately two-thirds of the way along the column, from the top.

Results are shown for several decades of elapsed simulation time:

 

chalcopyrite
Chalcopyrite dissolution within fracture zones above the water table.
brochiantite
Brochantite alteration produced within fracture zones above the water table.

 

dolomite
Modeled distribution of Ca, Mg-carbonate phase (as disordered dolomite); note leaching from oxidized zone and precipitation in reduce zone. As this is just a demonstration, slow precipitation kinetics of dolomite are ignored.
goethite
Goethite precipitation in fracture zones in the oxidized portion of model domain above the water table.

 

gypsum
Modeled gypsum precipitation, a reaction product of pyrite and chalcopyrite dissolution, within fracture zones.
secondaryCuSulfides
Zone of secondary copper ore enrichment in fracture zones beneath water table. Phases include bornite, chalcocite, and covellite.