- Email: [email protected]
Relation between the earthquake statistics and fault patterns, and fractals and percolation
Physica A 191 (1992) North-Holland
57-68
Relation between the earthquake statistics and fault patterns, and fractals and percolation Muhammad Sahimi”‘“, Michelle C. Robertsonb and Charles G. Sammisb “Department of Chemical Engineering, University of Southern California, Los Angeles, CA 90089, USA’ bDepartment of Geological Sciences, University of Southern California, Los Angeles, CA 90089, USA ‘HLRZ Supercomputer Center, clo KFA Jiilich, W-51 70 Jiilich 1, Germany
Using computer simulations and analyzing experimental data, we investigate the structure of fault patterns in heterogeneous rock, and the distribution of earthquake hypocentral locations. At large length scales, the fracture networks of rock are fractal with a fractal dimension D = 1.9 and 2.5, in 2D and 3D, respectively. We present a computer simulation model that explains these results in terms of percolation fractals. Three-dimensional fractal analysis of regional hypocentral locations yields a fractal dimension of about 1.8. This result can be explained if earthquakes are distributed on the backbone of the fractal network of fault structure with a fractal dimension equal to that of the percolation backbone. Thus, percolation provides a unified theory for the geometry of fault patterns and the spatial distribution of earthquakes.
1. Introduction
Two related phenomena in most natural rock masses are earthquakes, and the nucleation and propagation of fractures. Most rock masses contain fractures, in the form of a complex and interconnected network. The presence of such fractures is crucial to the economics of oil recovery processes from underground reservoirs, generation of heat and vapor from geothermal reservoirs for use in power plants, and the development of groundwater resources, as the fractures provide high permeability paths for fluid flow in reservoirs that are otherwise of very low permeabilities and porosities. Despite the obvious significance of fractures, the characterization of fractured rock is not as well developed as that of porous media (without fractures) and remains an active research field (for a comprehensive review, see ref. [l]). A clear consensus about the structure of fracture networks, which is particularly important for modeling of fluid flow in fractured rock has yet to emerge. 1Permanent
address.
0378-4371/92/$05.00
0
1992 - Elsevier
Science
Publishers
B.V. All rights
reserved
58
M. Sahimi et al. I Earthquake
statistics and fault patterns
Earthquakes are the result of a series of complex phenomena involving the interaction between stress concentrations and fluid flow. Study of the dynamics of earthquake faults is also a very active area of research [2-81. It has been suggested [3-51 that earthquakes are related to self-organized critical phenomena (SOCP), in that they are the product of a dynamical many-body system that reaches a critical state without a need to fine-tune its parameters. Such systems drive themselves into a statistically stationary state characterized by spatial and temporal correlation functions that follow power-law behavior, which means that they have no intrinsic length or time scale and, therefore, are at a critical state. Despite such proposals, a clear picture of how earthquakes develop, and a geometrical interpretation of the distribution of earthquake hypocenters, has not emerged. In this paper, we use computer simulations, rock fracture patterns, and regional distributions of seismicity to argue that a single phenomenon, namely, a percolation process [9], provides a unified theory for both the geometry of fault patterns and the spatial distribution of earthquakes. In particular, we argue that the fault patterns in rock form a fractal percolation cluster, while earthquake hypocenters are distributed on the backbone of the fault network, i.e. the multiply connected part of the cluster that transmits stress throughout the network. Thus, both fault patterns and earthquakes may be given a purely geometrical interpretation in terms of the well-understood percolation processes. The plan of this paper is as follows. In the next section we summarize the previous work in the area of fault patterns. In section 3 we briefly analyze the data on fractured reservoirs. In section 4 we outline a discrete model of fracture that can help explain the experimental data on fracture networks. Section 5 presents our analysis of the spatial distribution of earthquake hypocenters. The paper is summarized in section 6 where we comment on the possible relation between our interpretation of earthquake data and other models of this phenomenon. A brief account of our work was recently given as a letter [lo].
2. Review of past work on fracture networks in rock Over the past decade several authors have discussed the application of percolation theory to modeling of fracture networks in rock. However, these authors assumed that percolation describes the statistics of fracture networks in rock without actually proving it. For example, Chelidze [ll] considered both finite and infinite percolation networks, and suggested a percolation model for fracturing processes. His model contained several parameters such as the total
M. Sahimi et al.
I Earthquake
statistics and fault patterns
59
number of fractures, the coordination number of the fracture network and, for finite systems of fractures, their length. Englman et al. [12] ConsideredJlow in fractured rocks by modeling the fracture network as a random percolation cluster of fractures in a plane, in which fractures were assumed to extend infinitely in the direction orthogonal to the plane. Long and Witherspoon [13] developed two-dimensional models of fracture networks in which fractures were represented by long sticks to which hydraulic permeabilities, selected at random from a probability distribution, were assigned. They pointed out that their results may be compared with the percolation results of Englman et al. [12]. Charlaix et al. [14] represented a fractured rock by a random array of flat disks with a broad distribution of apertures. Using the ideas of Ambegaokar et al. [15], they argued that the permeability of the fracture network is determined by only that portion of the network which is made of large fractures. Thus, the permeability of the smaller fractures can be set to be zero, which reduces the problem to a percolation process. Similar ideas had already been used by Katz and Thompson [16] for estimating the permeability of porous media. The first hint about the applicability of percolation to modeling of fractured rocks can be found in the work of Long et al. [17]. These authors used a simulated annealing method, similar to those used in optimization processes and statistical mechanics of spin glasses [18], to find an “optimized” structure for fracture networks in rock, given a limited amount of data on the rock properties. They found that in all cases they examined, the optimized fracture network was “barely” connected, which implies that the fracture network was a percolation cluster very near the percolation threshold. However, Long et al. [17] did not attempt to explain this result in terms of percolation or any other statistical theory (in fact, they did not even mention percolation), nor did they attempt to explain why brittle fracture processes, which presumably give rise to fracture networks of rocks, should have any similarity with percolation. This is one of the major goals of our paper.
3. Analysis of fracture network data We start with the analysis of fracture patterns. We focus here on the structure of the fracture network of the Geysers geothermal field in Northeast California. This field, from which heat and vapor are extracted for use in power plants, covers an area of more than 35 000 acres and represents one of the most significant geothermal fields in the world. To investigate the structure of the fracture network, we first look at the surface of an outcrop of the reservoir rock (fig. 1) using the pavement method [19,20] of clearing a
M. Sahimi et al. I Earthquake
60
Fig.
1. A typical
feature
statistics and fault patterns
surface
of the Geysers
field.
subplanar surface and mapping the fracture surface in order to measure connectivity, density and other properties of the network. We find that the fractured pavements have a fractal structure, and that in a box of length L the number N of fractures is given by
N-L-D, where D is the fractal dimension of the fracture network. Typical results are shown in fig. 2. We find that for large length scales, D = 1.9. Qualitatively similar results are found for three-dimensional samples, namely, for small scales we find D = 1.8, whereas as L increases, the effective value of D also increases. For the largest sample we have examined, D = 2.1. On the other hand, somewhat different results were found for Yucca Mountain, Nevada [19], which was investigated for the U.S. Department of Energy as a potential underground repository for high-level radioactive wastes, and for the Monterey
M. Sahimi et al. I Earthquake
statistics and fault patterns
61
Geysers Regional
100 10
100
1000
length (meters) Fig. 2. Number
of fractures
N vs. length
scale L for fracture
surfaces.
Formation [20] in Central California, which contains a considerable amount of oil. These studies found that for small samples, D = 1.7. We point out that several previous studies of fault patterns around the world had found results similar to ours (for a detailed discussion, see ref. [21]). It has been generally found that at large length scales (order of 100 km) the fault pattern is fractal with a fractal dimension of about 1.9 in 2D, whereas at smaller length scales the fractal dimension has been found to be about 1.6-1.7. For example, Hirata [22] studied the self-similarity of fault geometry in various parts of Japan. For a set of faults read from maps a fractal dimension of about 1.6, in the size range 2-20 km, was obtained. One can also interpret coda measurements in terms of a fractal geometry. The coda is the late acoustic or seismic signal that one measures in a given remote seismic station, several seconds after one receives the first seismic signal that has emanated from an earthquake. It is generally believed that the coda reflects scattering of the seismic wave by inhomogeneities within the lithosphere. It is clear that such measurements provide information about the three-dimensional structure of rock. Based on this interpretation, Wu and Aki [23] and Sato [24] pointed out evidence of a fractal structure of inhomogeneities in the lithospace, with a fractal dimension of about 2.5, which is the same as that of 3D percolation. Finally, at smaller length scales, those that are of the order l-10 km, a single major fault is not a straight line, but can be characterized by a fractal dimension of about 1.2 [25,26]. It is not, however, clear to us why the fractal dimension of the fracture surfaces we analyzed is not simply one less than that of the three-dimensional fracture network. We now propose a model that can explain all these data.
62
M. Sahimi et al. I Earthquake
statistics and fault patterns
4. A discrete model of fracture
What is the theoretical explanation for these results? At small length scales, rock, although microscopically disordered, is macroscopically homogeneous. The nucleation and propagation of fractures in rock do not take place at random, but depend on the stress field in the system. Once a crack nucleates in rock, stress enhancement at its tip is larger than at any other point of the medium and, therefore, the next microcrack almost surely develops at its tip. Hence, at small length scales the microcracks are all connected and cluster together with almost no dead-end microcracks. This is consistent with the experimental observations [19,20,27] that fractures often seem to be clustered together at different parts of rock, separated by large intervals of unfractured rock. However, at larger length scales rock is macroscopically heterogeneous with spatially varying properties. In fact, it has been shown that [1,28] the permeability and other properties of heterogeneous rock are spatially and fractally distributed quantities. The growth of fractures in a heterogeneous medium is very different from that in a homogeneous medium described above. Crack growth continues only as long as a much stronger region of the rock is not encountered. But because of the spatially varying properties of rock, the growing crack will eventually encounter a strong region of the rock and its growth stops, in which case another microcrack in another region of the rock nucleates and starts to grow. The growth of the new crack also stops when it encounters another strong region of the rock, and so on. This has the consequence that, viewed from a large scale, the nucleation and growth of cracks in different regions of the rock take place essentially at random, since they are strongly dependent on the spatial distribution of the rock properties. Therefore, the fractal dimension of the fracture network at large length scales should be the same as that of percolation clusters, since in this process microcracks are created at random (i.e., bonds of a network are broken at random). Indeed, our 2D experimental result for large length scales, D = 1.9, is precisely the same as that of 2D percolation clusters, D = 91148 = 1.9, and our 3D experimental results with larger and larger samples indicate clearly a trend towards D = 2.5, the fractal dimension of 3D percolation clusters. To further test our ideas, we employed our discrete model of brittle fracture developed recently [29-311 to simulate fracture of a macroscopically heterogeneous medium. It has already been established that fracture in microscopically disordered, macroscopically homogeneous media can be successfully simulated with this model. Indeed, in two dimensions this model yields [29,31,32] fractal fracture patterns. If we analyze the statistics of all cracks of such patterns, we find that they have a fractal dimension, D 2: 1.7, in complete agreement with
M. Sahimi et al. I Earthquake statisticsand fault patterns
63
the experimental results. If, on the other hand, we analyze only the fractal behavior of the single fracture that spans the system, we find that it has a fractal dimension of about 1.2, consistent with the results of Okubo and Aki [25] and Aviles et al. [26] for a single major fault. To simulate fracture in a heterogeneous medium, we consider a 90 X 90 triangular network with periodic boundary condition in one direction. Every site of the network is characterized sites are by the displacement vector ui = (uiX, uiY), and nearest-neighbor connected by springs. We consider here the case of brittle fracture for which a linear approximation is valid up to a threshold (defined below). The nodal displacement ui is computed by minimizing the elastic energy E with respect to ui, where
E =
’ 2 2 (ij)
[(u; - ~j)
.R,]$, +
’
2
ZZ(Sqik)‘gijgik . (jik)
Here LYand /? are the central and bond-bending or angle-changing force constants, respectively, R, is a unit vector from site i to site j, g, is the elastic constant of the spring between i and j, and ( jik) indicates that the sum is over all triplets in which the bonds j-i and i-k form an angle whose vertex is at i. Although one may need a more complex form of the elastic energy 1321 for describing all properties of fracture networks, we believe that eq. (2) is adequate for studying the fructul properties of the system and, in particular, its fractal dimension, since D is often independent of microscopic features of E. To make the network macroscopically heterogeneous, we divide it into 25 blocks of size 18 x 18. Within each block g, was set to be the same for all springs, but for each block g,j was selected from the distribution Pi(g) = where 0 < y1 < 1. This is a broad distribution and, thus, the elastic (1 - r,)P, moduli and other properties of the entire network are also spatially dependent quantities and are broadly distributed, as required [1,28]. We then introduce a threshold value I, for the length of a spring, which is selected according to the distribution, P,(I,) = (1 - y,)lLy2, with 0 < y2 < 1. We found that our results are insensitive to the values of y1 and K. We then initiate the failure process by applying a fixed external strain to the network in a given direction, calculating the nodal displacements ui, and breaking the spring for which l- I, is maximum, where 1 is the current length of the spring in the strained network. After a spring is broken, we recalculate the nodal displacements ui for the new configuration of the network, select the next spring to break, and so on. If the external strain is not large enough to break any new spring, we gradually increase it. The simulation continues until the network finally becomes macroscopically disconnected and a macroscopic fracture network is formed. We then
64
M. Sahimi et al, I Earthquake
I
OO Fig. 3. Number of fracture.
of fractures
statistics and fault patterns
I
1
N vs. the radius
I
3
J 4
I&
r. obtained
with the two-dimensional
discrete
model
calculate the fractal dimension of the fracture network by plotting log N vs. log r, where N is the number of broken springs and r is the radius of the growing crack. The results are shown in fig. 3 (where yI = -yZ= 0.6, and /3/a = 0.04), from which we obtain D 2: 1.85, in very good agreement with our results discussed above, which also strongly supports our view that at large length scales the fracture network in heterogeneous rocks should have the structure of a percolation cluster. The fracture network cannot have the structure of a percolation cluster above p, because, as soon as the fracture network becomes sample-spanning, large scale stress concentrations are relieved and the stress distribution becomes localized. As a result, large scale fracture growth can no longer occur. The fact that the model does reproduce the data has two implications. The first is that the measured fractal dimensions are not merely effective values representing possible crossovers between small, intermediate and large length scales. These fractal dimensions do have clear physical interpretations, as discussed above. The second is that, in a highly heterogeneous medium, the generally accepted idea that stress enhancement at the tip of a growing crack always causes it to grow further may no longer be valid. Thus, crack nucleation and growth may become independent of the stress field in a strongly heterogeneous system. 5. Analysis
of seismic data
We now analyze seismic data and interpret the results. Four seismic data sets from four different regions in Southern California were selected for our
M. Sahimi et al. I Earthquake
65
statistics and fault patterns
Table I The parameters of the four sets of data used in the analysis, where m represents the magnitude of the earthquakes. Region
Volume (km)
Number of events
m
Time span
Location error (km)
SA-EL Parkfield Whittier Uoland
100 x 10 x 17 x 18 x
1859 717 224 291
>2 >O.Ol >I >l
1981-1990 1988-1990 10/1987-1211987 2/1990-6/1990
4-5 0.1 l-2 l-2
125 x 20 18 x 15 18 x 15 18 x 15
analysis. Table I compares the parameters for each data set. Each event represents an earthquake of size m. For each set, the fractal dimension was estimated 1331 by counting the minimum number of cubes N of edge length L required to completely cover the data set at various length scales L. For each set, the minimum L was taken to be the location-error (in kilometers) of the set. Fig. 4 presents the results, and it is clear that in all cases the spatial distribution of hypocenters is fractal over at least one order of magnitude, with a fractal dimension of about 1.8. What is the interpretation of this value of D? Geologic features are rarely simple planes, and in fact most faults are comprised of many shear fractures which are visible at a number of scales. Therefore, given our results for fracture networks, it is clear that, at large length scales, the fault patterns should make a percolation cluster with a fractal dimension, D = 2.5 in 3D. The locations of earthquake hypocenters are of course on the network of faults. However, these locations have to belong to that part of the network that is active, i.e. the part through which stress and deformation are transmitted. That is, the locations of earthquake hypocenters have to belong to the backbone of the fault network. Indeed, our result, D = 1.8 is perfectly consistent with the fractal dimension D, of the backbone of 1000
1
Fig. 4. Variations of the minimum number of occupied cubes N of edge length L to completely cover the data set at various length scales L for all four sets of data. The slopes of all data sets are about =1.80.
66
M. Sahimi et al. I Earthquake
statistics and fault patterns
percolation clusters [34], which is also about 1.8. This is the first time that the statistics of earthquakes has been given a clear geometrical interpretation. We remark that Kagan and Knopoff [35] and Andrews [36] had already pointed out the relevance of fractal concepts to the distribution of earthquakes (see below). The fractal dimension of the distribution of hypocentral locations appears to be universal, at least for the four sets of data that we analyzed. Kagan and Knopoff [35] studied the spatial distribution of earthquakes and found, D 2: 2, which they interpreted as meaning that earthquakes occur on planes. But they also found that if they explore larger depths and areas, the estimated fractal dimension decreases. On the other hand, if the hypocentral locations are plotted in 3D space, and rotated around all three axes, there does not seem to be any one plane about which the earthquakes cluster. The San AndreasElsinore (SA-EL) earthquakes are located over a very wide zone of deformation, between and inclusive of two major faults, whereas the Parkfield earthquakes are located over a relatively thin zone of deformation along the SA fault, yet both have a fractal dimension of about 1.8. This universality is also consistent with the fact that the fractal dimensions of the backbone of percolation clusters is also universal, and lends further support to our interpretation of the data. 3D
6. Discussion and conclusions We conclude the paper by making a few remarks. The number N, of earthquakes of magnitude greater than m is related to the energy release E, through, dN,ldE, - rnm7, where 1< r < 2. It has been argued that [37] 7 is determined by the geometric structure of the fault network. In view of our results, it is likely that T is related to the fractal dimension D, of the backbone of the fault network or percolation clusters. It has also been argued [4,37] that the occurrence of earthquakes and their spatial distribution is closely related to the propagation of fractures and the structure of their network in the rock, and our results establish this connection firmly. Kagan and Knopoff [35,37] also proposed a model of earthquake occurrence which appears to provide relatively accurate predictions. In their model it is assumed that each earthquake generates additional shocks with a probability that depends on time as tFH. The analysis of the data showed that a value of 19close to 1.8 seems to simulate the statistical properties of intermediate and deep shocks, and it is interesting to note that 8 = D, = 1.8, the same as the fractal dimension of the network of hypocentral locations. As mentioned above, several authors [3-51 have attempted to explain earthquakes based on the theory of SOCP. However, the
M. Sahimi et al. I Earthquake statisticsand fault patterns
67
application of this theory to the interpretation of earthquake statistics is appropriate if the system is macroscopically homogeneous whereas, as our results indicate, it may be more appropriate to develop the theory for a macroscopically heterogeneous system: the analysis of fracture data discussed above indicates that the fractal dimensions of the fracture networks in homogeneous and heterogeneous rocks can be quite different. Simulation of earthquake statistics in a macroscopically homogeneous system and the interpretation of the results based on SOCP yielded [3] r = 1.35. It may be that the same simulation in a heterogeneous system could yield an exponent 7 closely related to 4~1.8. Finally, in a previous paper [38] we proposed that the universality classes of fracture of disordered media can be characterized by universal fixed points at the fracture point. We suggested that there are two such fixed points, one of which corresponds to systems in which fracture depends on the stress field, while the other one is for systems in which fracture is independent of the stress field. Our results in the present paper indicate clearly that these fixed points correspond to fracture in microscopically disordered and macroscopically disordered media, respectively.
Acknowledgements
This work was supported in part by the U.S. Department of Energy, the Petroleum Research Fund, administered by the American Chemical Society, the USC Center for the Study of Fractured Reservoirs, and the San Diego Supercomputer Center. This paper was prepared while one of us (M.S.) was visiting the HLRZ Supercomputer Center at KFA Julich, as an Alexander von Humboldt Foundation Research Fellow. He would like to thank the Center for warm hospitality, Hans Herrmann for useful discussions, and the Foundation for financial support.
References [l] M. Sahimi, Rev. Mod. Phys., to be published. [2] S.R. Brown, C.H. Scholz and J.B. Rundle, Geophys. Res. Lett. 18 (1991) 215. [3] P. Bak, C. Tang and K. Wiesenfeld, Phys. Rev. Len. 59 (1987) 381; P. Bak and C. Tang, J. Geophys. Res. 94 (1989) 15635. [4] A. Sornette and D. Sornette, Europhys. Lett. 9 (1989) 197. [5] K. Ito and M. Matsuzaki, J. Geophys. Res. 95 (1990) 6853. [6] J.M. Carlson and J.S. Langer, Phys. Rev. Lett. 62 (1989) 2632; Phys. Rev. A 40 (1989) 6470. [7] H. Nakanishi, Phys. Rev. A 43 (1991) 6613. [8] Z. Olami, H.J.S. Feder and K. Christensen, Phys. Rev. Lett. 68 (1992) 1244.
68
M. Sahimi et al. I Earthquake statistics and fault patterns
[9] D. Stauffer and A. Aharony, Introduction to Percolation Theory, 2nd ed. (Taylor and Francis, London, 1992). [lo] M. Sahimi, M.C. Robertson and C.G. Sammis, Phys. Rev. Lett., to be published. [ll] T.L. Chelidze, Phys. Earth Planet Int. 28 (1982) 93. [12] R. Englman, Y. Gur and Z. Jaeger, J. Appl. Mech. 50 (1983) 707. [13] J.C.S. Long and P.A. Witherspoon, J. Geophys. Res. 90 (B4) (1985) 3087. [14] E. Charlaix, E. Guyon and S. Roux, Trans. Porous Media 2 (1987) 31. [15] V. Ambegaokar, B.I. Halperin and J.S. Langer, Phys. Rev. B 4 (1971) 2612. [16] A.J. Katz and A.H. Thompson, J. Geophys. Res. 92 (1987) 599. [17] J.C.S. Long, K. Karasaki, A. Davey, J. Peterson, M. Landsfeld, J. Kemeny and S. Martel, Int. J. Rock. Mech. Min. Sci. Geomech. Abstr. 28 (1991) 121. [18] S. Kirkpatrick, C.D. Gelatt and M.P. Vecchi, Science 220 (1983) 671. [19] C.C. Barton and P.A. Hsieh, Physical and Hydrological-Flow Properties of Fractures, Guidebook T385 (American Geophysical Union, Las Vegas, NV, 1989). [20] C.G. Sammis, L. An and I. Ershaghi, unpublished. [21] D. Sornette, in: Spontaneous Formation of Space-Time Structures and Criticality, T. Ristc, ed., NATO AS1 (Kluwer, Dordrecht, 1991). [22] T. Hirata, Pure Appl. Geophys. 131 (1989) 157. [23] R.-S. Wu and K. Aki, Pure Appl. Geophys. 123 (1985) 806. [24] H. Sato, Pure Appl. Geophys. 128 (1988) 43. [25] P.G. Okubo and K. Aki, J. Geophys. Res. 92 (1987) 345. [26] C.A. Aviles, C.H. Scholz and J. Boatwright, J. Geophys. Res. 92 (1987) 331. in: Reservoir Assessment of the Geysers [27] R.P. Thomas, R.H. Chapman and H. Dykstra, Geothermal Field, Publication TR27 (California Department of Conservation, 1981). [28] T.A. Hewett and R.A. Behrens, SPE Formation Evaluation 5 (1990) 217. [29] M. Sahimi and J.D. Goddard, Phys. Rev. B 33 (1986) 7848. [30] S. Arbabi and M. Sahimi, Phys. Rev. B 41 (1990) 772. [31] M. Sahimi and S. Arbabi, Phys. Rev. B, to be published. [32] L. de Arcangeli, A. Hansen, H.J. Herrmann and S. Roux, Phys. Rev. B 40 (1989) 877. [33] C.G. Sammis and R. Biegel, Pure Appl. Geophys. 131 (1989) 255. [34] M. Sahimi, J. Phys. A 17 (1984) 3073. [35] Y.Y. Kagan and L. Knopoff, J. Geophys. Res. 86 (1981) 2853. [36] D.J. Andrews, J. Geophys. Res. 86 (1981) 10821. (371 Y.Y. Kagan and L. Knopoff, Science 236 (1987) 1563. [38] M. Sahimi and S. Arbabi, Phys. Rev. Lett. 68 (1992) 608.
57-68
Relation between the earthquake statistics and fault patterns, and fractals and percolation Muhammad Sahimi”‘“, Michelle C. Robertsonb and Charles G. Sammisb “Department of Chemical Engineering, University of Southern California, Los Angeles, CA 90089, USA’ bDepartment of Geological Sciences, University of Southern California, Los Angeles, CA 90089, USA ‘HLRZ Supercomputer Center, clo KFA Jiilich, W-51 70 Jiilich 1, Germany
Using computer simulations and analyzing experimental data, we investigate the structure of fault patterns in heterogeneous rock, and the distribution of earthquake hypocentral locations. At large length scales, the fracture networks of rock are fractal with a fractal dimension D = 1.9 and 2.5, in 2D and 3D, respectively. We present a computer simulation model that explains these results in terms of percolation fractals. Three-dimensional fractal analysis of regional hypocentral locations yields a fractal dimension of about 1.8. This result can be explained if earthquakes are distributed on the backbone of the fractal network of fault structure with a fractal dimension equal to that of the percolation backbone. Thus, percolation provides a unified theory for the geometry of fault patterns and the spatial distribution of earthquakes.
1. Introduction
Two related phenomena in most natural rock masses are earthquakes, and the nucleation and propagation of fractures. Most rock masses contain fractures, in the form of a complex and interconnected network. The presence of such fractures is crucial to the economics of oil recovery processes from underground reservoirs, generation of heat and vapor from geothermal reservoirs for use in power plants, and the development of groundwater resources, as the fractures provide high permeability paths for fluid flow in reservoirs that are otherwise of very low permeabilities and porosities. Despite the obvious significance of fractures, the characterization of fractured rock is not as well developed as that of porous media (without fractures) and remains an active research field (for a comprehensive review, see ref. [l]). A clear consensus about the structure of fracture networks, which is particularly important for modeling of fluid flow in fractured rock has yet to emerge. 1Permanent
address.
0378-4371/92/$05.00
0
1992 - Elsevier
Science
Publishers
B.V. All rights
reserved
58
M. Sahimi et al. I Earthquake
statistics and fault patterns
Earthquakes are the result of a series of complex phenomena involving the interaction between stress concentrations and fluid flow. Study of the dynamics of earthquake faults is also a very active area of research [2-81. It has been suggested [3-51 that earthquakes are related to self-organized critical phenomena (SOCP), in that they are the product of a dynamical many-body system that reaches a critical state without a need to fine-tune its parameters. Such systems drive themselves into a statistically stationary state characterized by spatial and temporal correlation functions that follow power-law behavior, which means that they have no intrinsic length or time scale and, therefore, are at a critical state. Despite such proposals, a clear picture of how earthquakes develop, and a geometrical interpretation of the distribution of earthquake hypocenters, has not emerged. In this paper, we use computer simulations, rock fracture patterns, and regional distributions of seismicity to argue that a single phenomenon, namely, a percolation process [9], provides a unified theory for both the geometry of fault patterns and the spatial distribution of earthquakes. In particular, we argue that the fault patterns in rock form a fractal percolation cluster, while earthquake hypocenters are distributed on the backbone of the fault network, i.e. the multiply connected part of the cluster that transmits stress throughout the network. Thus, both fault patterns and earthquakes may be given a purely geometrical interpretation in terms of the well-understood percolation processes. The plan of this paper is as follows. In the next section we summarize the previous work in the area of fault patterns. In section 3 we briefly analyze the data on fractured reservoirs. In section 4 we outline a discrete model of fracture that can help explain the experimental data on fracture networks. Section 5 presents our analysis of the spatial distribution of earthquake hypocenters. The paper is summarized in section 6 where we comment on the possible relation between our interpretation of earthquake data and other models of this phenomenon. A brief account of our work was recently given as a letter [lo].
2. Review of past work on fracture networks in rock Over the past decade several authors have discussed the application of percolation theory to modeling of fracture networks in rock. However, these authors assumed that percolation describes the statistics of fracture networks in rock without actually proving it. For example, Chelidze [ll] considered both finite and infinite percolation networks, and suggested a percolation model for fracturing processes. His model contained several parameters such as the total
M. Sahimi et al.
I Earthquake
statistics and fault patterns
59
number of fractures, the coordination number of the fracture network and, for finite systems of fractures, their length. Englman et al. [12] ConsideredJlow in fractured rocks by modeling the fracture network as a random percolation cluster of fractures in a plane, in which fractures were assumed to extend infinitely in the direction orthogonal to the plane. Long and Witherspoon [13] developed two-dimensional models of fracture networks in which fractures were represented by long sticks to which hydraulic permeabilities, selected at random from a probability distribution, were assigned. They pointed out that their results may be compared with the percolation results of Englman et al. [12]. Charlaix et al. [14] represented a fractured rock by a random array of flat disks with a broad distribution of apertures. Using the ideas of Ambegaokar et al. [15], they argued that the permeability of the fracture network is determined by only that portion of the network which is made of large fractures. Thus, the permeability of the smaller fractures can be set to be zero, which reduces the problem to a percolation process. Similar ideas had already been used by Katz and Thompson [16] for estimating the permeability of porous media. The first hint about the applicability of percolation to modeling of fractured rocks can be found in the work of Long et al. [17]. These authors used a simulated annealing method, similar to those used in optimization processes and statistical mechanics of spin glasses [18], to find an “optimized” structure for fracture networks in rock, given a limited amount of data on the rock properties. They found that in all cases they examined, the optimized fracture network was “barely” connected, which implies that the fracture network was a percolation cluster very near the percolation threshold. However, Long et al. [17] did not attempt to explain this result in terms of percolation or any other statistical theory (in fact, they did not even mention percolation), nor did they attempt to explain why brittle fracture processes, which presumably give rise to fracture networks of rocks, should have any similarity with percolation. This is one of the major goals of our paper.
3. Analysis of fracture network data We start with the analysis of fracture patterns. We focus here on the structure of the fracture network of the Geysers geothermal field in Northeast California. This field, from which heat and vapor are extracted for use in power plants, covers an area of more than 35 000 acres and represents one of the most significant geothermal fields in the world. To investigate the structure of the fracture network, we first look at the surface of an outcrop of the reservoir rock (fig. 1) using the pavement method [19,20] of clearing a
M. Sahimi et al. I Earthquake
60
Fig.
1. A typical
feature
statistics and fault patterns
surface
of the Geysers
field.
subplanar surface and mapping the fracture surface in order to measure connectivity, density and other properties of the network. We find that the fractured pavements have a fractal structure, and that in a box of length L the number N of fractures is given by
N-L-D, where D is the fractal dimension of the fracture network. Typical results are shown in fig. 2. We find that for large length scales, D = 1.9. Qualitatively similar results are found for three-dimensional samples, namely, for small scales we find D = 1.8, whereas as L increases, the effective value of D also increases. For the largest sample we have examined, D = 2.1. On the other hand, somewhat different results were found for Yucca Mountain, Nevada [19], which was investigated for the U.S. Department of Energy as a potential underground repository for high-level radioactive wastes, and for the Monterey
M. Sahimi et al. I Earthquake
statistics and fault patterns
61
Geysers Regional
100 10
100
1000
length (meters) Fig. 2. Number
of fractures
N vs. length
scale L for fracture
surfaces.
Formation [20] in Central California, which contains a considerable amount of oil. These studies found that for small samples, D = 1.7. We point out that several previous studies of fault patterns around the world had found results similar to ours (for a detailed discussion, see ref. [21]). It has been generally found that at large length scales (order of 100 km) the fault pattern is fractal with a fractal dimension of about 1.9 in 2D, whereas at smaller length scales the fractal dimension has been found to be about 1.6-1.7. For example, Hirata [22] studied the self-similarity of fault geometry in various parts of Japan. For a set of faults read from maps a fractal dimension of about 1.6, in the size range 2-20 km, was obtained. One can also interpret coda measurements in terms of a fractal geometry. The coda is the late acoustic or seismic signal that one measures in a given remote seismic station, several seconds after one receives the first seismic signal that has emanated from an earthquake. It is generally believed that the coda reflects scattering of the seismic wave by inhomogeneities within the lithosphere. It is clear that such measurements provide information about the three-dimensional structure of rock. Based on this interpretation, Wu and Aki [23] and Sato [24] pointed out evidence of a fractal structure of inhomogeneities in the lithospace, with a fractal dimension of about 2.5, which is the same as that of 3D percolation. Finally, at smaller length scales, those that are of the order l-10 km, a single major fault is not a straight line, but can be characterized by a fractal dimension of about 1.2 [25,26]. It is not, however, clear to us why the fractal dimension of the fracture surfaces we analyzed is not simply one less than that of the three-dimensional fracture network. We now propose a model that can explain all these data.
62
M. Sahimi et al. I Earthquake
statistics and fault patterns
4. A discrete model of fracture
What is the theoretical explanation for these results? At small length scales, rock, although microscopically disordered, is macroscopically homogeneous. The nucleation and propagation of fractures in rock do not take place at random, but depend on the stress field in the system. Once a crack nucleates in rock, stress enhancement at its tip is larger than at any other point of the medium and, therefore, the next microcrack almost surely develops at its tip. Hence, at small length scales the microcracks are all connected and cluster together with almost no dead-end microcracks. This is consistent with the experimental observations [19,20,27] that fractures often seem to be clustered together at different parts of rock, separated by large intervals of unfractured rock. However, at larger length scales rock is macroscopically heterogeneous with spatially varying properties. In fact, it has been shown that [1,28] the permeability and other properties of heterogeneous rock are spatially and fractally distributed quantities. The growth of fractures in a heterogeneous medium is very different from that in a homogeneous medium described above. Crack growth continues only as long as a much stronger region of the rock is not encountered. But because of the spatially varying properties of rock, the growing crack will eventually encounter a strong region of the rock and its growth stops, in which case another microcrack in another region of the rock nucleates and starts to grow. The growth of the new crack also stops when it encounters another strong region of the rock, and so on. This has the consequence that, viewed from a large scale, the nucleation and growth of cracks in different regions of the rock take place essentially at random, since they are strongly dependent on the spatial distribution of the rock properties. Therefore, the fractal dimension of the fracture network at large length scales should be the same as that of percolation clusters, since in this process microcracks are created at random (i.e., bonds of a network are broken at random). Indeed, our 2D experimental result for large length scales, D = 1.9, is precisely the same as that of 2D percolation clusters, D = 91148 = 1.9, and our 3D experimental results with larger and larger samples indicate clearly a trend towards D = 2.5, the fractal dimension of 3D percolation clusters. To further test our ideas, we employed our discrete model of brittle fracture developed recently [29-311 to simulate fracture of a macroscopically heterogeneous medium. It has already been established that fracture in microscopically disordered, macroscopically homogeneous media can be successfully simulated with this model. Indeed, in two dimensions this model yields [29,31,32] fractal fracture patterns. If we analyze the statistics of all cracks of such patterns, we find that they have a fractal dimension, D 2: 1.7, in complete agreement with
M. Sahimi et al. I Earthquake statisticsand fault patterns
63
the experimental results. If, on the other hand, we analyze only the fractal behavior of the single fracture that spans the system, we find that it has a fractal dimension of about 1.2, consistent with the results of Okubo and Aki [25] and Aviles et al. [26] for a single major fault. To simulate fracture in a heterogeneous medium, we consider a 90 X 90 triangular network with periodic boundary condition in one direction. Every site of the network is characterized sites are by the displacement vector ui = (uiX, uiY), and nearest-neighbor connected by springs. We consider here the case of brittle fracture for which a linear approximation is valid up to a threshold (defined below). The nodal displacement ui is computed by minimizing the elastic energy E with respect to ui, where
E =
’ 2 2 (ij)
[(u; - ~j)
.R,]$, +
’
2
ZZ(Sqik)‘gijgik . (jik)
Here LYand /? are the central and bond-bending or angle-changing force constants, respectively, R, is a unit vector from site i to site j, g, is the elastic constant of the spring between i and j, and ( jik) indicates that the sum is over all triplets in which the bonds j-i and i-k form an angle whose vertex is at i. Although one may need a more complex form of the elastic energy 1321 for describing all properties of fracture networks, we believe that eq. (2) is adequate for studying the fructul properties of the system and, in particular, its fractal dimension, since D is often independent of microscopic features of E. To make the network macroscopically heterogeneous, we divide it into 25 blocks of size 18 x 18. Within each block g, was set to be the same for all springs, but for each block g,j was selected from the distribution Pi(g) = where 0 < y1 < 1. This is a broad distribution and, thus, the elastic (1 - r,)P, moduli and other properties of the entire network are also spatially dependent quantities and are broadly distributed, as required [1,28]. We then introduce a threshold value I, for the length of a spring, which is selected according to the distribution, P,(I,) = (1 - y,)lLy2, with 0 < y2 < 1. We found that our results are insensitive to the values of y1 and K. We then initiate the failure process by applying a fixed external strain to the network in a given direction, calculating the nodal displacements ui, and breaking the spring for which l- I, is maximum, where 1 is the current length of the spring in the strained network. After a spring is broken, we recalculate the nodal displacements ui for the new configuration of the network, select the next spring to break, and so on. If the external strain is not large enough to break any new spring, we gradually increase it. The simulation continues until the network finally becomes macroscopically disconnected and a macroscopic fracture network is formed. We then
64
M. Sahimi et al, I Earthquake
I
OO Fig. 3. Number of fracture.
of fractures
statistics and fault patterns
I
1
N vs. the radius
I
3
J 4
I&
r. obtained
with the two-dimensional
discrete
model
calculate the fractal dimension of the fracture network by plotting log N vs. log r, where N is the number of broken springs and r is the radius of the growing crack. The results are shown in fig. 3 (where yI = -yZ= 0.6, and /3/a = 0.04), from which we obtain D 2: 1.85, in very good agreement with our results discussed above, which also strongly supports our view that at large length scales the fracture network in heterogeneous rocks should have the structure of a percolation cluster. The fracture network cannot have the structure of a percolation cluster above p, because, as soon as the fracture network becomes sample-spanning, large scale stress concentrations are relieved and the stress distribution becomes localized. As a result, large scale fracture growth can no longer occur. The fact that the model does reproduce the data has two implications. The first is that the measured fractal dimensions are not merely effective values representing possible crossovers between small, intermediate and large length scales. These fractal dimensions do have clear physical interpretations, as discussed above. The second is that, in a highly heterogeneous medium, the generally accepted idea that stress enhancement at the tip of a growing crack always causes it to grow further may no longer be valid. Thus, crack nucleation and growth may become independent of the stress field in a strongly heterogeneous system. 5. Analysis
of seismic data
We now analyze seismic data and interpret the results. Four seismic data sets from four different regions in Southern California were selected for our
M. Sahimi et al. I Earthquake
65
statistics and fault patterns
Table I The parameters of the four sets of data used in the analysis, where m represents the magnitude of the earthquakes. Region
Volume (km)
Number of events
m
Time span
Location error (km)
SA-EL Parkfield Whittier Uoland
100 x 10 x 17 x 18 x
1859 717 224 291
>2 >O.Ol >I >l
1981-1990 1988-1990 10/1987-1211987 2/1990-6/1990
4-5 0.1 l-2 l-2
125 x 20 18 x 15 18 x 15 18 x 15
analysis. Table I compares the parameters for each data set. Each event represents an earthquake of size m. For each set, the fractal dimension was estimated 1331 by counting the minimum number of cubes N of edge length L required to completely cover the data set at various length scales L. For each set, the minimum L was taken to be the location-error (in kilometers) of the set. Fig. 4 presents the results, and it is clear that in all cases the spatial distribution of hypocenters is fractal over at least one order of magnitude, with a fractal dimension of about 1.8. What is the interpretation of this value of D? Geologic features are rarely simple planes, and in fact most faults are comprised of many shear fractures which are visible at a number of scales. Therefore, given our results for fracture networks, it is clear that, at large length scales, the fault patterns should make a percolation cluster with a fractal dimension, D = 2.5 in 3D. The locations of earthquake hypocenters are of course on the network of faults. However, these locations have to belong to that part of the network that is active, i.e. the part through which stress and deformation are transmitted. That is, the locations of earthquake hypocenters have to belong to the backbone of the fault network. Indeed, our result, D = 1.8 is perfectly consistent with the fractal dimension D, of the backbone of 1000
1
Fig. 4. Variations of the minimum number of occupied cubes N of edge length L to completely cover the data set at various length scales L for all four sets of data. The slopes of all data sets are about =1.80.
66
M. Sahimi et al. I Earthquake
statistics and fault patterns
percolation clusters [34], which is also about 1.8. This is the first time that the statistics of earthquakes has been given a clear geometrical interpretation. We remark that Kagan and Knopoff [35] and Andrews [36] had already pointed out the relevance of fractal concepts to the distribution of earthquakes (see below). The fractal dimension of the distribution of hypocentral locations appears to be universal, at least for the four sets of data that we analyzed. Kagan and Knopoff [35] studied the spatial distribution of earthquakes and found, D 2: 2, which they interpreted as meaning that earthquakes occur on planes. But they also found that if they explore larger depths and areas, the estimated fractal dimension decreases. On the other hand, if the hypocentral locations are plotted in 3D space, and rotated around all three axes, there does not seem to be any one plane about which the earthquakes cluster. The San AndreasElsinore (SA-EL) earthquakes are located over a very wide zone of deformation, between and inclusive of two major faults, whereas the Parkfield earthquakes are located over a relatively thin zone of deformation along the SA fault, yet both have a fractal dimension of about 1.8. This universality is also consistent with the fact that the fractal dimensions of the backbone of percolation clusters is also universal, and lends further support to our interpretation of the data. 3D
6. Discussion and conclusions We conclude the paper by making a few remarks. The number N, of earthquakes of magnitude greater than m is related to the energy release E, through, dN,ldE, - rnm7, where 1< r < 2. It has been argued that [37] 7 is determined by the geometric structure of the fault network. In view of our results, it is likely that T is related to the fractal dimension D, of the backbone of the fault network or percolation clusters. It has also been argued [4,37] that the occurrence of earthquakes and their spatial distribution is closely related to the propagation of fractures and the structure of their network in the rock, and our results establish this connection firmly. Kagan and Knopoff [35,37] also proposed a model of earthquake occurrence which appears to provide relatively accurate predictions. In their model it is assumed that each earthquake generates additional shocks with a probability that depends on time as tFH. The analysis of the data showed that a value of 19close to 1.8 seems to simulate the statistical properties of intermediate and deep shocks, and it is interesting to note that 8 = D, = 1.8, the same as the fractal dimension of the network of hypocentral locations. As mentioned above, several authors [3-51 have attempted to explain earthquakes based on the theory of SOCP. However, the
M. Sahimi et al. I Earthquake statisticsand fault patterns
67
application of this theory to the interpretation of earthquake statistics is appropriate if the system is macroscopically homogeneous whereas, as our results indicate, it may be more appropriate to develop the theory for a macroscopically heterogeneous system: the analysis of fracture data discussed above indicates that the fractal dimensions of the fracture networks in homogeneous and heterogeneous rocks can be quite different. Simulation of earthquake statistics in a macroscopically homogeneous system and the interpretation of the results based on SOCP yielded [3] r = 1.35. It may be that the same simulation in a heterogeneous system could yield an exponent 7 closely related to 4~1.8. Finally, in a previous paper [38] we proposed that the universality classes of fracture of disordered media can be characterized by universal fixed points at the fracture point. We suggested that there are two such fixed points, one of which corresponds to systems in which fracture depends on the stress field, while the other one is for systems in which fracture is independent of the stress field. Our results in the present paper indicate clearly that these fixed points correspond to fracture in microscopically disordered and macroscopically disordered media, respectively.
Acknowledgements
This work was supported in part by the U.S. Department of Energy, the Petroleum Research Fund, administered by the American Chemical Society, the USC Center for the Study of Fractured Reservoirs, and the San Diego Supercomputer Center. This paper was prepared while one of us (M.S.) was visiting the HLRZ Supercomputer Center at KFA Julich, as an Alexander von Humboldt Foundation Research Fellow. He would like to thank the Center for warm hospitality, Hans Herrmann for useful discussions, and the Foundation for financial support.
References [l] M. Sahimi, Rev. Mod. Phys., to be published. [2] S.R. Brown, C.H. Scholz and J.B. Rundle, Geophys. Res. Lett. 18 (1991) 215. [3] P. Bak, C. Tang and K. Wiesenfeld, Phys. Rev. Len. 59 (1987) 381; P. Bak and C. Tang, J. Geophys. Res. 94 (1989) 15635. [4] A. Sornette and D. Sornette, Europhys. Lett. 9 (1989) 197. [5] K. Ito and M. Matsuzaki, J. Geophys. Res. 95 (1990) 6853. [6] J.M. Carlson and J.S. Langer, Phys. Rev. Lett. 62 (1989) 2632; Phys. Rev. A 40 (1989) 6470. [7] H. Nakanishi, Phys. Rev. A 43 (1991) 6613. [8] Z. Olami, H.J.S. Feder and K. Christensen, Phys. Rev. Lett. 68 (1992) 1244.
68
M. Sahimi et al. I Earthquake statistics and fault patterns
[9] D. Stauffer and A. Aharony, Introduction to Percolation Theory, 2nd ed. (Taylor and Francis, London, 1992). [lo] M. Sahimi, M.C. Robertson and C.G. Sammis, Phys. Rev. Lett., to be published. [ll] T.L. Chelidze, Phys. Earth Planet Int. 28 (1982) 93. [12] R. Englman, Y. Gur and Z. Jaeger, J. Appl. Mech. 50 (1983) 707. [13] J.C.S. Long and P.A. Witherspoon, J. Geophys. Res. 90 (B4) (1985) 3087. [14] E. Charlaix, E. Guyon and S. Roux, Trans. Porous Media 2 (1987) 31. [15] V. Ambegaokar, B.I. Halperin and J.S. Langer, Phys. Rev. B 4 (1971) 2612. [16] A.J. Katz and A.H. Thompson, J. Geophys. Res. 92 (1987) 599. [17] J.C.S. Long, K. Karasaki, A. Davey, J. Peterson, M. Landsfeld, J. Kemeny and S. Martel, Int. J. Rock. Mech. Min. Sci. Geomech. Abstr. 28 (1991) 121. [18] S. Kirkpatrick, C.D. Gelatt and M.P. Vecchi, Science 220 (1983) 671. [19] C.C. Barton and P.A. Hsieh, Physical and Hydrological-Flow Properties of Fractures, Guidebook T385 (American Geophysical Union, Las Vegas, NV, 1989). [20] C.G. Sammis, L. An and I. Ershaghi, unpublished. [21] D. Sornette, in: Spontaneous Formation of Space-Time Structures and Criticality, T. Ristc, ed., NATO AS1 (Kluwer, Dordrecht, 1991). [22] T. Hirata, Pure Appl. Geophys. 131 (1989) 157. [23] R.-S. Wu and K. Aki, Pure Appl. Geophys. 123 (1985) 806. [24] H. Sato, Pure Appl. Geophys. 128 (1988) 43. [25] P.G. Okubo and K. Aki, J. Geophys. Res. 92 (1987) 345. [26] C.A. Aviles, C.H. Scholz and J. Boatwright, J. Geophys. Res. 92 (1987) 331. in: Reservoir Assessment of the Geysers [27] R.P. Thomas, R.H. Chapman and H. Dykstra, Geothermal Field, Publication TR27 (California Department of Conservation, 1981). [28] T.A. Hewett and R.A. Behrens, SPE Formation Evaluation 5 (1990) 217. [29] M. Sahimi and J.D. Goddard, Phys. Rev. B 33 (1986) 7848. [30] S. Arbabi and M. Sahimi, Phys. Rev. B 41 (1990) 772. [31] M. Sahimi and S. Arbabi, Phys. Rev. B, to be published. [32] L. de Arcangeli, A. Hansen, H.J. Herrmann and S. Roux, Phys. Rev. B 40 (1989) 877. [33] C.G. Sammis and R. Biegel, Pure Appl. Geophys. 131 (1989) 255. [34] M. Sahimi, J. Phys. A 17 (1984) 3073. [35] Y.Y. Kagan and L. Knopoff, J. Geophys. Res. 86 (1981) 2853. [36] D.J. Andrews, J. Geophys. Res. 86 (1981) 10821. (371 Y.Y. Kagan and L. Knopoff, Science 236 (1987) 1563. [38] M. Sahimi and S. Arbabi, Phys. Rev. Lett. 68 (1992) 608.