A fractal model for crustal deformation

Tectonophysics, 132 (1986) 261-269 Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

261

A fractal model for crustal deformation D.L. TURCOmE Department of Geological Sciences, Cornell University, Ithaca, NY 14853-1504 (U.S.A.) (Received August 26,1985; accepted February 20,1986)

Abstract Turcotte, D.L., 1986. A fractal model for crustal deformation. In: B. Johnson and A.W. Sally (Editors), Intraplate Deformation: Characteristics, Processes and Causes. Tectonophysics, 132: 261-269. We hypothesize that crustal deformation occurs on a scale-invariant matrix of faults. For simplicity we consider a two-dimensional pattern of hexagons on which strike-slip faulting occurs. The behavior of the system is controlled by a single parameter, the fractal dimension. Deformation occurs on all scales of faults. The fractal dimension determines the fraction of the total displacement that occurs on the first-order or primary faults. The value of the fractal dimension can be obtained from the f~quen~-mag~tude relation for earthquakes. Our results are applied to the San Andreas fault system in central California. Earthquake studies give D = 1.90. We associate the main strand of the San Andreas fault with the primary faults of our fractal system. We predict that the relative velocity across the main strand is 2.93 cm/yr. The remainder of the relative velocity of 5.5 cm/yr between the Pacific and North American plates occurs on higher-order faults. The predicted value is in reasonably good agreement with the value 3.39kO.29 cm/yr obtained from geological studies.

Introduction

It has long been recognized that many geological processes are extremely complex. However, much progress has been made in the last 20 yrs in simplifying and understanding a variety of these processes. Plate tectonics is one example. Although the hypothesis of rigid plates with well-defined bounda~es is only an approbation, results based on plate tectonics explain many geological processes. Another example of a successful hypothesis in geology is lithospheric flexure. Analytical solutions for flexure explain the basic morphology of ocean trenches, of island and seamount chains and many sedimentary basins. In this approach to geological problems boundary value problems are formulated, and solutions are obtained, either analytical or numerical. In principal any geological problem can be solved by forward modeling. However, the actual solution for a complex fault system is impossible ~~~951/86/$03.50

@ 1986 Elsevier Science Publishers B.V.

to obtain. As a specific example consider the San Andreas system in California. In order to formulate a boundary-value problem all the faults in the system would have to be specified, clearly this is impossible. Solutions can be obtained using finite-element techniques if only the largest faults are included and if other assumptions are made (Bird and Baumgardner, 1984; Turcotte et al., 1984). However, the validity of the results are open to serious question. An alternative approach to complex geological problems is to take a statistical approach. Many geological problems are scale invariant over a wide range of scales. The concept of fractals has been introduced in order to quantify a variety of scale invariant physical processes (M~delbrot, 1982). The original basis for the concept of fractals (fractional dimension) was a rocky coastline (Mandelbrot, 1967). One definition of a fractal is: p _ 11-D

0)

262

where P is the length of the coastline using a measuring rod of length 1. If P as a function of I satisfies eqn. (l), then D is the fractal dimension of the coastline. As the length of the measuring rod approaches zero, the length of the coastline approaches infinity. A typical value For D is 1.2. The concept of fractals can also be used to generate synthetic topography that looks remarkably realistic (Mandelbrot, 1982). An alternative definition of a fractal is: n-A-D/2

(2)

where n is the number of islands with area greater than A; on a world-wide basis it is found that I) = 1.3 (M~delbrot, 1975). It is recognized, both from the distribution of earthquakes and surface geology, that displacements on a major fault zone are distributed over a region. It is possib1.e to use the fractal concept to better underst~d this ~st~bution. As a specific example consider the San Andreas fault in central California. This fault is recognized as a major boundary between the Pacific and North American plates. The relative velocity between the two plates in the region of this fault is 5.5 cm/yr. An important question concerning the San Andreas fault is the fraction of this relative velocity that is accommodated on the fault. The mean velocity on the fault can be obtained in a variety of ways. Creep measurements in central California give an instantaneous velocity. However, it is not clear that the present velocity is a measure of the velocity on geological time scales. Offsets across the fault of dated geological units give estimates of mean displacements over longer times. These observations consistently give velocities that are less than the relative plate velocities. The relative velocity across the fault is probably in the range 3-3.5 cm/yr (Sieh and Jahns 1984). It is not surprising that only a fraction of the relative plate velocity is a~o~odat~ on the San Andreas fault. The region adjacent to the fault has extensive tectonics and many secondary faults. A fraction of the relative velocity between the plates may be accommodated in the rest of the western United States. In this paper it is hypothesized that the zone adjacent to the San Andreas fault is made up of a

fractal dist~bution of faults on all scales. Thus the primary fault, the San Andreas, accommodates only a fraction of the displacement that is accommodated on the system. The largest earthquakes occur on the San Andreas fault although smaller earthquakes also occur. The largest, or characteristic, earthquakes on secondary faults are hypothesized to have a fractal or scale-inv~ant distribution. In many geographic regions earthquakes have a well-defined statistical distribution. The number of earthquakes II’ with magnitude greater than YYE is related to the magnitude by: log N= -b??Zia (3) where a and b are constants. The b-value varies from region to region but generally lies in the range 0.5
(4) where p is the shear modulus, A the area of the fault break, and S the mean displacement on the fault break. The moment of an earthquake can be related to its magnitude by: log M=cpn+d

(5) where c and d are also constants. Kanamori and Anderson (1975) have established a theoretical basis for taking c = 1.5. These authors have also shown that it is a good approximation to take: M= w3

(6)

where Y= A1j2 is the linear ~mension of an earthquake. Combi~ng eqns. (31, (5), and (6) gives: log~=-~log~+lo~~

(71

where: log p = y

t a- :

log cy

And eqn. (7) becomes:

Compa~son with eqn. (2) gives:

as the fractal dimension of the earthquakes in a region.

263

Formulation of the model

We will now formulate a model for the continuum deformation of the crust. The formulation is based upon the work of King (1983). This author applies the concept of fractals to triple junctions where three strike-slip faults intersect. He shows that a fractal distribution of faults yields the frequency magnitude relation given in eqn. (3). In this paper the concept of deformation associated with a fractal matrix of triple junctions is extended to obtain the distribution of deformation on faults of various orders. We hypothesize a uniform matrix of vertical strike-slip faults. All deformations are in the horizontal plane. The basic fault pattern is a symmetric triple junction between Nn = 3 first-order, strike slip faults of length h, as illustrated in Fig. la. The depth of each first-order fault is also taken to be h,. This triple junction forms the interior of a hexagon as shown; the hexagon is our basic building block. The relative velocity across each of the first-order faults is ui. If all the faults are right lateral (or left lateral), the displacements are instantaneously compatible but will result in a deformation at the triple junction that will be discussed in a later section. Assuming complete symmetry, i.e. an infinite set of similar hexagons, the distance between parallel faults is wi = 3’/*h,. Neglecting complex two-dimensional interactions, the earthquake displacement 6, on a first-order fault is related to the stress drop uO by: (11) For this simple model the time interval TVbetween first-order earthquakes is: (12) And there will be N,, = 3 earthquakes in this interval. The total moment of the three earthquakes is: M, = 3ph;S, = 33’2uoh;

03)

In this idealized model the three earthquakes on the first-order faults are assumed to occur simulta-

(b)

Cd) Fig. 1. Illustration

of the fractal

a. Three first-order, form a symmetric each

strike-slip

model for crustal faults of length

triple junction.

The relative

fault is ct. b. Twelve second-order,

length

velocity

strike-slip

There are actually

h,/2

third-order,

strike-slip

of second-order h, - h,/4

faults of length

form a series of symmetric

triple junctions.

four layers of faults, each with a depth number

of third-order fault

is

+/4.

third-order

faults

equilateral

triangles

of

triple

two layers of faults. each with a

so that the total number

24. The relative velocity across each fault is 13/2.

each

across

faults

h 2 = h, /2 (solid lines) form a series of symmetric

junctions. depth

deformation.

hi (solid lines)

h1/4

faults is

c. Forty-eight (solid lines)

There are actually so that the total

faults is 192. The relative velocity across d. The

is illustrated.

pattern All

with side lengths

of faults

first,

second

lie on

and

a set of

h, /4.

neously in order to maintain geometrical compatibility. These results are summarized in Table 1. We hypothesize scale invariance and assume that there is a second-order set of hexagons with side lengths h, = h,/2. The set of second-order hexagons is illustrated in Fig. lb. The sides of second-order hexagons that lie entirely within the primary hexagon are attributed to the primary cell while those on the boundary are divided equally between the primary cell we are studying and adjacent primary cells. Thus there are 12 faults of length h,. However, in order to maintain scale invariance we assume that the second-order hexagons have a depth h, so that there are two layers of second-order hexagons. Thus there are Nn = 24 faults each with an area h$ = h:/4. The assumption that h, = h,/2 appears to be

264

TABLE 1 Summary of fault properties at orders 1, 2, 3, and II

arbitrary. However, the results are basically independent of this choice. Alternative assumptions such as h, = hi/3 or h, = h,/4 gives similar results. In the earth there will be a continuous distribution of fault sizes. By making the choice h, = hi/2 we have determined the bin size into which the continuous distribution is divided. The relative velocity across each of the secondorder faults is v/2 where v1 is the relative velocity across the first-order hexagon that is attributed to displacements on the second-order faults. Each fault is taken to be right lateral. Following the same procedure as for the first-order faults, the earthquake displacement on a second-order fault is: 3r/=a h

3’j2 u h

&=L!&Ty

faults is equal to the total moment on the firstorder faults if v2 = ~1~. The set of third-order hexagons is illustrated in Fig. lc. There are 48 faults of length h, = h,/4 associated with each layer of the primary hexagon. For the four layers there are N, = 192 faults each with an area h: = h:/f6. The velocity across each of the third-order faults is v3/4 where us is the relative velocity across the first-order hexagon that can be attributed to displacements on third-order faults. Again all displacements are right lateral. Following the same procedure, the earthquake displacements on the third-order faults are: 31’2a h

31i2 aoh,

s3w--.u=T_

P

08)

P

The interval between third-order earthquakes is:

lJ

where the stress drop a0 is assumed to be independent of order. The interval r2 between secondorder earthquakes on a given fault is:

262

31/z -e&1

earthquakes

in the

PV2

The number of second-order the period pi is given by:

earthquakes

Ne2 in

(16) The total moment of the earthquakes on the second-order faults in the period or is: M2 - 24phi6,

The number of third-order period 7r is given by:

05)

72=-= v2

(19)

:

= 33/2aoh: ;

(17)

As expected the total moment on the second-order

and the total moment of the earthquakes on the third-order faults in the period or is: M3 = 192ph&

;

= 3’/“uoh; ;

(21)

It is seen in Fig. Id that the first, second and third-order faults lie on a set of equilateral triangles with side lengths h,/4. The above results are easily generalized to the n th-order. There are N,, = 3 X 23(n-1) faults each

265

with an area hi = ht/2 *ln- l). The velocity across each of the n th-order faults is u,/2(“-‘) where u, is the relative velocity across the first-order hexagon that can be attributed to displacements on nth-order faults. The earthquake displacements on n th-order faults are: a,=

3’/*a h

0 “=__

3’/* u,h,

2”_’

P

(22)

CI

The interval between n th-order earthquakes is: (23) The number of nth-order faults is given by: Nfn = 3 X 23(n-1)

(24)

The number of n th-order earthquakes in the period T, is given by: NC, = N,,

:

= 3 x

23(n-1) :

and the total moment of the earthquakes nth-order faults in the period TVis:

on the

All these results are summarized in Table 1. So far the velocity ratio u,,/u, has not been specified. In order to do this we introduce the fractal concept. We define the fractal relation for the number of earthquakes of order n to be: N en N,,

(2”Pl)D=

2D(n-1)

(26)

We will show that this relation is equivalent to eqn. (2). Comparison of eqn. (24) with eqn. (26) requires that: _0, = 2-(3-DKn-1)

(27)

1’1

The total velocity u across the primary hexagon is the sum of the velocities at all orders and is given by: u=

E n-1

““=

VI

E

2-(3-D)(n-l)=

1 _

2:(3_DI

n-1

(28) Combining eqns. (27) and (28) gives: :

= (1 _ 2-‘3-D’)2-‘3-D”“-”

(29)

The fraction of the total velocity u that is accommodated on the primary (first-order) fault is: Ul -

=

1 _

2-(3-D’

U

(30)

If D = 2, then ut/u = i and one half of the total velocity on the fault system is on the primary faults. Solutions can be obtained 0 -ZD -c 3. As D approaches 3 virtually all the displacement occurs on small faults. Applications

As a particular example we will consider deformation associated with the San Andreas fault system in central California. The b-value for earthquakes in central California has been obtained by Shi and Bolt (1982), taking their mean value b = 0.95 we obtain D = 1.90 from eqn. (10) using c = 1.5. We assume that the full relative velocity between the Pacific and North American plates occurs on this system and take u = 5.5 cm/yr. We associate the main branch of the San Andreas fault with our first-order faults. From eqn. (30) we find that the predicted relative velocity across the primary fault is ur = 2.93 cm/yr. The remainder of the relative motion is taken up on higher-order faults. The average slip-rate on the San Andreas fault at Wallace Creek on the Carrizo Plain of central California has been determined by Sieh and Jahns (1984). For the past 3700 yrs they give a value of 3.39 f 0.29 cm/yr which is in reasonably good agreement with our predicted value. Since the trend of the San Andreas fault in Central California is very close to the predicted small circle of rotation required for a strike-slip fault, it is unlikely that a substantial fraction of the relative motion is accommodated elsewhere in the western United States (Turcotte, 1977). We argue that the higher-order faults that accommodate the remainder of the relative velocity between the Pacific and North American plates are adjacent to main strand of the San Andreas fault. Since our model assumes scale invariance it is not appropriate to consider great earthquakes on the San Andreas fault. The rupture lengths associated with these earthquakes are much greater than

1.37

1 1 1 2.93

Cumulative number of faults. Ncrn

Number of earthquakes, N,,

Cumulative number of earthquakes N,,,,

Total velocity, on (cm/yr)

64

88 5.5 x lo=

4.4 x 10’5 6.43

Earthquake interval rn (yr)

Earthquake moment, M,. (dyne cm)

Earthquake magnitude. me,,

5.83

0.685

2.93 41

Velocity on each fault urn (cm/yr)

5

5.23

6.88 x 10”

188

0.160

0.639

19

14

8

1

Number of faults, Nt,,

30

2.15

4

60

120 73

5.5

11

Fault dimension, h, (km)

Earthquake displacement, S, (cm)

3

9

2

1

Order:

4.62

8.59 x 10”

403

0.0372

0.298

71

52

585

512

15

1.375

4

Properties of faults of order 1 to 10 for the San Andreas fault system in central California

TABLE 2

4.02

1.07x10”

863

0.00869

0.139

266

195

4,681

4,096

1.5

0.6875

5

3.41

1.34x IO”’

1,847

0.00203

0.0649

993

727

37,449

32,768

3.75

0.3438

6

2.81

1.68X10””

3,964

4.73x10

0.0303

3,704

2.711

299,593

262,144

1.875

0.1719

7

‘l

2.21

2.10 x 10”

8.523

l.lox1o-4

0.0141

13,792

10,088

2,396,745

2.097.152

0.9375

0.0859

8

1.61

2.62 x 10’”

18,241

2.57x10

0.00657

51,502

37,710

19J73.961

16,777,216

0.4688

0.0429

9

’

1.Ol

3.28 x 10”

39,067

6.00x 10mh

0.00307

192,361

140,859

153,391,689

134,217,728

0.2344

0.0215

10

267

the rupture depths. Thus the first-order faults would have a different geometry than the higherorder faults. For this reason and because extensive information is available we will apply our model to the June 27, 1966 Parkfield earthquake. This earthquake occurred on the primary strand of the San Andreas fault so that it corresponds to an earthquake on a first-order fault in our model. Trifunac and Udwadia (1974) have studied the Parkfield earthquake in considerable detail and give a mean fault displacement 6, = 120 cm, a moment J4i = 4.4 X 1O25 dyne cm, and a shear modulus ,u = 3 X 1O1’dyne cm*. From eqn. (4) the fault area is A, = 122 km2 so that h, = 11 km. Applying our idealized model we find from eqn. (11) that the corresponding stress drop is a0 = 20 bars which is quite a reasonable value. Taking S, = 120 cm and or = 3.39 cm/yr the time interval between earthquakes from eqn. (12) is TV= 35.4 yrs. Bakun and McEvilly (1984) argue that previous earthquakes on the Parkfield section of the San Andreas fault occurred in 1857, 1881, 1901, 1922,1934. The interval between 1934 and 1966 is 32 yrs which is in quite good agreement with the interval 35.4 yrs given above, clearly the earlier intervals were considerably shorter. Using eqn. (5) with c = 1.5 and d = - 16 we find m = 6.43, Magnitude estimates for the 1966 earthquake range from 5.3 to 6.5 (mb = 5.3 USGS, mL = 5.5 Berkeley, mL = 5.8 Pasadena, m, = 6.5 Palisades). As an example of the distribution of seismicity on faults of various orders, we consider the Parkfield example in Table 2. The fault dimensions h, and displacements S, are scaled by factors of 2 according to our basic hypothesis. Assuming a single primary or first-order fault, the number of higher-order faults N,, scales by factors of 8. The cumulative number of faults with a linear dimension greater than h, Ncfn, is also tabulated and is given as a function of h in Fig. 2. To a good approximation NC,- h3. It is seen that the number of small faults is very large. The total velocity on faults of order n, u,, as obtained from eqn. (29) with D = 1.9 is tabulated. Also given is the velocity on each fault un, = u,/N,,. The interval between primary (first-order) earthquakes is obtained from eqn. (12) using S, = 210 cm and u1 = 2.93 cm/yr, with the result 71 =

N cf

IO-*

Id h, km

Fig. 2. The predicted number of faults NcI that make up the San Andreas fault system in central California with a linear dimension greater than h.

41 yrs. This differs from the 35.4 yrs given above because the observed ui = 3.39 cm/yr was used in that case. The intervals between earthquakes on higher-order faults r,, is tabulated and is given as a function of h in Fig. 3. To a good approximation 7 - h-l,‘. The interval between earthquakes on small faults is very long. The number of earthquakes on higher-order faults in the period 7, N,,, is tabulated. The moment of each individual earthquake on an n th-order fault M,, is given as well as the corresponding magnitude men. The cumulative number of earthquakes N,, in the period r with a magnitude greater than m is given as a function of 172in Fig. 4. Good agreement with eqn. (3) is obtained with b = 0.95 as expected. It should be emphasized that our classification of faults is an approximation to a continuous

IQ-’

I

IO

h, km

Fig. 3. The interval between earthquakes 7 as a function of the linear dimension of the fault h for the San Andreas fault system in central California.

268

IO _-II_.

‘0

I

2

3

4

5

6-7

Fig. 5. Illustration

m

Fig. 4. The predicted California

section

time interval

number

right lateral

of earthquakes

of the San Andreas

of void space (dilatancy)

for

N,, on the central

fault system during

r1 = 41 yrs with a magnitude

of the creation

strain with yxJ = 0.05.

the

m or larger.

by: * = &X”

(32)

distribution of sizes. It is in essence a bin classification. The fractal, dimension specifies the fraction of the total displacement that occurs on faults with a characteristic size greater than h. Once the fractal dimension has been specified the number

The dilatant volume can be filled with exolved minerals such as quartz or it can result in sag ponds at the surface.

of faults with a characteristic size greater than h is specified as shown in Fig. 4. Also the number of earthquakes and the mean interval between earthquakes are specified as shown in Figs. 2 and 3.

We have hypothesized that crustal deformation occurs on a scale-invariant matrix of faults. For

As discussed above the deformation of the hexagons is instantaneously compatible for infinitesimal deformations. However, for finite deformations

the corners

of the hexagons

become

offset. This can result either in the removal of material or in the creation of void space. For shallow levels in the crust the creation of void space may result in dilatancy. The creation of void space is illustrated in Fig. 5 for right lateral strain with yXY= 0.05. A triangular void with an area S*/2 is created at each corner of the hexagons. Thus the volume dilatancy is: * =

W2) _ 31/z

2

6 *

3112 ( W 1

(31)

yW2

This result is valid at all orders of deformation so that the dilatancy can be related to the total strain

Conclusions

simplicity we have assumed a two-dimensional pattern of hexagons on which strike-slip faulting occurs. Each higher-order set of faults has a characteristic dimension that is reduced by a factor of two. The crust of the earth is undoubtedly permeated with a much more random distribution of faults. However, experience with other fractal systems, such as porous media, indicate

that our scale

invariant matrix will be a good approximation truly random media.

to a

Having prescribed a relatively simple geometry and scale invariance, the behavior of the system is controlled by a single parameter, the fractal dimension. The fractal dimension determines the fraction of the total displacement that occurs on the primary faults. The value of the fractal dimension can be obtained from the frequency-magnitude relation for earthquakes. It is concluded that only a fraction of the displacement on the San Andreas fault system occurs on the main strand of the fault. The predicted velocity of 2.93 cm/yr is in reasonably good agreement with the value 3.39 + 0.29 cm/yr obtained from geological studies.

269

Acknowledgements

This research has been supported by the National Aeronautics and Space Administration under contract NAS 5-27340 and grant NAG 5-319. This is contribution 826 of the Department of Geological Sciences, Cornell University. The author would like to acknowledge a very helpful review from Toni Gangi. References Aki, K., 1981. A probabilistic synthesis of precursory phenomena. In: D.W. Simpson and P.G. Richards (Editors), Eartbquake Prediction. American Geophysical Union, Washington, D.C., pp. 556-574. Bakun, W.H. and McEviBy, T.V., 1984. Recurrence models and Parkfield, California, earthquakes, J. Geophys. Res., 89: 3051-3058. Bird, P. and Baumgardner, J., 1984. Fault friction, regional stress, and crust-mantle coupling in southern California from finite element models. J. Geophys. Res., 89: 1932-1944. Kanamori, H. and Anderson, D.L., 1975. Theoretical basis of some empirical relations in seismology. Bull. Seismol. Sot. Am., 65: 1073-1096.

King, G., 1983. The accommodation of large strains in the upper lithosphere of the earth and other solids by self-similar fault systems: the geometrical origin of b-value. Pure Appl. Geophys., 121: 761-815. Mandelbrot, B., 1967. How long is the coast of Britain? Statistical self-similarity and fractional dimension. Science, 156: 636-638. Mandelbrot, B.B., 1975. Stochastic models for the Earth’s relief, the shape and the fractal dimension of the coastlines, and the number-area rule for islands. Proc. Natl. Acad. Sci. U.S.A., 72: 3825-3828. Mandelbrot, B., 1982.‘The Fractal Geometry of Nature, Freeman, San Francisco, Calif. Shi, Y. and Bolt, B.A., 1982. The standard error of the magnitude-frequency b-value. Bull. Seismol. Sot. Am., 72: 1677-1687. Sieh, K.E. and Jahns, R.H., 1984. Holocene activity of the San Andreas fault at Wallace Creek, California. Geol. Sot. Am. Bull., 95: 883-896. Trifunac, M.D. and Udwadia, F.E., 1974. Parkfield, California, earthquake of June 27, 1966: A three-dimensional moving dislocation. Bull. Seismol. Sot. Am., 64: 511-533. Turcotte, D.L., 1977. Stress accumulation and release on the San Andreas fault. Pure Appl. Geophys., 115: 413-427. Turcotte, D.L., Liu, J.Y. and Kulhawy, F.H., 1984. The role of an intracmstal asthenosphere on the behavior of major strike-slip faults. J. Geophys. Res., 89: 5801-5816.