Go to ScienceDirect® Home Skip Main Navigation Links
 
Home
Browse
Search
My Settings
Alerts
Help
 Quick Search   Title, abstract, keywords   Author e.g.  j s smith
 Search tips (Opens new window)   Journal/book title   Volume   Issue   Page     Clear all fields    
International Journal of Rock Mechanics and Mining Sciences
Volume 39, Issue 7, October 2002, Pages 887-904
Result list |  previous  < 23 of 56 >  next 

SummaryPlus Full Text + Links PDF (527 K)  View thumbnail images | View full size images
Add to my quick links    Cited by    E-mail article    Save as citation alert    Export citation + link    Set up a citation RSS feed (Opens new window)   
View Record in Scopus
Cited By in Scopus (10)

doi:10.1016/S1365-1609(02)00067-9    How to Cite or Link Using DOI (Opens New Window)  
Copyright © 2002 Elsevier Science Ltd. All rights reserved.

Estimation of REV size and three-dimensional hydraulic conductivity tensor for a fractured rock mass through a single well packer test and discrete fracture fluid flow modeling

M. Wanga, P. H. S. W. KulatilakeCorresponding Author Contact Information, E-mail The Corresponding Author, a, J. Uma and J. Narvaizb
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 volumenext term (REV) size and three-dimensional (3-D) hydraulic conductivity tensor for a fractured rock mass. First, a 3-D 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 stochastic-deterministic 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 3-D 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; previous termRepresentative elementary volumenext term; Hydraulic conductivity tensor; Discrete fracture fluid flow modeling; Hydrogeology; Calibration; Field study; Numerical modeling



1. Introduction

Currently used methods for estimating three-dimensional (3-D) 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 second-rank, symmetric, positive-definite 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 stochastic-deterministic 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 3-D 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 3-D including a validation

The gneissic metamorphic rock mass dealt within the paper is composed of light-coloured quartz- and feldspar-rich layers inter-foliated with darker mafic-rich 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 C2-7 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 sub-vertical exposures (dip angle between 70° and 90°) were used to conduct fracture mapping to obtain the fracture data. It was decided to build a 3-D 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.


Full Size Image

Fig. 1. Relative locations of the borehole C2-7 and scanlines LS 7–LS 22 at the rock mass site.

Table 1. Information about packer tests

Image

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 2-D 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 3-D 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 1-D 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 1-D 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 3-D 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 3-D surrounding the packer test regions. Performed fracture trace length, 3-D size and 3-D intensity modeling resulted in Table 3. The gamma distribution was found to be the best distribution to represent fracture trace length and size in 3-D.


Full Size Image

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 Image

Table 3. Summary of 3-D fracture size and intensity statistics

Image

For a statistically homogeneous rock mass, to describe the fracture geometry pattern in 3-D, 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 3-D fracture intensity distribution for each fracture set with the calculated mean 3-D 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 3-D can be modeled using a uniform distribution. These statistical models were used to generate the fracture system in 3-D for the gneissic rock mass around the borehole.

The generated fracture system in 3-D was used to make predictions of fracture traces for each fracture set on a 2-D window. The size of the 2-D 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 2-D exposure that had highly irregular boundaries. To validate the generated discontinuity system, the predicted fracture traces on selected 2-D 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 2-D window was considered as the first criterion. As the second criterion, the 1-D intensity value along a scanline placed around the mid-level of the constructed 2-D window for the fracture set was expected to be in good agreement with the 1-D 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 2-D 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 four-fracture 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 mid-level on a vertical 2-D window that was placed in a 3-D volume used to generate fractures in 3-D. 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 3-D 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 3-D 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.


Full Size Image

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 2-D window placed in the 3-D 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 2-D window placed in the 3-D 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 3-D fracture size and intensity statistics of fracture sets used for generation—after calibration

Image

3. Features of the 3-D 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 one-dimensional pipes to represent conceptually the fluid flow in a 2-D 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 3-D. 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 steady-state flow. The intact rock permeability was ignored because it was negligible compared to the permeability of the fractures.


Full Size Image

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 Fc that is in a flow path and also belongs to the group of background fractures. This fracture is intersected by fractures Fb, Fd, Fe and Ff 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 bC is the width at one end of the flow path crossing the diameter of the fracture disc and bL 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 bC may be approximated using the following equation:


bC=(2D/NL)×(bL/bTm), (1)

where D is the diameter of the background fracture, NL is the total number of intersections between the fracture and other fractures or boundaries, and bTm is the average of lengths of all the intersections between the fracture and other fractures or boundaries. Note that by considering both in-flow to the background fracture and outflow from the fracture, (2D)/NL expresses the mean segment length at the diameter of the background fracture for each sub-flow domain between the fracture and each of the intersections. Also note that NLgreater-or-equal, slanted2. In Eq. (1), bL/bTm can be considered as a weight factor with the assumption that bC is proportional to bL.


Full Size Image

Fig. 5. Assumed features for background fractures: (a) intersections of a background fracture with other background fractures, (b) change in width of a flow sub-domain of a background fracture and, (c) a flow sub-domain of a background fracture represented by an equivalent mean width of flow.

Application of Darcian law to Fig. 5b gives


Q=K(−dh/dr)mb, (2)

where Q is the flow rate (L3/t) at a given cross-section of a fracture flow path, K the hydraulic conductivity (L/t), −dh/dr the hydraulic gradient (dimensionless), m the aperture of fracture (L), and b the width of cross-section of the flow region at a distance r from the centre of fracture disk to the specific cross-section. The width b can be written as

b=bC+(bLbC)r/R, (3)

where R is the length between the two ends of the flow region.

Substituting for b in Eq. (2) by Eq. (3), the following equation is obtained:


Q=K(−dh/dr)m (bC+(bLbC)r/R). (4)

By integrating Eq. (4) (r from 0 to R and h from HC to HL), the following equation is obtained:

Q=Km((HCHL)/R)×(bLbC)/ln(bL/bC), (5)

where HC is the water head at the end of the flow region with width bC and HL the water head at the end of the flow region with width bL.

In Fig. 5c a flow region with equivalent mean width bm 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:


Q=Km((HCHL)/R)bm. (6)

By comparing (5) and (6), bm can be written as a function of bC and bL as follows:

bm=(bLbC)/ln(bL/bC) (7)

or

bm=(bCbL)/ln(bC/bL). (8)

From the generated fracture network, all bC and bL values can be calculated for the background fractures. Then from (7) and (8) an average bm can be determined for a specific set of background fractures. Let this bm value be (bm)bg. Similarly, for a fluid flow domain which is not exactly trapezoid, the same (7) and (8) can be used to approximate the average (bm)bg.

Similarly, the following equation is suggested to calculate width bC for a fracture that intersects a borehole in an interval of packer test such as fracture F shown in Fig. 6:


bC=(πD/NB)(bL/bTm), (9)

where D is the diameter of the borehole intersected by fracture F (Fig. 6), NB is the total number of intersections between the fracture F and background fractures, and bTm is the average of the lengths for all the intersections in the fracture F. Note that (πD)/NB expresses the mean segment length of the borehole perimeter for each sub-flow domain between a fracture intersection and the borehole. In Eq. (9), bL/bTm can be considered as a weight factor with the assumption that bC is proportional to bL. As for each background fracture, for each of the fractures that intersects the borehole in a packer test region, an equivalent mean width of flow bm can be expressed through (7) and (8), where bC is obtained through Eq. (9). From the generated fracture network, all bC and bL values can be calculated for the fractures that intersect the borehole. Then from (7) and (8) an average bm can be determined for a specific set of fractures that intersect the borehole. Let this bm value be (bm)bh. It is expected that (bm)bg to be greater than (bm)bh for a borehole with a small diameter because the average widths of flow sub-domains at the borehole are small compared to those resulting from background fractures.


Full Size Image

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 sub-domains 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 (bm)bh, and an average width of flow paths for background fractures (bm)bg, for each fracture set. It was interesting to note that only small differences (less than 15%) were found between the calculated values of (bm)bg or (bm)bh for different fracture sets even though the size distributions were different for those four fracture sets.

In the 3-D 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 four-fracture 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.


Full Size Image

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 (HCm)bh, and the average flow rate per unit hydraulic gradient of a background fracture (HCm)bg, will be the same as the ratio between (bm)bh and (bm)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 3-D 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.


Full Size Image

Fig. 8. Conceptual equivalent pipe flow model used to build the discrete fracture fluid flow model in 3-D.

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 (HCm)bh and (HCm)bg should be calibrated from each packer test results. It should be pointed out that the calculation of the ratio of (HCm)bh to (HCm)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.


Full Size Image

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 3-D. 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, (bm)bh and (bm)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 (HCm)bh and (HCm)bg (i.e. the same as that between (bm)bh and (bm)bg) was calculated for each selected cube. A (HCm)bg value was selected and then using the calculated ratio between (HCm)bh and (HCm)bg, the corresponding (HCm)bh value was calculated. The packer test was simulated under steady-state 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 (HCm)bh and (HCm)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 (HCm)bg=3.3×10−8 m3/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.


Full Size Image

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 (HCm)bg=3.3×10−8 m3/s.

The resulting injected flow rate was then compared with the observed injected flow rate. This procedure was repeated using different selected (HCm)bg and corresponding calculated (HCm)bh values until excellent agreement was obtained between the predicted and observed injected flow rates. These resulting (HCm)bg and (HCm)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.


Full Size Image

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.


Full Size Image

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 (HCm)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.


Full Size Image

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 3-D

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 Kd(Image) is defined as follows:


Image (10)

where Image is a direction vector, Image is vector flux, and Δh is the gradient in direction Image.

The directional hydraulic conductivity along the gradient Δh can be expressed by


Image (11)

where Image is the magnitude of Δh.

The following equation relates Kd(Image) to Image incorporating the hydraulic conductivity tensor Image:


Image (12)

This is a canonical ellipsoid with semi-axes K1(−1/2), K2(−1/2) and K3(−1/2), where K1, K2 and K3 are the three principal hydraulic conductivities.

The following equation can be derived from Eq. (12):


Image (13)

where K11, K22, K33, K12, K23, and K13 are six independent components of the hydraulic conductivity tensor, Image is the ith directional hydraulic conductivity along the hydraulic gradient unit vector Image.

If N values of the directional hydraulic conductivity along N specific gradient directions are available (Nmuch greater-than6), the system of N linear algebraic equations with six unknown variables K11, K22, K33, K12, K23, and K13 can be represented in a matrix form by


AX=b, (14)

where

Image

By applying a linear regression model, the estimator Ke=(Ke11,Ke22,Ke33,Ke12,Ke23,Ke13)T of the regression parameters X=(K11,K22,K33,K12,K23,K13)T obtained from the method of least squares has the form [33]

Ke=(ATA)−1ATb. (15)

A computer program was developed to compute the hydraulic conductivity tensor along with the coefficient of determination using the above procedure and equations [32]. In addition, the three principal values of the hydraulic conductivity tensor and the corresponding three principal directions were determined.

6.2. Determination of 3-D 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 3-D space. Fig. 14 shows the variation of directional hydraulic conductivity on a number of 2-D 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 3-D for each depth range. Fig. 15 shows the integration of the information given in Fig. 14 in 3-D for the depth range 304.5–311 m. Both Fig. 14 and Fig. 15 show that the rock mass hydraulic conductivity is significantly anisotropic.


Full Size Image
Full Size Image

Fig. 14. Variation of directional hydraulic conductivity (in 10−8 m/s) on several 2-D 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.


Full Size Image

Fig. 15. Variation of directional hydraulic conductivity in 3-D 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 sub-horizontal and the third one is sub-vertical. Also, note that three of the fracture sets in the rock mass are sub-vertical and the fourth fracture set is sub-horizontal (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 sub-vertical direction. This compares well with the fracture system of the rock mass that has three sub-vertical fracture sets. Note that the ratio between the major principal hydraulic conductivity (K1) and minor principal hydraulic conductivity (K3) 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 0less-than-or-equals, slantRless-than-or-equals, slant3)

Image


Full Size Image

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 (L3/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 (bm)bh or (bm)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

Image

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, 3-D 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 3-D hydraulic conductivity tensor show a significant anisotropy of the hydraulic behaviour of the gneissic rock mass. The rock mass has three sub-vertical joint sets and a sub-horizontal foliation set. The major principal hydraulic conductivity direction was found to be in a sub-vertical direction. This agrees well with the fracture system that exist in the gneissic rock mass. It was shown how to estimate the 3-D 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 three-dimensional 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 three-dimensional 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 three-dimensional 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 large-scale 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, LBL-34834. 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 three-dimensional 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 1-2 (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 3-4 (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 three-dimensional 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 Contact Information Corresponding author. Tel.: +1-520-621-6064; fax: +1-520-621-8330; email: kulatila@u.arizona.edu



International Journal of Rock Mechanics and Mining Sciences
Volume 39, Issue 7, October 2002, Pages 887-904
Result list |  previous  < 23 of 56 >  next 
 
Home
Browse
Search
My Settings
Alerts
Help
Elsevier.com (Opens new window)
About ScienceDirect  |  Contact Us  |  Terms & Conditions  |  Privacy Policy
Copyright © 2007 Elsevier B.V. All rights reserved. ScienceDirect® is a registered trademark of Elsevier B.V.