
International
Journal of Rock Mechanics and Mining Sciences Volume 39, Issue 7, October 2002, Pages 887904 
Result list  previous < 23 of 56 > next 
 

 
Estimation of REV size and threedimensional hydraulic conductivity tensor for a fractured rock mass through a single well packer test and discrete fracture fluid flow modeling
M. Wang^{a},
P. H. S. W. Kulatilake^{}^{, }^{}^{, }^{a},
J. Um^{a}
and J. Narvaiz^{b}
^{a} Department of Mining and Geological
Engineering, University of Arizona, 1235 E. North Campus Drive, Mines Bldg, No.
12, Rm 229. P.O. 21 0012, Tucson, AZ 85721–0012, USA
^{b} Metropolitan Water District of Southern California,
Los Angeles, CA 90054, USA
Accepted 6 June 2002. Available online 8
October 2002.
Abstract
A new methodology is presented to determine the representative elementary volume (REV) size and threedimensional (3D) hydraulic conductivity tensor for a fractured rock mass. First, a 3D stochastic fracture network model was built and validated for a gneissic rock mass based on the fracture data mapped from scanline surveys at the site. This validated fracture network model was combined with the fracture data observed on a borehole to generate a stochasticdeterministic fracture network system in a cubic block around each packer test conducted at a different depth region in the same borehole. Each packer test was simulated numerically applying a developed discrete fracture fluid flow model to estimate the influenced region or effective range for the packer test. A cubic block of size 18 m, with the packer test interval of length about 6.5 m located at the centre of this block, was found to be suitable to represent the influenced region. Using this block size, the average flow rate per unit hydraulic gradient (defined as the transmissivity multiplied by mean width of flow paths) field for fractures was calibrated at different depth regions around the borehole by numerically simulating the packer tests conducted at different depth regions. The average flow rate per unit hydraulic gradient of the fractures that intersect the borehole was considered to be quite different to the average flow rate per unit hydraulic gradient of the fractures that do not intersect the borehole. A relation was developed to quantify the ratio between these two parameters. By studying the directional hydraulic conductivity behaviour of different cubic block sizes having the validated stochastic fracture network and calibrated hydraulic parameters, a REV for the hydraulic behaviour of the rock mass was estimated to be a block size of 15 m. The hydraulic conductivity tensor in 3D computed through regression analysis using the calculated directional hydraulic conductivity values in many directions was found to be significantly anisotropic. The principal directions of the hydraulic conductivity tensor were found to be agreeable with the existing fracture system of the site. Further, the geometric hydraulic conductivity calculated was found to be comparable to the hydraulic conductivity estimated through the radial flow assumption in continuum porous media.
Author Keywords: Fracture mapping; Single well packer tests;
Representative elementary volume; Hydraulic conductivity tensor; Discrete fracture fluid
flow modeling; Hydrogeology; Calibration; Field study; Numerical modeling
Article Outline
1. Introduction
Currently used methods for estimating threedimensional (3D) hydraulic conductivity tensor using either aquifer pumping test or packer test data are based on the assumption that the groundwater is flowing through a geologic continuum [1, 2, 3 and 4]. These methods can generally be applied when wells penetrate the porous geological media such as alluvial deposits, but may have limited applicability when the geological medium is dominated by a fracture system that has well defined fracture sets. The porous, continuum media assumption is based on an average flow within a representative elementary volume (REV) [1 and 2]. Even though the REV is small for porous media, for fractured rock masses it can be very large or in some cases may not exist [5]. If the REV does not exist, or is larger than the distance between the pumping and observation wells or the packer test hole and observation wells, it is not appropriate to use the equivalent continuum approach to analyze the aquifer pumping or packer test data for fractured rock masses. Before applying an equivalent continuum approach for a rock mass, one should investigate the REV for the hydraulic behaviour of the rock mass. The REV for the hydraulic behaviour is defined as the size beyond which the rock mass hydraulic conductivity tensor remains the same. For sizes larger than or equal to the REV size, the equivalent continuum approach can be applied for the rock mass with a secondrank, symmetric, positivedefinite hydraulic conductivity tensor. Estimation of the REV for hydraulic behaviour of a highly heterogeneous, anisotropic rock mass requires monitoring of groundwater level of a significant number of observation wells placed at different distances in different directions from the pumping well or the packer test hole. This is a very time consuming and expensive exercise. On the other hand, the discrete fracture flow approach [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17 and 18] along with a good set of fracture data and results available from a few hydraulic tests performed at the field site has the capability to investigate and capture the hydraulic behaviour of the fractured rock mass at any scale, irrespective of the distance between the pumping well and the observation wells or between the packer test hole and the observation wells, without worrying about the REV size of the fractured rock mass.
This paper shows how the following objectives were accomplished for a gneissic rock mass using fracture data obtained from several scanlines and a borehole, and single hole packer test results obtained from the same borehole: (a) to build and calibrate a stochasticdeterministic fracture network model; (b) to numerically simulate the conducted packer tests in the borehole using a discrete fracture fluid flow approach and to estimate the influenced region for the packer tests; (c) to numerically simulate the conducted packer tests in the borehole using a discrete fracture fluid flow approach to calibrate the flow rate per unit hydraulic gradient of fractures in the discrete fracture fluid flow model; (d) to investigate whether an REV for hydraulic properties exists for the rock mass around the borehole and if so, to estimate the REV size; (e) If an REV for hydraulic properties exists for the rock mass, to determine the equivalent continuum hydraulic conductivity tensor in 3D for the rock mass around the borehole which is valid for block sizes equal to or greater than REV size.
2. Fracture system modeling in 3D including a validation
The gneissic metamorphic rock mass dealt within the paper is composed of lightcoloured quartz and feldsparrich layers interfoliated with darker maficrich layers. It is located in the Arrowhead East Tunnel site in San Bernardino, USA. The Arrowhead East Tunnel construction is a part of the Inland Feeder Project currently in progress and managed by the Metropolitan Water District of Southern California [19]. The Inland Feeder Project is designed to more than double the Metropolitan Water District's water delivery capacity and to provide Southern California with up to 650 million gallons of additional water a day.
At the site a borehole of length 396 m (borehole C27 in Fig. 1) provided 831 fracture data through acoustic televiewer technique. Packer test data were available from various depths from the same borehole. For each packer interval, flow tests have been performed using several injection pressures. Table 1 provides the water head equivalent to one of the injection pressures used to conduct the packer test at each interval along with the corresponding flow rate. Sixteen scanlines (LS 7 through LS 22 in Fig. 1) that were within 750 m from the aforementioned borehole provided 859 more fracture data. Thirteen horizontal and three vertical scanlines drawn on subvertical exposures (dip angle between 70° and 90°) were used to conduct fracture mapping to obtain the fracture data. It was decided to build a 3D discrete fracture fluid flow model to simulate packer tests conducted at various depth intervals in the borehole. The packer test results should depend on the fractures that intersect the borehole, the fracture system around the borehole and the hydraulic boundary conditions. Even though all the important fractures that intersected the borehole were known through acoustic televiewer information, only a few samples of fractures were available around the borehole through scanline surveys conducted. Therefore, it was decided to build a stochastic fracture network to represent the fracture system around the borehole based on the scanline fracture data and a deterministic fracture system based on acoustic televiewer data to represent the fractures that intersected the borehole.
Fig. 1. Relative locations of the borehole C27 and scanlines LS 7–LS 22 at the rock mass site.
Table 1. Information about packer tests
Recently, Kulatilake [20] completed a software package named FRACNTWK based on information given in more than 10 journal papers that he published on the topic of fracture characterization and network modeling in the rock mechanics literature between 1984 and 1997. This package can be used to analyze discontinuity data obtained from boreholes, scanlines and 2D exposures such as rock outcrops, tunnel walls, tunnel roofs, etc., to perform fracture characterization and network modeling for discontinuous rock masses and to generate rock discontinuity systems in rock masses. This computer package provides procedures: (a) to identify statistically homogeneous regions in a rock mass [21, 22 and 23]; (b) to identify discontinuity clusters in a statistically homogeneous region [23 and 24]; (c) to apply corrections for sampling biases associated with orientation, spacing and trace length distributions of discontinuity clusters [23, 24, 25 and 26]; (d) to obtain best probability distributions for orientation, spacing, trace length and discontinuity size in 3D of discontinuity clusters [23, 24 and 27]; (e) to obtain a map of the discontinuity traces sampled through either scanline or area sampling surveys [20]; (f) to estimate 1D discontinuity frequency along mean normal vector directions of discontinuity clusters using discontinuity spacing data mapped from some other directions [23 and 24]; (g) to estimate 1D discontinuity frequency in any direction in the rock mass [20]; (h) to estimate distributions for block size, number of blocks per unit volume and number of discontinuities per unit volume for the rock mass [20, 23 and 24]; (i) to estimate fracture tensor parameters for each discontinuity cluster as well as for the rock mass [20 and 28]; (j) to generate discontinuities in 3D for the rock mass and to obtain discontinuity trace predictions on vertical and horizontal exposures [20 and 24]; and (k) to verify the used discontinuity system models [20 and 24]. The package has 26 calculation and 24 graphical programs. Some programs from this package were used to build the stochastic fracture network model around the borehole. To cut down the space, only a summarized account of the built fracture network is reported in this paper.
The scanline data indicated that the fracture data are coming from a statistically homogeneous region. Fig. 2 shows the orientation distribution of poles of fracture sets on an upper hemispherical polar equal area projection. Table 2 shows the summary of fracture set delineation results obtained for the gneissic rock mass based on the scanline data. Both the Bingham and hemispherical normal distributions were found to be unsuitable to represent the statistical distribution of orientation of fracture sets [29]. Therefore, empirical distributions obtained for orientation of fracture sets were used to generate orientation of fractures in 3D surrounding the packer test regions. Performed fracture trace length, 3D size and 3D intensity modeling resulted in Table 3. The gamma distribution was found to be the best distribution to represent fracture trace length and size in 3D.
Fig. 2. Orientation distribution of poles of fracture sets on an upper hemispherical polar equal area projection.
Table 2. Summary of fracture set delineation results for the gneissic rock mass
Table 3. Summary of 3D fracture size and intensity statistics
For a statistically homogeneous rock mass, to describe the fracture geometry pattern in 3D, it is necessary to specify the number of fracture sets, and the statistical distributions for the following fracture geometry parameters for each fracture set: (1) number of fractures per unit volume; (2) orientation; (3) diameter; and (4) location of fracture centres. Because the fracture spacing for each fracture set followed an exponential distribution [29], according to the statistical theory, the Poisson distribution can be used to model the 3D fracture intensity distribution for each fracture set with the calculated mean 3D fracture intensity value. The empirical distribution obtained for the orientation was used to model the orientation distribution for each fracture set. For each fracture set, the diameter was represented by a gamma distribution. Since the fracture spacing of each fracture set followed an exponential distribution, according to the statistical theory, the location of fracture centres in 3D can be modeled using a uniform distribution. These statistical models were used to generate the fracture system in 3D for the gneissic rock mass around the borehole.
The generated fracture system in 3D was used to make predictions of fracture traces for each fracture set on a 2D window. The size of the 2D window for each fracture set was selected in such a way to make a comparison with the fractures mapped using the scanline survey on the 2D exposure that had highly irregular boundaries. To validate the generated discontinuity system, the predicted fracture traces on selected 2D windows from fracture sets were compared with actual field line surveys. Three criteria were used to decide whether the predicted data are in close agreement with the field data. A close agreement of the apparent orientation distribution of the chosen fracture set between the fractures mapped on the exposure in the field and the traces predicted on the selected 2D window was considered as the first criterion. As the second criterion, the 1D intensity value along a scanline placed around the midlevel of the constructed 2D window for the fracture set was expected to be in good agreement with the 1D intensity value obtained along the scanline in the field survey. The third criterion used was the predicted mean trace length of the fracture set on the 2D window should be in agreement compared to the range between the observed and corrected mean trace lengths obtained from the size analysis performed on field fracture data for the same fracture set. Performed comparisons indicated that for all fourfracture sets, the apparent orientations are in good agreement [29]. Figs. 3a and c show the fractures from regional cluster nos. 2 and 4, respectively, that intersected a horizontal scanline that was drawn at the midlevel on a vertical 2D window that was placed in a 3D volume used to generate fractures in 3D. Figs. 3b and d show the fractures from the same two sets, respectively, obtained from field mapping on the same scanline. Figs. 3a–d compare quite well with respect to statistical properties of fracture traces. For some fracture sets, the diameter and the 3D fracture intensity parameter values that were used to generate fractures had to be changed slightly until the generated fractures satisfied the aforementioned three criteria quite well. Table 4 provides the modified values of the summary statistics for the diameter and 3D intensity parameters used for the generation. These parameter values are termed as “after calibration values”. Table 3 and Table 4 can be compared to see the changes made to the fracture network model to produce generated fracture traces that compare well with the fracture traces obtained from field measurements.
Fig. 3. Comparison between predicted and field observed fracture traces on scanline LS 14: (a) predicted fracture traces of regional cluster 2 on the scanline drawn on a vertical 2D window placed in the 3D volume used to generate fractures for the gneissic rock mass around the borehole, (b) fracture traces of regional cluster 2 obtained from field mapping on the scanline, (c) predicted fracture traces of regional cluster 4 on the scanline drawn on a vertical 2D window placed in the 3D volume used to generate fractures for the gneissic rock mass around the borehole, (d) fracture traces of regional cluster 4 obtained from field mapping on the scanline.
Table 4. Summary of 3D fracture size and intensity statistics of fracture sets used for generation—after calibration
3. Features of the 3D discrete fracture fluid flow approach
3.1. General features of conceptual pipe flow network model
Tsang and Tsang [30] suggested that fluid flow in a fracture occur through one or many channels. Thus, it is reasonable to use equivalent onedimensional pipes to represent conceptually the fluid flow in a 2D fracture in order to simplify and minimize the complexity of fluid flow in a fracture [5, 11, 13, 14, 15, 16 and 31]. In fact, it is very difficult or impossible to represent the actual flow situation in a fracture that exist in the subsurface because it depends very much on the highly complex, stress dependant, spatially variable aperture distribution of the fracture. One of the most important aspects in a fluid flow model for a fractured rock mass is to represent correctly, at least relatively with respect to the magnitude, the hydraulic properties of different fracture sets having different fracture geometry properties.
Fractures were generated as equivalent circular discs in 3D. Intersection lines between these fractures were determined (Fig. 4). Between any two intersected fractures, two equivalent single flow pipes as shown in Fig. 4 represent flow from one fracture to another fracture. A number of possibilities exist to place these pipes on the fractures. In the developed model, the pipes are placed between the centre of the fracture disc and the midpoint of the intersected line of the intersected fractures. Through this way a network of pipe flow was formed to represent the fluid flow between intersected fractures ( Fig. 4). Conceptually, the transmissivity multiplied by width of flow path is used as the hydraulic parameter to represent the hydraulic property of a linear flow pipe under the steadystate flow. The intact rock permeability was ignored because it was negligible compared to the permeability of the fractures.
Fig. 4. Conceptual pipe flow model in general used to build the discrete fracture fluid flow model.
A fracture that is not intersected at the borehole is termed as background fracture in this paper. Fig. 5a shows a particular fracture F_{c} that is in a flow path and also belongs to the group of background fractures. This fracture is intersected by fractures F_{b}, F_{d}, F_{e} and F_{f} that are also background fractures. Note that each background fracture that contributes to flow should intersect either at least two other fractures or at least one fracture and a boundary. The flow into the background fracture should occur through at least one intersection and flow out of the fracture should occur at least through another intersection. The region of flow from the centre of a background fracture disc to an intersection interface or vice versa can be considered as a trapezoid shown in Fig. 5b, where b_{C} is the width at one end of the flow path crossing the diameter of the fracture disc and b_{L} is the width at the other end of the flow path equal to the length of the intersection line between the background fracture and another fracture or boundary. The width b_{C} may be approximated using the following equation:
Fig. 5. Assumed features for background fractures: (a) intersections of a background fracture with other background fractures, (b) change in width of a flow subdomain of a background fracture and, (c) a flow subdomain of a background fracture represented by an equivalent mean width of flow.
Application of Darcian law to Fig. 5b gives
Substituting for b in Eq. (2) by Eq. (3), the following equation is obtained:
In Fig. 5c a flow region with equivalent mean width b_{m} and mean aperture m is assumed to produce the same flow rate under the same hydraulic boundary conditions between the two ends of the flow region. By applying the Darcy's law to Fig. 5c, the following equation is obtained:
Similarly, the following equation is suggested to calculate width b_{C} for a fracture that intersects a borehole in an interval of packer test such as fracture F shown in Fig. 6:
Fig. 6. Assumed features for fractures that intersect the borehole: (a) intersections between a fracture that intersects the borehole and background fractures, (b) change in width of flow subdomains of a fracture that intersects the borehole at the immediate vicinity of the borehole.
3.2. Specific features of conceptual pipe flow network model
A procedure and a computer program called FRACPIPES [32] was developed as a preprocessor to determine both an average width of flow paths for fractures intersected with a packer test interval in a borehole (b_{m})_{bh}, and an average width of flow paths for background fractures (b_{m})_{bg}, for each fracture set. It was interesting to note that only small differences (less than 15%) were found between the calculated values of (b_{m})_{bg} or (b_{m})_{bh} for different fracture sets even though the size distributions were different for those four fracture sets.
In the 3D conceptual linear pipe discrete fracture fluid flow model, the average value of the transmissivity for each fracture in a fracture set needs calibration. For the total borehole depth in which the packer tests were conducted, the fracture data obtained through acoustic televiewer technique exhibited a statistically homogeneous region with depth [29]. The in situ observation data for the fractures exposed on the borehole walls indicated that the average values of the aperture and the extent of the filling in fractures for different fracture sets with different orientations to be almost the same. It provided the support on the point that the difference between the average values of the transmissivity for different fracture sets is not significant. In addition the analysis of the packer test results for different depths in the borehole showed that no strong correlation exists between the hydraulic conductivity obtained from the equivalent isotropic continuum porous media approach and the fracture composition from different fracture sets intersecting the borehole ( Fig. 7). Thus, to simplify the model, it is reasonable to assume that the average value of the transmissivity for a fracture is the same for all fourfracture sets for a specific depth range of the packer tests. Note that for the investigated site, no data were available to incorporate the effect of stress resulting from depth on transmissivity of a fracture. Therefore an average value was used for the transmissivity of a fracture to represent a depth range of 6.5 m.
Fig. 7. Percentage of fractures belonging to each fracture set and hydraulic conductivity per fracture (calculated based on the equivalent continuum porous media approach) intersected with the packer test borehole with depth.
If the average transmissivity of a fracture that intersects the borehole is assumed to be the same as the average transmissivity of a background fracture, then the ratio between the average flow rate per unit hydraulic gradient (defined as the transmissivity multiplied by the mean width of flow paths) of a fracture that intersects the borehole (HC_{m})_{bh}, and the average flow rate per unit hydraulic gradient of a background fracture (HC_{m})_{bg}, will be the same as the ratio between (b_{m})_{bh} and (b_{m})_{bg}. Based on the aforementioned features of the fracture network model and the conceptual pipe flow model, a conceptual equivalent pipe flow network model was developed to build the discrete fracture fluid flow model (Fig. 8). In the conceptual equivalent pipe flow network model, an equivalent pipe is specified for the flow path between the centre of a fracture disc and the midpoint of the intersection between the fracture and each of the other adjacent fractures. Further, an average value of the flow rate per unit hydraulic gradient resulting from all the fractures in the background is assigned to each pipe that connects between background fractures and an average value of the flow rate per unit hydraulic gradient resulting from all the fractures that intersect the borehole is specified to each pipe that connects between the borehole and the background. In the 3D conceptual linear pipe discrete fracture fluid flow model, those two average values of the flow rate per unit hydraulic gradient should be calibrated. Section 4 deals with such calibrations.
Fig. 8. Conceptual equivalent pipe flow model used to build the discrete fracture fluid flow model in 3D.
4. Estimation of influenced region for packer tests and calibration of flow rate per unit hydraulic gradient
To numerically simulate each packer test, only two observed data were available: (a) the injection water pressure of the packer test and (b) the injection flow rate of the packer test. It was decided to use the injection water pressure value as the input to the numerical discrete fracture fluid flow model and predict the injection flow rate as the output. Only the injection flow rate of the packer test was available to compare with the predicted flow rate to calibrate the flow rate per unit hydraulic gradient field in the considered discrete fracture fluid flow model. On the other hand, at the minimum level (HC_{m})_{bh} and (HC_{m})_{bg} should be calibrated from each packer test results. It should be pointed out that the calculation of the ratio of (HC_{m})_{bh} to (HC_{m})_{bg} was dealt with in Section 3. The developed concepts about the flow rate per unit hydraulic gradient in Section 3 were used along with the packer test results in calibrating the average or representative flow rate per unit hydraulic gradient field in the numerical discrete fracture fluid flow model.
The validated stochastic fracture network was used to generate fractures in a cube of size 30.5 m. Cubic samples of sizes 6.6, 9.1, 12.2, 15.2, 18.3 m were selected from the centre of the cube of size 30.5 m to investigate the effect of sample size on calibration of flow rate per unit hydraulic gradient for fractures. It was necessary to conduct this investigation to determine the influenced region or effective region of the packer tests. In each of these samples, the borehole was placed vertically at the centre of the cube to have the centre of the packer interval to coincide with the centre of the cube (Fig. 9). Note that the vertical boundaries of the cube are denoted as B1, B2, B3 and B4, and the two horizontal boundaries are denoted as B5 (top) and B6 (bottom). The fractures from the stochastic fracture network scheme intersecting the borehole were replaced by the observed fractures on the borehole for the considered length. For fractures that intersected the borehole, the known dip, dip direction and the location of each fracture on the borehole were used in placing the observed fractures. However, no information was available about the size of these fractures. To obtain an estimate for this size, several imaginary vertical boreholes of the same size as the actual borehole were placed at a number of random locations in the 30.5 m cube that had the stochastic fracture system. Then the size of the diameter of each intersected fracture and the fracture set it belongs to were recorded. This information was used to estimate the mean value for the diameter of each fracture set that appeared on the borehole. These mean diameter values were used in representing the fracture size of the observed fractures on the borehole wall.
Fig. 9. Cubic sample box with the packer interval used to numerically simulate the packer test.
The fractured rock mass was considered as an unconfined aquifer in 3D. The matrix rock was considered as impermeable. To implement the conceptual and numerical model, the computer program FRACPIPES was used as a preprocessor to obtain all the necessary information including intersections between the fractures, intersections between the fractures and boundaries, (b_{m})_{bh} and (b_{m})_{bg} values, the coordinates (X,Y,Z) of the ends of pipes (nodes), the lengths of pipes and nodal identifications (constant head or constant flow rate and an internal node or boundary node) for each node. Then a modified version of the TRINET computer program [13] was used to calculate the flow through the fracture network.
The constant head boundary condition was applied to all the boundaries of each selected cube. The ratio between (HC_{m})_{bh} and (HC_{m})_{bg} (i.e. the same as that between (b_{m})_{bh} and (b_{m})_{bg}) was calculated for each selected cube. A (HC_{m})_{bg} value was selected and then using the calculated ratio between (HC_{m})_{bh} and (HC_{m})_{bg}, the corresponding (HC_{m})_{bh} value was calculated. The packer test was simulated under steadystate hydraulic condition using the injected water pressure as the head difference between the borehole and the constant head boundaries of the selected cube and using (HC_{m})_{bh} and (HC_{m})_{bg} as the flow rate per unit hydraulic gradient values for all the fractures that intersected the borehole and the background fractures, respectively. Fig. 10 shows the variation of predicted injected flow rates at different boundaries with the block size for the packer test conducted at the depth range 388.3–394.5 m using (HC_{m})_{bg}=3.3×10^{−8} m^{3}/s. The horizontal flow rate can be considered to be almost the same for block sizes greater than or equal to about 10 m. On the other hand, the vertical flow rate from boundaries 5 and 6 becomes almost constant for a block size greater than or equal to 15 m. Note that the length of the packer is about 6.6 m. This means that if the boundaries of the cube are placed more than about 4.2 m from each edge of the packer, then the flow rate from each boundary as well as the total injection flow rate stay almost the same with the block size. This block size of 15 m can be considered as the influenced region corresponding to the packer test simulation in the depth range for the gneissic rock mass. Similarly, the influenced region for other depth ranges turned out to be about 15 m.
Fig. 10. Variation of flow rates at different boundaries of the cubic block with block size for the depth range 388.3–394.5 m using (HC_{m})_{bg}=3.3×10^{−8} m^{3}/s.
The resulting injected flow rate was then compared with the observed injected flow rate. This procedure was repeated using different selected (HC_{m})_{bg} and corresponding calculated (HC_{m})_{bh} values until excellent agreement was obtained between the predicted and observed injected flow rates. These resulting (HC_{m})_{bg} and (HC_{m})_{bh} values are called calibrated values. Fig. 11 shows how the calibrated flow rate per unit hydraulic gradient for background fractures vary with the block size for the same depth range. This plot also confirms that the influenced region corresponding to the packer test simulation for this rock mass is around 15 m. To be on the conservative side and also to accommodate higher injection effective water pressures at other depth ranges, for this rock mass, the influenced block size value corresponding to packer test simulation was selected as 18 m.
Fig. 11. Variation of calibrated mean flow rate per unit hydraulic gradient for background fractures with block size for the depth range 388.3–394.5 m.
The calibrated flow rate per unit hydraulic gradient values for both background and borehole fractures resulting from a block size of 18 m and different packer tests conducted at different depth regions are shown in Fig. 12. These flow rate per unit hydraulic gradient values were obtained by matching the predicted injection flow rates with observed values almost perfectly.
Fig. 12. Calibrated mean flow rates per unit hydraulic gradient for background fractures and the fractures that intersect the borehole at different depth regions of packer tests.
5. Estimation of REV size
The validated stochastic fracture network was used to generate fractures in a cube of size 30.5 m. Cubic samples of sizes 6.6, 9.1, 12.2, 15.2, 18.3 m were selected from the centre of the cube of size 30.5 m to estimate the REV size for the gneissic rock mass around the borehole. Each cubic block had two vertical faces perpendicular to the North–South direction. The calibrated (HC_{m})_{bg} value for the depth range 388.3–394.5 m was assigned for the fractures in each block. To calculate the directional hydraulic conductivity in a chosen direction, a unit hydraulic gradient was applied in that direction using the two perpendicular planes of the cube corresponding to the chosen direction. The other boundaries of the block were specified as no flux boundaries in performing the flow calculations under the steady state using the finite element method. The directional hydraulic conductivities obtained in the directions N–S, E–W and vertical for each block are shown in Fig. 13. It is clear from the figure that the hydraulic conductivities in different directions do not change with the block size (in practical sense) for block sizes greater than about 12.5 m. The same was found to be applicable for other depth ranges. This implies that the REV size for the gneissic rock mass can be considered to be around 12.5 m. However, to be on the conservative side, the REV size to express the hydraulic behaviour of the gneissic rock mass was selected as 15 m.
Fig. 13. Directional hydraulic conductivity in horizontal and vertical directions vs. block size for the depth range 388.3–394.5 m.
6. Calculation of hydraulic conductivity tensor in 3D
6.1. General procedure
The following provides the procedure and equations to obtain the hydraulic conductivity tensor for a REV based on the calculated directional hydraulic conductivity along several directions of hydraulic gradient. The directional hydraulic conductivity K_{d}() is defined as follows:
The directional hydraulic conductivity along the gradient Δh can be expressed by
The following equation relates K_{d}() to incorporating the hydraulic conductivity tensor :
The following equation can be derived from Eq. (12):
If N values of the directional hydraulic conductivity along N specific gradient directions are available (N6), the system of N linear algebraic equations with six unknown variables K_{11}, K_{22}, K_{33}, K_{12}, K_{23}, and K_{13} can be represented in a matrix form by
6.2. Determination of 3D hydraulic conductivity tensor
To calculate the hydraulic conductivity tensor for each depth range a cube of size 15 m was chosen from the centre of the cube of size 30.5 m having the corresponding fracture network. For each depth range, this 15 m block was selected in a number of different orientations. The first block was selected in such a way to have two vertical boundaries of the cube to be perpendicular to the north direction. As mentioned earlier, to calculate the directional hydraulic conductivity in the chosen direction, a unit hydraulic gradient was applied in that direction using the two perpendicular planes of the cube corresponding to the chosen direction. The other boundaries of the block were specified as no flux boundaries in performing the flow calculations using the finite element method. The calibrated flow rate per unit hydraulic gradient value for the background fractures at the corresponding depth range was assigned for the fractures in the block. Similarly, by rotating the block at 30° increments around the vertical axis on the horizontal plane and applying the aforementioned procedure, the directional hydraulic conductivities were calculated at every 30° increment on the horizontal plane starting from the north direction. Then vertical planes were selected having strikes N–S, N30°E, N60°E, E–W, S60°E and S30°E. On each of these vertical planes, the directional hydraulic conductivities were calculated in the directions that made angles 15°, 30°, 60°, 120°, 150° and 165° with the vertical using the aforementioned procedure. Through this way, for each depth range, the directional hydraulic conductivity was calculated in more than 80 directions covering the 3D space. Fig. 14 shows the variation of directional hydraulic conductivity on a number of 2D planes for the depth range 304.5–311 m. These values were then used in the least square regression analysis procedure mentioned above in calculating the principal directions and principal values for the hydraulic conductivity tensor in 3D for each depth range. Fig. 15 shows the integration of the information given in Fig. 14 in 3D for the depth range 304.5–311 m. Both Fig. 14 and Fig. 15 show that the rock mass hydraulic conductivity is significantly anisotropic.
Fig. 14. Variation of directional hydraulic conductivity (in 10^{−8} m/s) on several 2D planes for the depth range 304.5–311 m: (a) on the vertical plane striking N–S, (b) on the vertical plane striking N30°E, (c) on the vertical plane striking N60°E, (d) on the vertical plane striking E–W, (e) on the vertical plane striking S60°E, (f) on the vertical plane striking S30°E and (g) on the horizontal plane.
Fig. 15. Variation of directional hydraulic conductivity in 3D space (scale: 1 cm=2.1×10^{−7} m/s) for the depth range 304.5−311 m: (a) view along the direction having TREND=320° and downward PLUNGE=20°, (b) view along the direction having TREND=20° and downward PLUNGE=20°.
Table 5 gives the principal directions for the hydraulic conductivity tensor for the gneissic rock mass around the borehole. Note that two of these directions are subhorizontal and the third one is subvertical. Also, note that three of the fracture sets in the rock mass are subvertical and the fourth fracture set is subhorizontal (Table 2). It is interesting to compare between the directions of fracture sets ( Table 2) and principal hydraulic conductivity directions ( Table 5). Fig. 16 shows how the magnitude of hydraulic conductivity tensor varies with depth around the borehole. Note that the hydraulic conductivity of the rock mass is significantly anisotropic. The highest hydraulic conductivity is obtained along a subvertical direction. This compares well with the fracture system of the rock mass that has three subvertical fracture sets. Note that the ratio between the major principal hydraulic conductivity (K_{1}) and minor principal hydraulic conductivity (K_{3}) ranges between 3.5 and 4 for the rock mass around the borehole. Fig. 16 also compares geometric mean hydraulic conductivity obtained through the principal hydraulic conductivity values calculated against the hydraulic conductivities obtained through the radial flow continuum porous media assumption for the rock mass around the borehole. These two hydraulic conductivities compare quite well.
Table 5. Ranges for directions of principal hydraulic conductivities obtained from different depth regions (note 0R3)
Fig. 16. Principal hydraulic conductivities for different depth regions or intervals of packer test.
7. Discussion
The equivalent linear pipe fluid flow network model has the capability to represent the characteristics of the fluid flow in a fracture system with the model simplicity by applying an equivalent hydraulic parameter such as the flow rate per unit gradient (L^{3}/t) for a fracture. At the same time, the model needs less computer resources and could solve field problems at large scales. Moreover, because of the linear distribution of the head in pipes with the constant hydraulic parameter for a pipe segment between two intersection nodes of fractures, it is easier to calibrate the fluid flow model. In addition, the parameters for a fracture fluid flow system represented by a fluid flow pipe network can be adjusted easily in implementing the fluid flow model. Note that for the investigated site, no data were available to incorporate the effect of stress resulting from depth on fluid flow rate per unit hydraulic gradient of each equivalent pipe. Therefore an average value was used for the fluid flow rate per unit hydraulic gradient to represent a depth range of 6.5 m.
If the ratio of transmissivities for different sets of fractures can be determined, it is apparent that the proposed methodology to estimate the hydraulic conductivity tensor using a single well packer test is still valid. The transmissivity could be evaluated through an assessment of apertures of fractures, filling status, and hydraulic tests for different sets of fractures. Moreover, a similar procedure to assess the hydraulic behaviour can be implemented in cases where (b_{m})_{bh} or (b_{m})_{bg} is different for each fracture set.
So far, for each depth region, the hydraulic conductivity tensor was calculated based on one realization of the stochastic fracture network around the borehole. To determine whether the particular realization has a significant influence on the hydraulic conductivity tensor calculated, for the depth range 388.3–394.5 m the hydraulic conductivity tensor was calculated based on a second realization of the stochastic fracture network. Table 6 provides the magnitudes and directions for the principal hydraulic conductivities resulting from the two realizations. It is clear that from Table 6, the calculated hydraulic conductivity tensor is not affected in practical sense by the particular realization used for the stochastic fracture network.
Table 6. (a) Comparison of principal hydraulic conductivities, and (b) comparison of directions of principal hydraulic conductivities for the depth range, 388.3–394.5 m between two realizations of the stochastic fracture network
The plots for the reciprocal of directional hydraulic conductivity have shown approximately a shape of ellipse on different vertical planes. The coefficient of determination for a performed linear regression using the above procedure and equations in determining the hydraulic conductivity tensor for each of the packer test interval was about 0.95. It indicates that the selected cube of the fractured media provided a necessary property of REV for hydraulic behaviour and the determined hydraulic conductivity tensor had an acceptable accuracy on representation of calculated directional hydraulic conductivities.
As demonstrated above, 3D hydraulic conductivity tensor could be estimated using a single well packer test without other observation points. In fact, addition of an observation point or points which provide hydraulic response information for the region beyond the packer test depth range could improve the estimation through calibrations to the observed data from both packer test well and the observation point or those points.
8. Conclusions
The procedures used for fracture characterization, generation and validation
showed the capability of validating the developed stochastic fracture network
model for the gneissic rock mass of the site dealt with. The developed discrete
fracture fluid flow model had the capability of simulating the packer tests
conducted in the rock mass to estimate the influenced region for the packer test
and to determine the REV size for the rock mass with respect to hydraulic
behaviour. The discrete fracture fluid flow model also exhibited the capability
of calibrating the flow rate per unit hydraulic gradient field for the fractures
in the rock mass through simulation of packer tests. The average flow rate per
unit hydraulic gradient of the fractures that intersected the borehole was
considered to be quite different to that of the fractures at a significant
distance from the borehole. A conservative estimate of 18 m was obtained
for the block size as the influenced region to simulate each packer test that
had a length of about 6.5 m. Conservatively, a block size of 15 m was
estimated as the REV size to represent the hydraulic behaviour of the rock mass.
The calculated directional hydraulic conductivities and the 3D hydraulic
conductivity tensor show a significant anisotropy of the hydraulic behaviour of
the gneissic rock mass. The rock mass has three subvertical joint sets and a
subhorizontal foliation set. The major principal hydraulic conductivity
direction was found to be in a subvertical direction. This agrees well with the
fracture system that exist in the gneissic rock mass. It was shown how to
estimate the 3D hydraulic conductivity tensor in a rock mass using single well
packer test results along with the fracture information for the site. The
geometric hydraulic conductivity calculated at different depths was found to be
comparable to the hydraulic conductivity estimated through the continuum porous
media radial flow assumption.
Acknowledgements
Financial support received for this project from the Metropolitan Water District of Southern California is very much appreciated.
References
1. J. Bear Dynamics of fluids in porous media, American Elsevier, New York (1972).
2. Marsily GDe. Flow, transport in fractured rocks: connectivity and scale effect. Proceedings of the International Association of Hydrogeologists. Hydrogeology of Rocks of Low Permeability. Tucson, memoirs, vol. XVII, 1985. p. 267–77.
3. P.A. Hsieh and S.P. Neuman , Field determination of the threedimensional hydraulic conductivity tensor of anisotropic media, 1, theory. Water Resour Res 21 11 (1985), pp. 1655–1665. View Record in Scopus  Cited By in Scopus
4. P.A. Hsieh, S.P. Neuman, G.K. Stiles and E.S. Simpson , Field determination of the threedimensional hydraulic conductivity tensor of anisotropic media, 2, methodology and application to fractured rocks. Water Resour Res 21 11 (1985), pp. 1667–1676. View Record in Scopus  Cited By in Scopus
5. P.H.S.W. Kulatilake and B.B. Panda , Effect of block size and joint geometry on jointed rock hydraulics and REV. ASCE J Eng Mech 126 8 (2000), pp. 850–858. Full Text via CrossRef  View Record in Scopus  Cited By in Scopus
6. F.W. Schwartz, W.L. Smith and A.S. Crowe , A stochastic analysis of microscopic dispersion in fractured media. Water Resour Res 19 5 (1983), pp. 1253–1265. View Record in Scopus  Cited By in Scopus
7. Robinson PC. Connectivity, flow and transport in network models of fractured Media. PhD dissertation, Oxford University, Oxford, 1984.
8. J.C.S. Long and P.A. Witherspoon , The relationship of the degree of interconnection to permeability in fracture networks. J Geophys Res 90 B4 (1985), pp. 3087–3097.
9. J. Andersson and B. Dverstorp , Conditional simulations of fluid flow in threedimensional networks of discrete fractures. Water Resour Res 23 10 (1987), pp. 1876–1886. View Record in Scopus  Cited By in Scopus
10. A. Rouleau and J.E. Gale , Stochastic discrete fracture simulation of ground water flow into an underground excavation in granite. Int J Rock Mech Min Sci Geomech Abstr 24 2 (1987), pp. 99–112. Abstract  Abstract + References  PDF (1295 K)  View Record in Scopus  Cited By in Scopus
11. M.C. Cacus, E. Ledoux, G.De. Marsily, B. Tille, A. Barbreau, E. Durand, B. Feuga and P. Peaudecerf , Modeling of fracture flow with a stochastic discrete fracture network: calibration and validation—1, the flow model. Water Resour Res 26 3 (1990), pp. 479–489.
12. E.A. Sudicky and R.G. McLaren , The Laplace Transform Galerkin technique for largescale simulation of mass transport in discretely fractured porous formations. Water Resour Res 28 2 (1992), pp. 499–514. Full Text via CrossRef  View Record in Scopus  Cited By in Scopus
13. Segan S, Karasaki K. TRINET: a flow and transport code for fracture networks—user's manual and tutorial, LBL34834. Lawrence Berkeley Laboratory, Berkeley, 1993.
14. B.B. Panda and P.H.S.W. Kulatilake , Influence of discontinuity geometry parameters and transmissivity on hydraulic behavior of discontinuous rock. ASCE J Eng Mech 125 1 (1999), pp. 41–50. Full Text via CrossRef  View Record in Scopus  Cited By in Scopus
15. B.B. Panda and P.H.S.W. Kulatilake , Relations between fracture tensor parameters and permeability tensor parameters for discontinuous rock. ASCE J Eng Mech 125 1 (1999), pp. 51–59. Full Text via CrossRef  View Record in Scopus  Cited By in Scopus
16. W.S. Dershowitz and C. Fidelibus , Derivation of equivalent pipe network analogues for threedimensional discrete fracture networks by the boundary element method. Water Resour Res 35 9 (1999), pp. 2685–2691. View Record in Scopus  Cited By in Scopus
17. K. Karasaki, B. Freifeld and A. Cohen , A multidisciplinary fractured rock characterization study at Raymond field site, Raymond, CA. J Hydrol 236 12 (2000), pp. 17–34. SummaryPlus  Full Text + Links  PDF (1178 K)  View Record in Scopus  Cited By in Scopus
18. S. Nakao, J. Najita and K. Karasaki , Hydraulic well testing inversion for modeling fluid flow in fractured rocks using simulated annealing: a case study at Raymond field site, California. J Appl Geophys 45 3 (2000), pp. 203–223. SummaryPlus  Full Text + Links  PDF (3249 K)  View Record in Scopus  Cited By in Scopus
19. J. Gallanes, V. Romero and T. Tellers , On target: the arrowhead east and west tunnels. ASCE Civil Eng Mag 66 (1996), pp. 50–53. View Record in Scopus  Cited By in Scopus
20. Kulatilake PHSW. Software manual for FRACNTWK—a computer package to model discontinuity geometry in rock masses. University of Arizona, Technical Report, Metropolitan Water District of Southern California, 1998, submitted for publication.
21. P.H.S.W. Kulatilake, D.N. Wathugala, M. Poulton and O. Stephansson , Analysis of structural homogeneity of rock masses. Int J Eng Geol 29 (1990), pp. 195–211. Abstract  Abstract + References  PDF (685 K)  View Record in Scopus  Cited By in Scopus
22. P.H.S.W. Kulatilake, R. Fiedler and B.B. Panda , Box fractal dimension as a measure of statistical homogeneity of jointed rock masses. Int J Eng Geol 48 34 (1997), pp. 217–230.
23. P.H.S.W. Kulatilake, J. Chen, J. Teng, X. Shufang and G. Pan , Discontinuity geometry characterization for the rock mass around a tunnel close to the permanent shiplock area of the Three Gorges dam site in China. Int J Rock Mech Min Sci 33 (1996), pp. 255–277. SummaryPlus  Full Text + Links  PDF (2351 K)  View Record in Scopus  Cited By in Scopus
24. P.H.S.W. Kulatilake, D.N. Wuthagala and O. Stephansson , Joint network modelling, including a validation to an area in Stripa Mine, Sweden. Int J Rock Mech Min Sci 30 5 (1993), pp. 503–526. Abstract  Abstract + References  PDF (1937 K)  View Record in Scopus  Cited By in Scopus
25. P.H.S.W. Kulatilake and T.H. Wu , Estimation of mean trace length of discontinuities. Rock Mech Rock Eng 17 (1984), pp. 215–232. Full Text via CrossRef  View Record in Scopus  Cited By in Scopus
26. D.N. Wathugala, P.H.S.W. Kulatilake, G.W. Wathugala and O. Stephansson , A general procedure to correct sampling bias on joint orientation using a vector approach. Comput Geotech 10 (1990), pp. 1–31. Abstract  Abstract + References  PDF (1111 K)  View Record in Scopus  Cited By in Scopus
27. Kulatilake PHSW, Wu TH. Relation between discontinuity size and trace length. Proceedings of the 27th US Symposium on Rock Mechanics, Tuscaloosa, AL, 1986. p. 130–3.
28. P.H.S.W. Kulatilake, S. Wang and O. Stephansson , Effect of finite size joints on deformability of jointed rock at the threedimensional level. Int J Rock Mech Min Sci 30 5 (1993), pp. 479–501. Abstract  Abstract + References  PDF (1407 K)  View Record in Scopus  Cited By in Scopus
29. Kulatilake PHSW, Um J, Wang M, Gao H, Chen J. Stochastic discontinuity characterization for the metamorphic rock mass of the arrowhead east tunnel site, Inland Feeder Project. University of Arizona, Technical Reports, vols. 1 and 2, Metropolitan Water District of Southern California, 1998, submitted for publication.
30. Y.W. Tsang and C.F. Tsang , Channel model of flow through fractured media. Water Resour Res 23 3 (1987), pp. 467–479. View Record in Scopus  Cited By in Scopus
31. M. Wang, P.H.S.W. Kulatilake, B.B. Panda and M.L. Rucker , Groundwater resources evaluation case study via discrete fracture fluid flow modeling. Int J Eng Geol 62 (2001), pp. 267–291. SummaryPlus  Full Text + Links  PDF (467 K)  View Record in Scopus  Cited By in Scopus
32. Wang M. Discrete fluid flow modeling, field applications in fractured rocks. PhD dissertation, The University of Arizona, Tucson, AZ, 2000.
33.
J.V. Beck and K.J. Arnold Parameter estimation in engineering
science, Wiley, New York (1977).
^{} Corresponding author. Tel.: +15206216064; fax: +15206218330; email: kulatila@u.arizona.edu
International
Journal of Rock Mechanics and Mining Sciences Volume 39, Issue 7, October 2002, Pages 887904 
Result list  previous < 23 of 56 > next 

