anisotropy

Current research on fracture hydrology at Lawrence Berkeley Laboratory is focused on two scales: the single fracture and fracture networks. In the area of single frac tures, we have established a research program aimed at characterizing the geometry of the void spaces. This in volves the use of metal casts and mercury injection to elu cidate details of the complex geometry that result as a fracture closes under stress. With this geometry, the Bow regime is being examined using an approach based on Navier-Stokes equations. Three techniques are being used to investigate the geometry of fracture networks: (1) statistical analysis of observed field data, (2) geophysical techniques for seeing into the rock, and (3) prediction of fracture patterns through the application of rock mechanics. Statistical analysis has centered on developing a parent-daughter model to represent fracture swarms. Techniques for comparing the results of seismic tomography with the fracture properties controlling flow are being investigated. Shear zones have become a major focus of our efforts because they appear to control the hydrology of many sites. New models for. fluid flow include the development of codes to analyze flow and transport in a network of channelized frac tures.


INTRODUCTION
Modeling the behavior of fluid flow in fractured rock is largely a problem 'of defining the geometry of the fractures. The calculation of flow, once the geometry is known, is a far simpler problem. However, a major problem that must be faced is that of determining the location of all the fractures and describing their properties. We can at most hope to describe the fractures and derive the properties of the network in a _stochastic sense. To create a stochastic model requires that one have a conceptual model that can reveal how the fractures are located in space and what relationship they have to each other. Given the conceptual model, the field data can be used to derive the parameters of the model, and the parameterized model can then be used to estimate the behavior of the field site.
Recent research on fracture hydrology at Lawrence Berkeley Laboratory has been focused on two different scales: the single fracture and fracture networks. Investigations of single fractures have been concerned with the problem of identifying the details of the geometry of the flow spaces within the fracture. Once the geometry of the flow paths is known, one can examine the nature of this flow regime. The primary problem in dealing with fracture networks is characterizing the location of the fractures and this involves a statistical sample from which a stochastic model can be inferred.

SINGLE FRACTURES
A variety of laboratory and analytical techniques are being developed to define the po're space in a fracture (Myer et al., 1985;Myer et al., 1986). Concurrently, a numerical method for solving the Navier-Stokes equations in a complex, multiply connected domain is being developed by Muralidhar and Long (1987a, 198ib). This is an attempt to investigate the theoretical relationship between flow and fracture void geometry and will also be used to analyze flow in voids of a real fracture.

Void Geometry
Efforts to characterize the void geometry of a fracture have involved several experimental techniques (Brown and Scholtz, 1985;Tsang and Tsang, 1987;Gentier, 1986). A new experimental technique that is described in Myer et al. (1986) involves injection of a low melting point metal into the fracture under controlled stress and fluid pressure conditions. Metal injections can be performed on fractures under different normal loads. When the metal cools, the fractures are taken apart and photographed at several magnifications. The metal is brighter and shows up as white in the photograph. The contacts are black. Some metal sticks to both surfaces, so the phot05 must be superimposed to obtain the true distribution of contacts and voids. An example from the work of Myer et al. (1986) of such a photo prepared from a fracture in crystalline rock at 85 ~1Pa effective stress is shown in Figure 1. The photographs can be analyzed with. the help of an image analyzer. The image analyzer ... \ r v resolves each photo into 512 by 512 pixels which are either black or white according to a cut-off value for brightness. Choosing an appropriate cut-off value involves an error. Superimposing the pictures also involves an error because it requires aligning the pictures. We are studying the patterns formed by the pixels for different cutoffs and different alignments. One can use a two-point correlation function to demonstrate insensitivity to alignment. The two-point correlation function is the conditional probability of finding a pixel on the flow path at an arbitrary distance from a pixel on the network of conducting pixels (i.e., voids). Thus, the two-point correlation function is a measure of connectivity. Preliminary results show that at some resolution the twopoint correlation function does not depend on the way the photographs are superimposed. Obviously, changing the cutoff changes the picture of the percentage of contact area and connectivity.
. The aperture distribution for the voids can also be studied by injecting mercury at different pressures from which one obtains a pore-size distribution function. Pore sizes found in this way can be assigned to the white pixels. However, we have little information about where a particular size of void should be assigned. Most probably voids are spatially correlated. That is, voids are smaller near con tacts or near other small voids and larger near large voids. Performing the metal injection at different states of stress should help to determine this.

Flow Calculations
Given the void geometry at some state of stress, one can proceed to examine the hydraulic characteristics of the fracture. From the point of view of fluid mechanics, a fracture is a.n irregular three-dimensional channel with irregularly spaced and sized obstacles. In such channels we expect four possible How regimes. If the fra'cture is sufficiently open, the How behavior may be essentially the same as that of How bet~een parallel plates. As the fracture closes the distribution of the contact areas will become a controlling factor in determining the How characteristics. If the contact area is not well distributed but occurs in sparse clumps, the How may be non-linear, even at low Reynolds numbers. Thus Hux will not be proportional to gradient. If the contact area is finely distributed, the flow will be very much like that of porous media. Flux will be linearly related to gradient and potential theory can be used to describe head distributions. If the fracture is sufficiently closed, the How will only occur in tortuous channels. In this case, flow in the fracture plane becomes essen tially a two-dimensional network of one dimensional conductors. We have developed a three-dimensional grid generation procedure which maps the interior of a complex region onto the interior of a rectangular parallelepiped (~1uralidhar and Long,198ia,198ib). In the transformed region (com pu tational space) the boundaries lie along lines which are parallel to one of the new coordinate axes. The flow computation is performed in this newly defined region. and an inverse transformation rule extracts values of flow velocities and pressure in the original system of coordinates. The formulation has been developed in three dimensions. suitable for flow analysis in a single rock fracture.
Because of the transformed spatial coordinates, one must also transform the governing equations for flow before they can be solved. In doing this. each term in the three-dimensional governing equations becomes 81 terms. Thus. we have traded complexity in geometry for a large number of unknowns. The equations governing grid generation in three dimensions, including a method to improve orthogonality of the generated coordinates, and the techniques for solving flow in the transformed region have been fully derived in Muralidhar and Long {1987a}.
A real fracture is also expected to have a number of contact areas serving as obstacles to Buid Bow. In this approach, the contact areas are mapped to rectangles, of about the same total area. An example of this procedure is shown in Figure 2. The actual contact area is circular ( Figure 2a) and in the transformed plane, it is a square (Figure 2b). The dimensionless duct is a size 6 x 6 x 2 dimensionless units and the obstacle is of size 1.2 x 1.2 x 2. The duct is closed at the top and the sides. The grid in the physical domain is gradually distorted to accommodate the boundary shape which is not coincident with a rectangular mesh pattern.
Once the grid generation has been performed, the flow field can be solved. Flow enters on the left side and leaves the right side, in response to a prescribed pressure drop. Results have been obtained for specified Reynolds numbers, Re, where Re is based on the pressure drop in the channel. Figure 3 shows a velocity vector plot at the mid-plane of the duct, at Re = 25. A symmetric roll pattern is clearly visible on the rear side of the obstruction, and is a consequence of assuming the Bow to be steady. As the Reynolds number is raised, the steady state computation retains the symmetric form of the recirculation region, except that it grows in size. It is known from experiments in flow past cylinders and spheres that the symmetric bubble is stable only for low Re. At higher flow rates, the Bow becomes unsteady. with a trailing vortex area being produced behind the cylinder, which in turn destroys symmetry. The onset of unsteadiness in Figure 3 would have to be determined from a fully transient calculation. More than one obstacle in the flow field provides a stabilizing inBuence to the Bow. Many obstacles would lead to the development of Darcy flow.
The next step will be to analyze the data set from the real fracture as described above and use it to create an input file for the grid generation and flow calculation models. One can then obtain the geometry for different states of stress and observe the effect on changes in the fi·ow regime.
3 STUDIES OF NETWORK GEOMETRY.

Parent-Daughter Model
The primary problem in developing a model of fractured rock from field data is to characterize the location of the fractures. One approach to the problem is to treat the fractures that can be observed as a statistical sample and try to infer a stochastic model for the fracture system (Long and Billaux, 1986). This is the approach that we have taken.
Fractures are often observed to occur in swarms or zones. In order to model fractures in swarms, one can use a statistical description of this type of pattern called a "parent-daughter" model. In this model, the fracture swarms, or daughters are nucleated around seeds, or parents. The location of the parents may be purely random, i.e., a fixed rate Poisson process, or there may be a regional variation in the density of the parents. Once the parents have been determined, the daughters are found at some distance from the parents where this distribution may, for example, be Gaussian. The location of parents, location of daughters relative to their parents, and number of daughters are taken to be independent random quantities. Thus, implementation of the model requires that we know the distribution of: (a) the de~ sity of the parents, (b) the number of'daughters per parent, (c) the daughters around the parents.
There are two major difficulties in obtaining these distributions from field data. The first problem is that it is not necessarily clear to which parent a fracture daughter belongs. This is further complicated by regional changes in the density of fractures. The second problem i.s that the parents and daughters are distributed in three-dimensional space but we can only observe traces. The centers of the traces are not the same as the centers of the fractures and the parent is a point that does not necessarily lie in the plane of observation.
The first problem requires that we propose a parent-daughter model and see how well it fits. Given a trial example of a three-dimensional parent-daughter model with regional variation in density, Deverly {1984} has worked out the statistics which would be observed on a plane. J. P. Chiles (Long' et al.,198i) modified this derivation to account for the case where the daughters are disc shaped fractur.es rather than points. The result of his work is that when the distributions for the parameters of the regionalized parent-daughter model can be assumed. we can calculate the theoretical variogram of fracture trace density on tht> plane. Then, one can compare the theoretical variogram with that derived from the fracture mapping of drift walls. Thus, we are now able to select an appropriate regionalized parent-daughter model by trial and error.
Application of this statistical model to data from the Fanay-Augeres mine in France is now underway. Fanay-Augeres is a uranium mine owned by COGEMA, a French public company, where various investigations for nuclear waste storage have been carried out. A drift in the mine is the focus of our invt>stigation. Two sections of the drift have been mapped. One section of the drift is wet and was used as the site of an inflow measurement. The other section is dry. Models of both sections are being developed in an effort to determine why one section is wet and the other dry.
The rock surrounding the wet section of the drift is penetrated by a reverse fault called the Faille de Recette. There is strong evidence that in fact this fault and not the fracture network, is the reason that this section is wet (J. L. Bles, personal communication, 1987). If we do not see any major difference in the fracture networks of the two zones, this will be strong evidence that the shear zone is dominating the hydrology.
The other two approaches that we are taking to characterize the geometry of fracture networks are both concerned with characterizing such shear zones because they apparently are so important. In this regard, the seismology approach is particularly useful for shear zone identification. Further, we are trying to develop shear zone genesis models as a basis for fracture hydrology models.

Seismic Tomography Studies for Fracture Detection and Characterization
Geophysical methods provide one of the best potential tools for obtaining the geometric data that are needed in analyzing fracture network behavior. In particular. we are investigating the use of three-component vertical seismic profiling (VSP) for fracture detection and characterization. It is important to record the three components of the signals because of the phenomenon of shear wave splitting (Majer et aI., 1986). When a shear wave passes through a fracture, the shear component splits into two parts, one faster than the other, (Crampin, 1978(Crampin, , 1981(Crampin, , 1984a(Crampin, , 1984b(Crampin, , 1985. The wave parallel to the fracture is unaffected by the fracture. The wave perpendicular to the fracture is attenuated because the stiffness across the fracture is less than that of the rock matrix. Schoenberg (1980Schoenberg ( . 1983 has shown that the ratio of the velocity of a seismic wave perpendicular and parallel to a set of discontinuities is a function of the spacing of the discontinuities as well as the stiffness. Knowing the stiffness and the velocity anisotropy, one should be able to determine the average fracture spacing, or density. Alternatively, given independent information on fracture density, one could determine the' fracture stiffness and hopefully relate this stiffness to actual fracture properties. One might then discriminate between filled and open fractures and the effect on hydraulic conductivity, Such information can be invaluable in constructing a fracture ne'twork for hydrologic analysis of a particular field site.
In order to calibrate geophysical tools as indicators of hydrologic behavior. one needs a basis for comparing the hydrologic response to the geophysical response, Th is is not a straigh tforward problem.
Seismologic profiles are done by creating a source and monitoring the response at a receiver. The ray path between the two points can be estimated and thus the tomographic inversion of the seismic data is possible, The inversion yields a pixel map of the average rock properties which affect the propagation of the wave. However, even with a complete hydrologic description of the fracture network, it is difficult to predict the resulting behavior of a seismic wave passing through the rock. The reasons for this are that we do not know the relationship between fracture stiffness and permeability, and there are many calculational problems associated with solving the wave equation at each fracture boundary. Seismologists can invert their field data but, they are not proficient at forward modeling.
Hydrologic data, on the other hand, are notoriously hard to invert. Fluid can pumped into one zone in a well and the response monitored in another zone. However, particularly in fractured rock, the path the fluid has taken between the signal and the monitoring points is difficult to know exactly. It could have traveled' through anyone of a large number of paths. Thus, inversions tend to be highly nonunique. Some recent progress on this problem has been made by Doyen (1987), who uses maximum entropy methods to identify likely network parameters. Conditional ge08tatistical simulation also holds promise for limiting the number of possibilities. On the forward modeling side, hydrology is well advanced. Given the properties of the fracture network, there are a number of techniques for modeling flow and transport in two-and three-dimensions (Long et al., 1982;Long, 1983;Long and \Vitherspoon, 1985;Dershowitz, 1984;Schwartz et al., 1983;Smith and Schwartz, 1984). There is good forward modeling and poor inverse modeling in the field of hydrology, whereas the reverse is true in the field of seismology.
In order to find a basis for comparing the methods of hydrology and seismology,we have begun a numerical study of seismic wave transmission through the same twodimensional fracture networks used for the hydrologic parameter studies. In this way, the seismic response can be compared to the hydrologic response. In the programs that we have developed, the effects of fracture stiffness and bulk rock properties have been incorporated into the model (Myer et al., 1985). We are also studying the relationship between fracture void geometry and both hydraulic behavior and stiffness. Such an approach will produce a fracture network with values of fracture stiffness and permeability that both conform to the same geometry. A simple example is shown in Figures 4, 5 and 6. Figure 4 shows a two-dimensional fracture network representing a highly fractured shear zone in a less fractured matrix. We have calculated the behavior of seismic rays as they pass through the network in a cross-hole tomographic pattern from two horizontal wells assumed to be at the top and bottom of the pattern. Figure 5 shows the attenuation field resulting from inversion of the synthetic seismic data. This particular inversion was based only on amplitude data. The amplitude of the seismic wave is attenuated as it passes through the fracture. More sophisticated inversions could include such phenomena as shear wave splitting.
For the same fracture network, we have calculated How in different directions. These How data have been resolved into directional permeability res.ults and plotted in Figure 6 as a permeability ellipse. This plot shows the anisotropy in permeability that results from the shear zone. Although this is avery simple example, it demonstrates the principle we would like to use in comparing seismic results with hydrologic behavior. It has become apparent at many sites throughout the world that shear zones dominate the hydrology of the site. Such highly fractured zones conduct most o( the water and must be characterized for site evaluation. In general these shear zones are not simple hydrologic features. They are complex zones including shear fractures, extension fractures, brecciated zones, stress dissolution features, various fillings and alteration materials. Such a highly heterogeneous zone might not be well characterized as an equivalent continuum. This is especially the case when one considers transport of contaminants in a zone with this kind of heterogeneity. For this reason, we have begun to look more closely at the genesis of shear zones. Specialists in geology and rock mechanics have studied shear zones, and several geologic models are available (Segal and Pollard, 1983;Bles and Feuga, 1986). New work by Kemeney (1987) has opened the possibility of modeling the growth of shear zones through consideration of failure criteria. Geologic and mechanical information has not been explicitly included in fracture network models for hydrology. Kemeney's aim is to find a way to generate fractures that make sense from the point of view of mechanics and geology and use these models to study hydrology. A new three-dimensional model called CHANGE has been developed by D. Billaux (personal communication, 1986) which introduces channelized flow into our threedimensional disc model for fracture flow. This model simply introduces a twodimensional network of conductors that are imbedded in the fracture discs. The channels in each fracture can be generated indepe.ndently using a two-dimensional fracture mesh generator. Any distribution can be used. The channels can be parallel or have any random distribution of orientation in the plane of the fracture. The length and conductivity distributions of the channels can likewise be controlled. Figure 7a shows an example of a three-dimensional network of disc-shaped fractures with sub-networks of channels in the fractures. In Figure 7b the outlines of the discs have been eliminated, showing only the conductive channels in space.
The fracture intersections are treated as a separate class of channels in this model. Currently, we give these intersection channels a much larger conductivity than the channels in the fractures as this seems to be geologically most reasonable. However, these intersections can be treated otherwise if that is desired. The generation of channels within a fracture disc can be confined to any sub-area of the disc. This effectively eliminates the constraint of having disc shaped fractures. The disc shape is only used to find the intersections between fractures. After that, the shape of the disc is irrelevant.
The most important feature of this model is that it essentially reduces the threedimensional network problem to a quasi two-dimensional problem. Flow through the network is calculated with exactly the same algorithm used to calculate flow through a two-dimensional network. We call it "quasi two-dimensional" because the issue of connectivity in the network is both two-and three-dimensional. There is a twodimensional connectivity in the plane of each fracture, and three-dimensional connectivity between the fractures.

Two-Dimensional Transport
A new model for transport in two-dimensional fracture networks has been developed by Karasaki (1987). The code can be run on meshes generated with CHANGE as described above. Thus, we are able to calculate conservative tracer transport in three-dimensional networks of channelized fractures. This code incorporates a Lagrangian and Eulerian scheme with adaptive gridding and is capable of handling transport under steady or transient conditions. Previous models. such as that developed by Endo et aI., (1984) have treated only the mechanical component of the transport by assuming a Poisseuille distribution of Bow within each fracture. In other words, Endo assumed that Bow in the fractures was essentially like that between parallel plates. Models developed by Long and Shimo (1987) or Robinson (1982) assumed most of the dispersion was due to the geometry of the network. These models neglect diffusion and dispersion within the fracture. In the new model of Karasaki, transport in the fracture is assumed to behave like transport in porous media. In other words, there is an assumption of one-dimensional, advection-dispersion in each fracture channel. This approach is indicated by the fact· that fractures are not necessarily open channels. Contact area and infillings will make Bow in the fracture behave more like that of a porous medium.
The code is designed to model tracer tests in fractured rock which are being conducted either under ambient (steady) or forced (transient) Bow conditions. One expected use of the code is to help determine likely fracture patterns from a combination of flow and tracer tests. Borehole fracture mapping can be used to obtain fracture network parameters such as fracture frequency and orientation distribution. Hydraulic conductivity of individual fractures and the network might be estimated from small and large scale pumping tests, respectively. This kind of data would be needed to hypothesize network parameters and simulate tracer tests in the models. By comparing the tracer test behavior in the field with that of the various fracture models, we propose to learn more about the likely properties of actual fracture networks.