FRACTURE-DISTRIBUTION MODELING IN ROCK MASS USING

BOREHOLE DATA AND GEOSTATISTICAL SIMULATION

Katsuaki Koike 1 , Kazuya Komorida 2 , and Yuichi Ichikawa 3

1: Department of Civil Engineering, Faculty of Engineering, Kumamoto University, 2-39-1, Kurokami,

Kumamoto 860-8555, Japan. e-mail: [email protected]
2: Graduate School of Science and Technology, Kumamoto University, 2-39-1, Kurokami, Kumamoto 860-8555, Japan.
3: PASCO Co. Ltd., 1-1-2, Higashiyama, Meguro-ku, Toyko 153-0043, Japan. e-mail: [email protected]

ABSTRACT- Understanding spatial distribution of rock fractures is significant to various fields in

geosciences. It is, however, difficult to obtain a realistic fracture model because the amount of fracture data
is small and their ol cations are strongly biased in a study area. Data are usually taken by borehole
investigations and then, spatial correlation structures of fractures among boreholes with different directions
are needed to be clarified. For this problem, we focus on the fracture density along a borehole (number of
fractures per 1-meter interval) and appearance relation of azimuths (strikes and dips) between a fracture pair.
Semivariogram of the density, 1 (h), and cross-semivariograms of the two indicators for two d irections, 2 (h),
are produced. Firstly, fracture density map is produced using the 1 (h) and a sequential gaussian simulation.
Next, a direction of each simulated fracture is assigned using the data set, 2 (h), and ordinary kriging
combined with the principal component analysis. Finally, length of each fracture is defined considering the
distance and difference of azimuths between a fracture pair located closely. This proposed method was
applied to 629 investigation data by several boreholes in a granitic site situated in northeastern Japan. The
size of study area was defined as 6060 m2 . Horizontal distribution of fractures and continuities of them
along the vertical direction were successfully estimated. The most noteworthy features are found in that
high-density zones are located in no data areas and continuous fractures striking along the boreholes are
inferred. In addition, permeability of rock mass was calculated to be about 10-17 m2 from this distribution
model and a permeability tensor analysis.

Key Words : Fracture, Granite, Semivariogram, Stochastic simulation, Permeability tensor

INTRODUCTION in fractured reservoirs (e.g., National Research

Council, 1996; Coward et al., 1998; Adler and
Fractures with diverse origins and different Thovert, 1999). In special, characterization of
scales from microcracks to faults are generally hydraulic properties of rock masses has become
developed in rock masses. They were formed very important recently, because it is closely related
attributable to emplacement and cooling of rock to environmental problems.
bodies, crustal movements, and regional stress fields For this purpose, a general view of fracture
at different geologic stages. Understanding spatial distribution in a study area is required to be drawn
distribution of rock fractures is significant to various using stochastic methods. There are many papers
fields in geosciences, including hydrogeology for that have tried to clarify scaling laws of fractures
fracture-affected flow channels, and resource using fractal theory (e.g., Vignes-Adler et al., 1991;
exploration for vein-type mineral deposits and fluids Davy et al., 1992; Dawers et al. 1993; Watterson et
al., 1996; Koike and Kaneko, 1999), but most of
N
N
them deal with length distribution over wide scales.
KDH-1
In addition to the lengths, locations of fractures
should be considered to reveal spatial correlation
KDT-1
structures of fracture attributes. Geostatistical

10
techniques that can incorporate the structures have KM-2

K D A-8

K D A-
also been applied to fracture-distribution modeling
(Long and Billaux, 1987; Young, 1987; Chils,
1988; Gringarten, 1996). It is, however, difficult to
KDK-1
KDT-2
obtain a realistic model because the amount of
KDH-3
fracture data is small and their locations are strongly 15 m
biased in a study area. Fracture attributes that are Figure 1. Arrangement of boreholes and segments that
applicable to distribution modeling considering represents the locations and strikes of the fractures.
length, appearance pattern, and azimuth of fracture
should be firstly identified. directions, N15W and N75E, at almost the same
As the first step to a fracture-distribution level, 250 m above the sea. Orientations (strikes and
modeling over various scales, this paper studies the dips), in addition to kinds of interstitial minerals,
fracture data obtained by boreholes in granitic rocks. width of alteration halo, and fracture width were
We focus on the fracture density and appearance measured for 629 fractures, by the examination of
relation of azimuths between a fracture pair, and the borehole TV images and drilling cores. Figure 1
aims at inferring three-dimensional distribution of shows the arrangement of boreholes and segments
joint-sized fractures. that represent the locations and strikes of the
fractures.
VARIOGRAM ANALYSIS
Semivariograms
Borehole Fracture Data Clarifying whether the spatial distributions of
The fracture data used in this study were fracture attributes are perfectly random or have
obtained by eight horizontal boreholes in the certain spatial correlation structures is the most
Kamaishi mine, northeastern Japan. The mine is fundamental process for fracture -distribution
located in an Early Cretaceous dioritic granite modeling. By producing experimental semi-
(partly by diorite). The boreholes were dilled over variograms of several attributes, it was found that
215-m length in total and toward two perpendic ular the widths of alteration halos have almost a pure
7 1.6
Semivariogram (h)
Semivariogram (h)

6 1.5

5 1.4

4 1.3

3 1.2

2 1.1
0 1 2 0 1 2 3 4
Distance between data h (10 m) Distance between data h (10 m)

Figure 2. Semivariogram of fracture densities, number of Figure 3. Semivariograms of first principal values for the
fractures per 1-m interval along each borehole. indicator sets of four strike sectors.
 nugget effect, whereas the widths of interstitial N
minerals contain a weak correlation structure. As
compared to these attributes, fracture densities,
which are defined by the number of fractures per
1-m interval along each borehole, showed a more
clear feature. It can be demonstrated by the
experimental semivariogram in Figure 2. A Kamaishi mine
spherical model drawn by the broken line is
applicable to approximate the semivariogram.
Another attribute considered here is the spatial
relation of fracture orientations which means
directional similarity between a fracture pair. We did
not use directly the difference of strike angles or dip
angles between two fractures, but attached great
importance to the distribution relation of strike or
dip data classified into several groups. This idea was 10 km
proposed to consider continuit ies of the same Figure 4. Lineaments around the Kamaishi mine
fracture set and corporate conjugate-pattern features extracted from the Landsat TM band 4 image using the
into fracture-distribution modeling. For this purpose, Segment Tracing Algorithm (STA).
the four directional sectors for the EW, NW-SE, NS,
and NE-SW are defined, and each fracture strike is experimental semivariograms are needed for the
transformed into an indictor set (I1 , I2 , I3 , I4 ) that first to third principal values. Figure 3 shows the
represents presence (Ii =1) or absence (Ii =0) of strike experimental semivariogram of the first principal
in the (EW, NW-SE, NS, NE-SW). For example, values. While the values are scattered, the
N45W strike defines the set as (0, 1, 0, 0 ). semivariogram can be approximated by an
This set requires four semivariograms for the exponential model (the broken line in Fig. 3). The
same sector and six cross-semivariograms for the same model is applicable to the semivariograms of
different sector pairs. To reduce calculation amount, second and third principal values.
the principal component analysis was adopted and To validate the spatial correlation structures
the set data were transformed into four principal clarified in the two semivariograms, satellite image-
values. Since the fourth principal values become derived lineaments were investigated as another
constant for the above indicator set, only three fracture-related geologic element. Figure 4 depicts

15 1.6
Semivariogram (h)

Semivariogram (h)

14
1.4
13
12 1.2
11
A 1.0 B
10
90 0.8
2 4 6 8 10 0 1 2 3 4 5 6
Distance between data h (10 4 m)
Distance between data h (103 m)
Figure 5. Semivariograms of A, the densities expressing number of lineament centers in
unit mesh; and B, the first principal values for indicator sets of four strike sectors.
 0.28
0.28 data were simply transformed into indicator data, 1
for northward dip and 0 for southward dip. The
0.26
0.26 experimental semivariogram of these indicators is
Semivariogram (h)

shown in Figure 6, which can be approximated by a

0.24
0.24 spherical model.

0.22
0.22 TWO-DIMENSIONAL MODELING OF
FRACTURE DISTRIBUTION
0.20
0.2

Production of Fracture-Density Map

0.18
0.18 The fracture density was chosen as an
00 110 220 330 440 550 important key to sow the seeds of fractures which
Distance between data h (10 m) are bases on simulating fracture distribution. It is
Figure 6. Semivariogram of the indicator-transformed dip necessary to extend the one-dimensional fracture
direction data. Indicators are 1 for northern dip and 0 densities along the boreholes to two-dimensional
for southern dip. densities over the study area. The study area,
defined as 6060 m2 in size, was superimposed on a
the lineaments around the Kamaishi mine extracted grid with 11 m2 meshes. Calculation of fracture
from the Landsat TM band 4 image. The Segment density in each mesh is the first step for the seeds, to
Tracing Algorithm (STA: Koike et al., 1995) was which the semivariogram of fracture densities (Fig.
adopted to extract the lineaments. Using the centers 2) is available.
of lineaments, two experimental semivariograms of We used ordinary kriging (OK) for the
the densities expressing number of centers in unit calculation, but could not obtain an appropriate
mesh and the first principal values for the indicator result as shown in Figure 7A. It is obvious that the
sets of strike sectors were produced as shown in densities are strongly affected by the nearest data
Figure 5. Since they can be approximated by the and there are many linear, unnatural boundaries of
same semivariogram models as the case of borehole densities. To overcome this problem, a sequential
fractures, the semivariograms shown in Figures 2 gaussian simulation (SGS: Journel, 1989) was
and 3 are considered to allow to estimate fracture examined. It consists of transformation of the data
distribution in the study area. into normal scores, determination of random path,
As for the fracture dips, only the dip directions Monte-Carlo simulation using kriging variance at
indicated a spatial correlation structure. The dip each grid point to draw an estimation value, addition

A B
Number of
N fractures in 1 m2

15 m
Figure 7. Comparison of two fracture-density distributions produced by A, OK; and B, SGS.
of the value to the data set, and repetition of the measured fractures at each fracture center. The
calculation. The validness of SGS can be confirmed interpolated value above 0.5 assumes the dip
from the result (Fig. 7B) in that high-density zones direction of fracture to be northward, otherwise
are estimated in no data areas in the northwest and southward.
unnatural boundaries are not formed. In addition, (4) Strike and dip angles at each fracture center are
the density distribution shows large spatial variation stochastically drawn by combining the
as compared to that produced through OK. cumulative distribution functions of those angles
in the selected strike sector and dip direction with
Fracture-Distribution Modeling Monte-Carlo method.
A method of fracture-distribution modeling (5) Connection of fracture centers considering both
proposed here utilizes the fracture-density map and the center distance and directional difference is
consists of the following five steps. carried out at the final step. If an arbitrary
fracture center (A) can find out another center (B)
(1) In each mesh, fracture centers are generated by whose differences of strike and dip angles from
the same number as the fracture density estimated the A are under 10 and 20, respectively, the two
for the mesh and their locations are given using centers are connectable. Then, the B is used as a
random numbers. start point for the next search. Isolated fracture
(2) The strike data expressed in the above indicator that has no connectable center is eliminated.
set are transformed into the principal values. OK
uses these values and interpolates them at each Figure 8A depicts the result of the fracture-
fracture center localized by the step (1). The distribution modeling. To clarify distribution
interpolated principal values are then converted to characteristics in detail, only fractures estimated
values in the original coordinate system longer than 10 m are extracted as shown in Figure
corresponding to the indicator set. The strike 8B. The conspicuous directions of continuous
sector with the highest value is chosen as the fractures (N20-40E and N20-50W) agree with
estimated fracture direction. For example, a those of the measured fractures. Usefulness of the
calculated set (0.8, 0.3, 0.1, 0.4) assumes the proposed method can be proved in that it can
strike to be EW. estimate the ENE-WSW striking fractures almost
(3) With respect to dip direction, OK interpolates parallel to the boreholes, which are hardly appeared
the indicator-transformed dip directions of the in the boreholes, in the northern area and continuous
N
A B

15
15 m
m

Figure 8. A, Result of t wo-dimensional fracture- distribution modeling; and

B, distribution of extracted fractures estimated longer than 10 m.
 A B

5 mm
x
C Figure 9. A, Linear features extracted from an
Inada granite thin-section, using a modified
version of the STA; B, intersections of the
supposed lines and linear features; and C,
estimated linear-feature distribution.

attributable to cracks and trace the portions along

the cracks.
Figure 9A shows the extracted 3559 linear
features that partly correspond to microcracks. On
the image, virtual lines were set by modeling after
the arrangement of boreholes in the study area (Fig.
fractures in the zones without data. Another 1), and the intersections of the lines and linear
interesting feature is existence of the zones that features (Fig. 9B) were used to reconstruct the
contain the concatenated fractures passing through linear-feature distribution. Unit mesh size for
the study area. calcula ting the linear-feature density was defined as
0.250.25 mm2 so that the number of meshes,
VALIDATION OF THE PROPOSED METHOD 6060, is equal to that used for the fracture-density
map in Figure 7.
The proposed method of fracture-distribution As seen in the final result (Fig. 9C), the
modeling must be validated whether it allows to reconstructed pattern is similar to the linear-feature
estimate practical distribution and characterize distribution in that it highlights three dominant
hydraulic properties of the rock mass. For this directions, about 30, 120, and 160 counter-
purpose, a rock specimen over which fracture clockwise from the x-axis. The validness of the
distribution is entirely observed is desirable. We proposed method is demonstrated because these
prepared a thin-section of the Inada granite with dominant continuous line features are estimated in
1.51.5-cm2 size. The Inada granite situated in the zones with no data.
central Japan is known to contain many microcracks.
A modified version of the STA suitable for the EXTENSION TO THREE-DIMENSIONAL
microcrack analysis was used to extract linear FRACTURE- DISTRIBUTION MODELING
features from the digital image of thin-section. This
non-filtering technique can specify the portions in For more practical modeling of fracture
an image where the brightness largely change distribution, it is necessary to extend the proposed
 252-m level 248-m level N

Figure 10. Fracture distributions estimated at the 248 and 252-m levels. 15 m

method to be available in the three-dimensional Schmidts net. There are two sets of the simulated
space. Information on vertical continuities of the continuous fracture planes, the most remarkable
fractures cannot be obtained because the boreholes NW-SE striking set and the NNE-SSW striking set
are located at almost the same level. Therefore, an subordinate to it.
assumption is needed for a three-dimensional
modeling. We assumed that the fractures observed at APPLICATION TO CHARACTERIZATION
the 250-m level may be continuous only within the OF HYDRAULIC PROPERTIES
5-m depth range.
The fracture locations at the 245 to 255-m The three-dimensional fracture-distribution
levels were given at 1-m depth interval by extending model can be linked to the permeability tensor
the observed fractures along their strikes and dips. analysis to estimate permeability of the rock mass.
At each level, virtual boreholes were set on the According to Oda et al. (1987), the permeability
same locations at the 250-m level, and the tensor, k ij of degree three is expressed by
intersection points of the virtual boreholes and the
extended fractures were used for calculating the (
k ij = Pkk ij Pij )
fracture densities. Two-dimensional distribution
(1)
Pij = r t n n E (n, r , t )ddrdt
2 3
modelings were at first executed by the above- i j
4 0 0
mentioned method and then, fracture planes were
constructed by connecting simulated fractures at
adjacent levels based on the similarity of strike and where ij : delta function, : a constant expressing
dip and the nearness of locations. For example, the the continuity of fractures, r: diameter of fracture, t:
distributions conje ctured at the 248 and 252-m hydraulic aperture of fracture, n: normal vector of
levels are shown in Figure 10, from which different fracture planeE(n, r, t): probabilistic density func-
patterns can be found with the level. tion with respect to n, r, and t, : volumetric
The modeling result is visualized in Figure 11, fracture density, and : solid angle. We defined the
which represents the locations of fracture planes as 1/12 and supposed that the t is 10-6 times the r.
longer than 5 m along the vertical direction. The The r was determined by converting a fracture plane
constructed planes are slightly undulate, but can be into a disc with the equivalent area.
considered to steeply dip at larger than 70. Figure Discretization of Eq. (1) was carried out and
12 shows the azimuthal frequencies of the planes three principal values (k1 , k 2 , k 3 ) and directions of
using the lower hemisphere projection of the the principal axes were calculated as shown in Table
 N
x
N

z y

5m
10 m

10 m

Perspective view from the east side.

10 m
N

5m Perspective view from the top side.

6 - 7 m Figure 11. Result of three-dimensional

fracture-distribution modeling, which
8m +
represents the locations of simulated
fracture planes longer than 5 m along
5m

the vertical direction.

y
x
10
m Table 1. Three principal values and directions of principal
m

axes calculated for the permeability tensor of Eq. (1).

10

z
Directions of principal axes
Principal values
Perspective view from the bottom side. (m2 )
Azimuth Dip
N
Frequency (%)
k1 1.0 10-17 N68 W 14
0
k2 4.3 10-18 N20 E 4
2.3
k3 2.4 10-18 N83 W 76
W E 4.7

7.0 compared to the present study area. Considering the

data, condition of the rock mass, and the
9.4+
permeability data reported for hard granitic rocks in
S many areas (e.g., Brace, 1980; Brace, 1984; Clauser,
Figure 12. Azimuthal frequencies of the planes using the 1992), the first principal value , 10-17 m2 can be
lower hemisphere projection of the Schmidts net. regarded as a pla usible magnitude.

1. The principal axes generally correspond to the CONCLUSION

two dominant strikes for the k 1 and k2 and the dip
direction of the predominant fractures for the k3. By analyzing the borehole fracture data in the
The magnitude of principal values cannot be Kamaishi mine, spatial correlations were clarified
validated because there is no measured permeability on the three attributes of fractures, fracture densities
data in the study area. However, hydraulic well tests along the boreholes, appearance relations of fracture
were executed in other areas in the Kamaishi mine strikes, and dip directions with respect to northward
where the fractures are much developed as or southward dip. The experimental semivariograms
of these attributes could be approximated by Geological Society, London, Special publications, 127 p.
spherical, exponential, and spherical models, Davy, P., Sornette, A., and Sornette, D. (1992)
respectively. The SGS was proved to be effective Experimental discovery of scaling laws relating fractal
for estimating fracture densities in the sparse data dimensions and the length distribution exponent of
zones. To consider the azimuthal correlation of fault systems, Geophys. Res. Letters, v. 19, no. 4, p.
fracture pairs, we transformed the strike and dip 361-363.
data into the indicators. Applying the principal Dawers, N. H., Anders, M. H., and Scholz, C. H. (1993)
component analysis and ordinary krig ing to the Growth of normal faults: Displacement-length scaling,
indicators, two- and three-dimensional models of Geology, v. 21, p. 1107-1110.
fracture distributions, which incorporate both the Gringarten, E. (1996) 3-D geometric description of
azimutha l and positional information of the fracture fractured reservoirs, Math. Geology, v. 28, no. 7, p.
data, could be constructed. 881-893.
The three-dimensional model was effectively Journal, A. G. (1989) Fundamentals of geostatistics in
linked to the permeability tensor analysis for five lessons, Short course in geology: v. 8, AGU, 40 p.
estimating the permeability and principal axes. Two Koike, K., Nagano, S., and Ohmi, M. (1995) Lineament
characteristics were detected by it in that the analysis of satellite images using a Segment Tracing
permeability along the first principal axis is about Algorithm (STA), Computers & Geosciences, v. 21, no.
10-17 m2 and the major axes correspond to the 9, p. 1091-1104.
dominant azimuths of simulated fracture planes. The Koike, K., and Kaneko, K. (1999) Characterization and
estimated permeability is plausible by considering modeling of fracture distribution in rock mass using
the condition of rock mass and the hydraulic test fractal theory, Geothermal Science & Technology, v. 6,
data obtained at other sites in the mine. p. 43-62.
Long, J. C. S., and Billaux, D. M. (1987) From field data
Acknowledg ments: The authors wish to express their to fracture network modeling: An example
grateful thanks to the Japan Nuclear Cycle Development incorporating spatial structure, Water Resources Res., v.
Institute for providing the in situ measurement results in 23, no. 7, p. 1201-1216.
the Kamaishi mine. National Research Council (1996) Rock fractures and
fluid flow, contemporary understanding and applica-
REFERENCES tions, National Academic Press, Washington, D.C.,
Adler, P. M., and Thovert, J. -F. (1999) Fractures and 551 p.
fracture networks, Kluwer Academic Publishers, Oda, M., Hatsuyama, Y., and Ohnishi, Y. (1987)
Dordrecht, The Netherlands, 429 p. Nu merical experiments on permeability tensor and its
Brace, W. F. (1980) Permeability of crystalline and application to jointed granite at Stripa mine, Sweden, J.
argillaceous rocks, Int. J. Rock Mech. Min. Sci. & Geophys. Res., v. 92, no. B8, p. 8037-8048.
Geomech. Abstr., v. 17, 221-251. Vignes-Adler, M., Le Page, A., and Adler, P. M. (1991)
Brace, W. F. (1984) Permeability of crystalline rocks: Fractal analysis of fracturing in two African regions,
New in situ measurements, J. Geophys. Res., v. 89, no. from satellite imagery to ground scale, Tectonophysics,
B6, p. 4327-4330. v. 196, p. 69-85.
Chils, J. P. (1988) Fractal and geostatistical methods for Watterson, J., Walsh, J. J., Gillespie, P. A., and EastonS.
modeling of a fracture network, Math. Geology, v. 20, (1996) Scaling systematics of fault sizes on a
no. 6, p. 631-654. large-scale range fault map, J. Struct. Geol., v. 18, nos.
Clauser, C. (1992) Permeability of crystalline rocks, EOS, 2/3, p. 199-214.
v. 73, no. 21, p. 233-240. Young, D. S. (1987) Indicator kriging for unit vectors:
Coward, M. P., Daltaban, T. S., and Johnson, H. (eds.) Rock joint orientations, Math. Geology, v. 19, no. 6, p.
(1998) Structural geology in reservoir characterization, 481-501.