University of Southern California Dissertations and Theses

Simulating nonplanar fault surfaces using stochastic branching

(USC Thesis Other)

Simulating nonplanar fault surfaces using stochastic branching

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content INFORMATION TO USERS This manuscript has been reproduced from the microfilm master. U M I films the text directly from the original or copy submitted. Thus, some thesis and dissertation copies are in typewriter face, while others may be from any type of computer printer. The quality of this reproduction is dependent upon the quality of the copy submitted. Broken or indistinct print, colored or poor quality illustrations and photographs, print bleedthrough, substandard margins, and improper alignment can adversely affect reproduction. In the unlikely event that the author did not send U M I a complete manuscript and there are missing pages, these will be noted. Also, if unauthorized copyright material had to be removed, a note will indicate the deletion. Oversize materials (e.g., maps, drawings, charts) are reproduced by sectioning the original, beginning at the upper left-hand comer and continuing from left to right in equal sections with small overlaps. Photographs included in the original manuscript have been reproduced xerographically in this copy. Higher quality 6" x 9” black and white photographic prints are available for any photographs or illustrations appearing in this copy for an additional charge. Contact U M I directly to order. ProQuest Information and Learning 300 North Zeeb Road. Ann Arbor, Ml 48106-1346 USA 800-521-0600 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. SIMULATING NONPLANAR FAULT SURFACES USING STOCHASTIC BRANCHING Copyright 2001 by Eric Jason Libicki A Thesis Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree MASTER OF SCIENCE (GEOLOGICAL SCIENCES) August 2001 Eric Jason Libicki Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. UMI Number: 1409597 UMI’ UMI Microform 1409597 Copyright 2002 by ProQuest Information and Learning Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. ProQuest Information and Learning Company 300 North Zeeb Road P.O. Box 1346 Ann Arbor, Ml 48106-1346 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. UNIVERSITY OF SOUTHERN CALIFORNIA T H E G R ADUATE SCHO O L U N IV E R S IT Y RARK LOS A N G ELES. C A L IF O R N IA SOOO? This thesis, written by Eric Jason Libicki under the direction of h ii Thesis Committee, and approved by a ll its members, has been pre sented to and accepted by the Dean of The Graduate School, in partial fulfillm ent of the requirements fo r the degree of Master of Science D ate A u ^ u s t_ 7 ,.J 0 0 L S S J L / Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Acknowledgments I would like to thank my advisor Yehuda Ben-Zion for pointing out the opportunity of this project to me, helping me along the way, and for his ability to come up with unexpected ways to tackle a problem. Thanks to Shoshana Levin for providing a program used within this thesis, for always being there to answer my never ending questions, and for continuously giving me insights into the world of fractals. Thanks to Zhigang Peng for teaching me every strange UNIX command I needed to know without which I never would have been able to finish my programs. Also thanks to Susan Owen, Charlie Sammis, Mike Tamada, Yan Kagan, David Okaya, and Cindy Waite who each had some part along the way in making this thesis a reality. Special thanks to Mom, Dad, Lisa, Susan, and Cuddles, who always showed support, and were always there for me. And finally, thanks to everybody else who deserves to be acknowledged here, but may have been accidentally left out. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. iii Table of Contents Page Acknowledgments ii List of Figures and Tables iv Abstract v 1 Introduction 1 2 Fractals 4 2.1 Fractal Analysis 5 2.2 Artifacts and Potential Errors 11 2.3 Fractal Studies of Earthquake Fault Geometry 16 3 Probability Distributions 23 3.1 Probability Terminology 23 3.2 Uniform Distributions 26 3.3 Power-Law Distributions 27 3.4 Cauchy Distributions 28 3.5 Poisson Distributions 35 3.6 Power Law Earthquake Statistics 37 4 Stochastic Branching Theory 39 4.1 Branching Terminology 40 4.2 Fault Branching Model and Assumptions 40 4.3 Criticality 44 4.4 Branching Studies of Earthquake Fault Geometry 48 5 Methods 50 6 Results and Discussion 58 7 Conclusion 8 1 Bibliography 83 Appendix A: Program Poipro18.f90 87 Appendix B: Program Slices5.nb 95 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. iv List of Figures and Tables Page Figures Figure 1: Generating a Koch Snowflake 7 Figure 2: An example of a probability mass function 25 Figure 3: An example of a cumulative distribution function 25 Figure 4: PDFs and CDFs for various x~ 2 distributions 29 Figure 5: A power-law and Cauchy PDF comparison with a= .005 & 1.0 32 Figure 6: A power-law and Cauchy CDF comparison with a= .0015 33 Figure 7: Comparing Cauchy and arctan Cauchy RNGs with a= 1.0 36 Figure 8: Comparing Cauchy and arctan Cauchy RNGs with a= 0.1 36 Figure 9: Graphical representation of mathematical tree definitions 4 1 Figure 10: A survival function for = 10,000 mathematical trees 47 Figure 11: Diagrams of the process that generates the synthetic fault 5 1 Figure 12: Fractal dimension plots when using box counting methods 62 Figure 13: A comparison of various slices of a synthetic fault surface 66 Figure 14: A comparison of fault slices with extreme parameter values 7 1 Figure 15: Time progression of a synthetic fault surface 73 Figure 16: Plots showing the change in fractal dimension over time 75 Tables Table 1: General data for synthetic fault surfaces 1 1 through 22 59 Table 2: Fractal dimensions of slices parallel to the y-z plane 60 Table 3: Fractal dimensions of slices parallel to the x-y plane 64 Table 4: General data for synthetic fault surfaces 31 through 42 67 Table 5: General data for synthetic fault surfaces 5 1 through 62 68 Table 6 : Fractal dimensions of slices of different size synthetic faults 69 Table 7: Fractal dimensions using a 3-D box counting technique 78 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. V Abstract This thesis builds off a model presented in the paper "Stochastic Model of Earthquake Fault Geometry” (Kagan, 1982) in order to create and study synthetic fault surfaces to better understand various fractal techniques. The idea of revisiting the work was motivated by the large increase in papers that claim to observe fractal features over the past two decades and for the ability to expand on the model due to advancements in theory and computer technology. It also allows for a reexamination of both the results of the model and its theoretical foundations. Kagan’s model provides an opposite end member to the standard Euclidean-continuum framework and as such it is important to understand the properties and consequences of the model. The model makes use of stochastic branching and probability theory to take a non-deterministic approach to simulating fault surfaces. Since fractals are characterized by power laws, the model incorporates various power-law probability distributions. Using power laws for determining fault geometry complements other observed power-law statistics associated with earthquakes such as the Gutenberg-Richter relationship and Omori’s law. After the surfaces are created, 2-D and 3-D fractal box-counting methods are applied to different aspects of the faults. The results suggest that stochastic branching fault surfaces are more complicated and quite distinct from the mathematical objects used to develop fractal theory and thus do not share all of the same generalizations. This implies that inferences such as those made about the fractal dimension of an entire fault surface from limited data such as fault traces may be incorrect. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 1 1 Introduction This study attempts to model earthquake fault geometry using fractals, probability distributions, and branching theory to create synthetic ruptures. The goal is to have a model that simulates spatial earthquake statistics that have been observed using a fractal framework. There have been many works which measure fault traces and hypocenter distributions using fractal dimensions, but little work has been done to simulate a fault “surface” which has the suggested geometrical properties. This is done in the present paper by building on the previous works of Kagan (1982) and Vere-Jones (1976, 1977). Earthquake modeling serves several functions. Models are used for seismic risk analysis, for gaining a better understanding of rupture propagation, and for trying to predict future earthquakes. Using random elements in the fractal framework model makes full deterministic predictions of the propagation of a fault surface impossible (Kagan, 1997). However, having random elements does not mean that there are not preferred orientations or locations for the fault surface so it can still be used as a tool for seismic risk and for understanding rupture growth. An unusual aspect of fractal models is the lack of clear definitions for concepts related to a single earthquake or fault. When simple fault planes are abandoned in favor of fractal faulting it is difficult to determine friction, surface area, slip, rupture length, and many other concepts associated with planar surfaces (Kagan, 1982). All models involve making underlying assumptions, and earthquake modeling studies often make a fundamental assumption that faults are Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 2 essentially planar. This assumption makes it possible to represent a complicated phenomenon in a very simple way. For example, sliding blocks could be used to represent plates, and any associated equations derived from the model are often less complicated. However, there have been many spatial earthquake observations that do not support the view that faults are planar. So even though a planar model may be useful for approximating stress and other quantities associated with earthquakes, they may not be able to properly resemble spatial characteristics such as fault traces and hypocenter distributions. Although moving from Euclidean geometry towards fractal geometry often adds significantly to the complexity it may be unavoidable in order to accurately model spatial characteristics of earthquakes. The purpose of this thesis is to expand on Kagan’s (1982) attempt to model spatial properties of earthquakes by using a stochastic branching model that utilizes power-law probability distributions which are associated with fractal geometry. Kagan (1982) argued that the stochastic branching model does simulate observed patterns by qualitatively comparing simulation results with “traces of cracks in sections of fracture specimens, with zones of shear fracture in mines [and] with traces of earthquake faults on tectonic maps.” At the time fractal concepts were just beginning to be applied to faults, but since then the topic has been researched extensively. This paper will try to compare calculated fractal dimensions of the synthetic faults to observed fractal data to see if there is a correlation. Additionally, this paper will examine if the fractal dimensions of the model have spatial or temporal changes. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 3 As with ideal models, the model tries to minimize the degrees of freedom because the more degrees of freedom that are involved, the more that non-uniqueness becomes a problem. This paper will also test the non uniqueness by observing how changes in output correspond to varying different assumptions. Trying to have a simple model with minimal parameters is difficult when trying to model the complexities of fractal behavior. There will be an attempt to justify all parameter choices, though when using probability there is not always an obvious physical interpretation associated with a specific distribution. This model does not attempt to focus on trying to fit one particular data set. The objective is to work with a fixed set of parameters and try to match the characteristics of the model to those in several sets of field and lab results. The model will create an irregular surface that shares many properties to those expected of a fault surface. Such a model allows the surface to be completely observed unlike with real fault surfaces. The model thus lends itself to more accurate and advanced measurements that cannot typically be performed with real data. Therefore several fractal measurements will be applied to the model with the goal of being able to use the results to better understand the relationship between various fractal methods and how well such methods work with natural objects. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 4 2 Fractals Many recent papers attempt to describe the geometrical patterns of surface traces and hypocenter distributions. In order to discuss spatial relationships a suitable framework is necessary. Euclidean geometry is often used for describing fault geometry because it provides simple equations for representing shapes such as planes, which are often how fault surfaces are idealized. Although Euclidean models can have a random component added to them to create a more irregular surface they will only be rough over a small scale and still be relatively planar (Power and Tullis, 1991). However, fractal geometry which has been used to study earthquake observations for the past two decades is able to generate irregularities that repeat over multiple scales. Thus they have been used extensively for describing earthquake fault patterns that are suggestive of having a nonplanar geometry. Fractals are ideal for measuring the non-linearity and roughness of mathematical objects such as Cantor dust, Koch curves, and Sierpinski objects (e.g. arrowheads, gaskets, and sponges). However, measuring the fractal dimension of objects found in nature, such as faults, is not as straightforward. The difficulties of this transition are due to the underlying assumptions of fractal analysis. As will be explained in the following section, while these assumptions can be perfectly satisfied by mathematical objects, they are not completely satisfied for natural ones. Unfortunately, papers from different scientific disciplines (e.g. Xu et al., 1993; Malcai et al., 1997) tend to agree that there is currently no definitive method to overcome the difficulties that arise when using fractal measurements with natural objects. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 5 The following background section defines fractals, reviews reasons why fractal analysis can cause errors or artifacts when applied to earthquake fault geometry and examines how these problems can affect results. The purpose of doing this is to examine whether fractals actually characterize earthquake fault geometry, and if so, provide a numerical range of fractal dimensions that have been observed to characterize earthquakes. These values can then be checked against values obtained by the model. 2.1 Fractal Analyaia A fractal can have slightly different meanings depending on its context, but in general it is a self-similar (equivalent to scale invariant! or self-affine object that has an underlying power-law equation and is characterized by a positive fractal dimension. Self-similar mathematical objects retain a pattern that is identical at an infinite number of scales making magnifications of the shape indistinguishable from the original shape. Self-similarity implies that a function f(x) is proportional to f(kx) for all k (where k is some constant). This leads directly to the fact that fractals have an underlying power law since only power-law functions fit this requirement. Self-similarity for natural objects is not the same since they can only u be made superimposable in a statistical sense” rather than geometrically superimposable (Goltz, 1997). This is an important difference that will be addressed later on. Self-similarity also implies that the scaling relationship is independent of direction unlike self affinity where different axes must be scaled by different amounts to have the object appear the same under magnification. Whereas a Euclidean dimension (also called topological dimension) is always an integer, such as D = 1 for a line, and D = 2 for a plane, a fractal dimension is used for Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 6 describing a shape that does not quite fit into such classifications and thus is rarely an integer. Additionally, unlike lines and planes, fractal shapes have a “rough” surface meaning that they are not differentiable. One of the complications of using fractal analysis is that fractal dimensions can be calculated in an infinite number of ways. However, the four most commonly used fractal dimensions are typically Ds, Do, D i, and D2 (using the notation found in Goltz (1997)). Ds is known as the similarity dimension and is defined as log b/log a where b is the number of smaller objects of length 1/a which are necessary to replace a larger object of the same shape. For example, since a 2x2x2 cube can be made up of eight 1x1x1 cubes, b = 8 (since there are eight smaller cubes) and a = 2 (since each small cube has half the length of the original cube). Hence Ds for a cube is log8/log2 = 3 which is expected since a cube fills 3-D space. All fractal dimension calculations for a cube will equal three. Determining the similarity dimension of a self-similar fractal object such as a Koch snowflake shown in Figure 1 is analogous. Note that after each iteration of the Koch snowflake every segment of unit length is replaced by 4 smaller lengths each 1/3 the size of the original. Consequently, b = 4 and a = 3 so the dimension for the perimeter of the Koch snowflake shown in Figure 1 is Iog4/log3 = 1.26. Since 1.26 falls between the values of 1 and 2, the perimeter will fill up more space than a line, but will not fill an entire plane. It is important to note that the similarity dimension can only be used for shapes that are perfectly self-similar (i.e. not only statistically self-similar) over all scales. Therefore similarity dimension cannot be used for natural objects. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced with permission o f th e copyright owner. Further reproduction prohibited without permission. A Length = L Length = L/3 Figure 1: The top panel (based on page 42 of Mandelbrot (1983)) shows a Koch snowflake after two Iterations. The perimeter of the object has a fractal dimension of log4Aog 3. This is demonstrated by the lower panel. Note that after each iteration, every side of the snowflake with length L is replaced by four segments each of length L/3. 8 Three of the most common fractal dimensions used for natural objects are Do. D i, and D2 which are the capacity dimension, information dimension, and correlation dimension respectively. Additionally, there is a corresponding dimension Dn for any integer n. The capacity dimension is essentially equivalent to the dimension derived from box-counting methods such as those used by Robertson et al. (1995). Box-counting techniques take a set of data and attempt to cover (enclose) the entire data set using a minimal number of adjacent, equal-sized, boxes. The data set may be points (e.g. hypocenters), lines (e.g. fault traces), or higher dimensional shapes (e.g. fault surfaces). The “boxes” used to cover the data are often squares to cover 2-0 data sets, or cubes to cover 3-0 data sets. To determine Do a data set is covered multiple times using varying box sizes. In terms of box-counting, Do is equal to the limit l- * o log(N)/log(1/L), where N is the number of boxes (which contain at least one piece of data) used to cover the set, and L is the length of the boxes used. A graph with log(N) plotted on the y-axis, and log(1 /L) plotted on the x-axis is used to find Do as L approaches zero. The resulting graph should theoretically plot as a straight line with the fractal dimension equal to the slope of the line. If log(L) is plotted instead of log(1/L) on the x-axis as is often done, Do will be equal to -1 times the slope. There are a few final points about Do worth mentioning. Note that by rewriting the equation for Do as D0*log(1/L) = log(N) - ► log(1/L)Do = log(N) _ ► (1/L)Do = N a power-law relation is obtained as is expected with fractals. Another Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 9 important detail is that a box which covers only one data point is treated the same as boxes which contain multiple data points, therefore the capacity dimension cannot be used to differentiate uniform from non-uniform point distributions. The information dimension takes the distribution of points into account by giving more weight to boxes that contain more points. Calculating D i is similar to computing Do except the numerator uses the average information provided by the covering rather than the number of filled boxes. The information i provided by an event E is defined as i(E ) = -log P(E). P(E), the probability that the event E occurs, is calculated by determining the chance that a random data point can be found in a particular box. Goltz (1997) notes that information can be thought of as the amount of “surprise.” When the probability of an event is certain (P(E) = 1) or impossible (P(E) = 0), the amount of surprise is zero (i(E ) = 0) or infinite (i(E ) = °°) respectively. Do is equivalent to D i when using a uniform distribution of points, otherwise D i is less than Do- Although comparing Do to D i will show if the distribution of a set of points is non-uniform, correlation dimension can also make this distinction and is often the preferred technique when working with earthquake fault geometry. The correlation dimension uses the correlation between pairs of points rather than box-counting techniques. Determining the correlation function C(r) is a multistep process. First, the distance between every two points in a data set is calculated. Then those distances are compared to a given radius r. The number of distances which are less than rare tallied and divided by Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 10 N2, the number of points squared. D2 is then found by the expression lim it r-*o logC(r)/log(r). Since there are (N)*(N-1)/2 distances to calculate for data sets with N points, having many points may pose a computational problem in terms of the time needed to perform the calculation. One method used to resolve this problem is to choose a group of selected reference points and only calculate distances which involve at least one of these points. This will change the final correlation function, both by reducing the number of distances that are compared to r and by dividing the tally by N*Nre f (where Nr e f is the number of reference points) instead of N2. However, using only a selected number of points to determine the dimension is only an estimation and therefore may be less accurate. All these different dimensional calculations are used together because an object is characterized by its set of fractal dimensions {D -c o ,.... D o..... Deo}. Objects that have all fractal dimensions equal to each other are monofractals. Otherwise the object is multifractal and has Dn s Dm for n > m (e.g. D2 s D i s Do). Multifractality is caused when the amount of clustering varies. It is currently unclear whether fault geometry is monofractal or multifractal, but due to faults being natural, as opposed to synthesized by mathematical equations, some authors favor a multifractal approach to data. However, Eneva (1994) points out that determining between the two is not as simple as theory suggests since synthetic monofractals tested by sampling a finite number of points appear multifractal. This uncertainty is unfortunate because distinguishing whether earthquake fault geometry falls into the category of Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 11 monofractals or multifractals can lead to different inferences. Some of the most significant differences between monofractals and multifractals are the inferences that may be made about the fractal dimension of the subsurface. For example the equation Dsubsurface = Dtrace +1 as mentioned in Okubo and Aki (1987) is only relevant to certain homogeneous fractals, thus only a monofractal scenario would allow surface results to be applied easily to the subsurface. 2.2 Artifacts and Potential Errors Miscalculations and artifacts can easily arise when using fractal dimensions. There are several general problems to consider that are not specific to any particular dimension. When measuring fractal dimensions for data sets it is important to clearly define the scaling region since it is impossible for a finite set to be self-similar at all scales. When a log-log plot is used to find the dimension, the plot will not be linear over all scales and thus upper and lower bounds of the self-similarity need to be defined. If this range is too small self-similarity may not be valid. Goltz (1997) suggests that spanning at least one order of magnitude is necessary for validity whereas Gonzato et al. (1998) proposes spanning “at least two or three orders of magnitude.” Spanning too large a range may also be a problem if it does not make sense physically. For example, extending an argument in favor of self similarity for faulting should not have an upper range much larger than the size of the brittle lithosphere, nor a lower range below that of grain size. Most of the literature tends to show self-similarity only over a small range of magnitudes. Avnir et al. (1998) took a sample of 96 papers and found all of them to report self-similarity for at most three orders of magnitude with the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 12 average being only about 1.3 orders of magnitude. Since fractal theory was designed to apply over an unlimited number of scales it is necessary to understand the reasons that most fractal studies are limited to a few orders of magnitude. One common reason for the scaling to be restricted is data resolution. Malcai et al. (1997) have also found that in certain cases the self similarity of natural systems can at most be tested for two to three orders of magnitude since fractal characteristics may deteriorate over time. Although they are referring to c lu s te rin g systems, such an argument could apply to the case of fault surfaces. Determining the range of self-similarity is difficult as Eneva (1994) points out that there has not been a standard way in the literature for determining how to choose it. A simple method is to pick a segment of the log-log plot that looks linear whereas a more thorough analysis may include making calculations which examine the change in slope between pairs of data points (Eneva, 1994). However, there is no agreed upon standard for distinguishing between acceptable and unacceptable deviations. There are other problems that also arise when trying to fit a line to log- log fractal dimension plots. One such concern is trying to fit a line to too few data points. Gonzato et al. (1998) suggests that at least twenty points are needed to test the robustness of a linear fit. They explain that having a lack of data points can arise when the box-counting method is performed by hand, rather than through the use of a computer, since not enough different box sizes are used. Even with enough data points there is still the potential problem of multiscaling (Goltz, 1997). It may be the case that the data are best fit with two or more lines of different slopes. This would suggest that the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 13 fractal dimension being sought after actually varies over different ranges. This, for example, could happen in faulting if different physical processes are dominant at different length scales. Lei et al. (1993) suggests that preexisting, nonfractal length scales such as crustal thickness may be threshold values for where fractal dimension can vary. If one tries to fit one line through multiscaling data then the dimension would be overestimated for some portions of the scale and underestimated for others. Another problem is the case where the data are anisotropic (Goltz, 1997). For example, hypocenter data appears to more often cluster parallel to the fault “plane” than perpendicular to it. If one tried to examine hypocenter data by looking at a 2-D slice or projection, such as using epicenter data, the apparent resulting dimension may depend on the angle between the Earth’s surface and the dip of the fault. So using 2-D fault measurements to infer the subsurface fractal dimension may not work. This is a problem with the previously mentioned equation “Dsubsurface = Dtrace + 1 ” that is often suggested in papers without being supporting by data. Robertson et al. (1995) shows that it does not even hold for all simple mathematical objects such as a Sierpinski gasket. In terms of using it to understand fault geometry it is even more problematic. Matsumoto et al. (1992) claims that fault traces are not representative of the fault at depth due to near surface heterogeneities. Okubo and Aki (1987) though seem to imply the opposite by summarizing that seismic data is just as complex at depth if not actually more complex. Although as discussed in a subsequent section, relocating seismic events has simplified the apparent geometry at depth. Either way the absence of lithostatic pressure changes the conditions at the surface from that at depth (Kagan, 1991). Therefore it Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 14 appears that geological conditions would tend to support a multifractal environment and the equation presented above could not be used to infer such a simple relationship between Dtrace and DS ubsurface The previous errors and artifacts mentioned were based on the interpretation of data. Other types of general errors caused in fractal calculations involve the data itself. Many papers (e.g. Okubo and Aki, 1987; Hirata, 1989) mention that trying to determine the fractal dimension of faults is only as good as the maps themselves. One potential difficulty with maps is that the geometry is dependent on the resolution. For example, if a map is on the scale of kilometers, the finer details at the meter or centimeter scale will be lost. Therefore two problems arise, first is that testing for scale-invariance at scales that are not featured on the map would be meaningless. Second, the dimension of a trace would be underestimated since the low resolution would smooth out the high resolution detail. A possible way to resolve this issue is to combine maps of different resolutions together. Unfortunately, even with multiple maps there is still the problem of accuracy. Faults are often inferred from features such as scarps, offsets, and sag ponds which are all affected by erosion (Okubo and Aki, 1987). Erosion will destroy the finer details of the fault, once again leading to an underestimation of the dimension. There are many specific considerations to take into account when using box-counting to determine Do. Since box-counting relies on the number of filled boxes, it is necessary to minimize the number of boxes used at every chosen length scale (Robertson et al., 1995). Arbitrarily covering the data with a randomly placed grid will rarely minimize the number of filled Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 15 boxes and thus make the numerator too high, and overestimate Do- A computer program can avoid this issue by attempting to place the grids in numerous positions until finding the optimal placement. Box-counting by hand will usually not be able to find the optimal placement when working with many data points. The finiteness of the data also causes an inherent saturation problem. Once the box length decreases to the point where no box contains more than one point, the capacity dimension reaches saturation since the numerator cannot continue to increase while the denominator still does. This underestimates the value of Do- It is suggested that the box counting method should not be applied once the number of filled boxes reaches more than 20% of the number of data points (Goltz, 1997). Another method that calculates fractal dimension is the ruler counting method (also known as the divider method) which can only be used for line segments. The technique involves trying to approximate a fractal segment using adjoining straight segments (rulers) of a particular length. Classic examples that use this technique involve trying to measure the coastlines of countries (Mandelbrot, 1983). It has also had limited success in measuring geological structures such as the perimeters of lava flows (Turcotte, 1997). The method relies on the observation that the smaller the ruler size used to measure coastlines or perimeters, the greater the apparent length. The slope of the graph of log(total length) vs. log(ruler size) equals 1 minus Do- There are two inherent problems with ruler counting. One problem is that part of the end of the jagged segment will never be covered so there will be a remainder of the segment that is less than a full ruler size. The other problem is that the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 16 segment cannot have gaps or branches thus making it an inadvisable method for using on fault traces which often do have such features. According to Goltz (1997) Do is difficult to determine for values greater than 2 citing that approximately 10°° to 40°° points are necessary for reliability whereas using the correlation dimension needs as little as 1/5 the number of data points for accuracy. However, there are also disadvantages when using D2- A correlation function, unlike box-counting, is altered by edge effects since it is not independent of boundaries. Since an earthquake region is finite, the correlation for points that lie near boundaries will be low. This is due to half the area within radius r of those points being empty since it will span outside of the data’s range. This will lead to an underestimation of D2- Correlation dimension is also sensitive to noise which is inherent to earthquake data. For example, if a hypocenter’s location is given with a possible error of +/-e meters, then determining the correlation for a radius r< e is likely to be significantly corrupted by noise. As just explained, using fractals to measure roughness contains many subtleties. Hence it is always necessary to understand which method a particular author uses when calculating the dimension, and to verify if proper precautions were taken. 2.3 Fractal Studiaa of Earthquake Fault Geometry The difficulty with using fractals to study earthquakes is the large variation in results derived from different techniques such as the ruler and box counting methods. Even using the same method can lead to different results when looking at a particular fault. Choosing to examine only the main Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 17 branch of a fault versus measuring the entire fault system will lead to different fractal values. Even different maps of the same area can vary the calculations for a particular fault. Gonzato et al. (2000) have attempted to quantify the amount that fractal dimensions vary when using 2-D box- counting techniques on digitized images. They have determined that the initial box size and placement both can cause significant errors. If the initial bounding box is larger than the width of the object being measured the dimension will usually be underestimated by a few percent. This problem can be resolved though by starting with a box length equal to the width of the object being measured (Gonzato et al., 2000). However, they found that a more significant problem is that the value of a fractal dimension will change if the boxes are set at a different angle with respect to the object being measured. This problem can cause overestimates of about 6 to 8 percent (Gonzato et al., 2000). Assigning a fractal dimension to an object can be difficult, but monofractals calculations have been applied to several types of earthquake statistics. The ones most important to this study are that of spatial characteristics such as the geometry of fault traces and distribution of hypocenters. Okubo and Aki (1987) and Aviles et al. (1987) independently attempted to characterize the fractal dimension of segments of the San Andreas fault. These two studies alone begin to show the difficulty in using fractals to measure fault traces. Okubo and Aki used a circle covering method over the main trace of the San Andreas and nearby branches to arrive at values of D that range from 1.1 to 1.4. Aviles et al. strictly focused on the main fault line and used the ruler technique to arrive at a value of D from Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 18 1.00 to 1.02. As explained before the ruler estimate should not work well for faulting. They handled the lack of connectivity of the San Andreas Fault by connecting it into one trace and disregarded all branching segments. These factors should tend to underestimate the value of 0 since it treats the San Andreas Fault as a smoother object than it actually is. There is clearly a discrepancy between the values though it must be taken into account that they include different parts of the San Andreas Fault and the results are complicated by the fact that they use different techniques to measure dimension. Therefore, it is hard to compare the studies. Hirata (1989) comments that self-similarity was not taken care of properly in the San Andreas papers since they focused solely around the main fault, and not on the sub-branches of branches. Working on faults in Japan and using the box-counting method, Hirata (1989) arrived at values of D that varied from 1.05 to 1.60 for the southern faults. However, Hirata’s analysis has in turn been critiqued by other authors. One such group is Matsumoto et al. (1992) which obtained values for D from 1.0 to 1.3 for Japan’s fault traces. Although they used the circle covering method rather than box-counting, they argue Hirata’s box-counting method may have produced artifacts since he got values of D < 1 for Japan’s northern faults which they say does not make intuitive sense. However, they also get values of O < 1 for fractal calculations they did of the fault system in the Philippines. They claim this is caused by the Philippine fault system being almost linear, but with varying gaps in the fault. Regardless of the reason for low values of D, Gonzato et al. (1998) explains that another problem with Hirata’s results is that lines are fit to data sets that have only five points and Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 19 therefore have little validity. To summarize, the results so far suggest that D varies widely both between and within fault systems, but it is unclear as to how much of these differences can be explained by the variations in the methods chosen. Hirata (1989) suggests that rock fracturing has an upper limit for D around 1.6 while Kagan (1991) suspects that D may be universal and the discrepancies between various results are due to the many causes of uncertainty in the data, but neither of these assertions currently appear to be sufficiently supported. Unlike that of fault traces, there has not been an extensive amount of work done measuring the fractal dimension of hypocenter data. One of few such studies was done by Robertson et al. (1995) who looked at five data sets in California. They obtained values of Do between 1.79 and 2.07 using the 3-D box counting technique. One important aspect of the results that will be referred to later is that these values were approached asymptotically as the density of events increased over time. All the studies mentioned so far have focused on using one type of fractal dimension to entirely characterize the dimension. As mentioned previously one value of D only characterizes homogeneous fractals where Do * D i * D2 = Deo. Thus it should not theoretically matter which Dn one calculates. However, if natural fractals are multifractal, i.e. heterogeneous fractals, different Dn calculations lead to different results. The greater the difference between the values of Dn, the larger the heterogeneity. Several studies characterizing earthquakes as multifractal processes include Hirabayashi et al. (1992), Lei et al. (1993), and Ouillon et al. (1996). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 20 Lei et al. (1993) looked at epicenter distributions of earthquakes in China and distributions of granite acoustic emissions in laboratory compression testing. For the epicenters they found three distinct scaling regions for D2 of 1.47,1.02, and 1.06 with changes at 80km and 300km length scales respectively. They compared values of D2 to D15 and found them to be significantly different thus concluding that epicenters are multifractal. They used epicentral data because of the error in the vertical direction when determining hypocenter locations. However, epicenters collapse 3-D distributions into 2-D and therefore are an artificially planar distribution that probably cannot be applied to hypocenters. For the granite sample they found two distinct values of D2 (2.1 and 2.5) with the transition occurring at the average grain size. They also found it to be multifractal, but acoustic emission data misses many small events which could alter the fractal dimension. Examining other studies, Hirabayashi et al. (1992) came to the conclusion of multifractality for epicenters of California and hypocenters of Greece and Japan while Ouillon et al. (1996) found multifractality in Saudi Arabia though only for scales greater than 6 km. These wide variations make it difficult to generalize the results. An important issue when measuring the dimension of hypocenters is the reliability of such records in earthquake catalogs. Calculating the fractal dimension of hypocenters only works under the assumption that the hypocenter data are accurate. Without being able to directly observe hypocenters, the determinations for the location and time for hypocenters are not necessarily correct. This is currently a prominent research topic and Rubin et al. (1999) are one of several groups going back to earthquake Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 21 catalogs to relocate the data using waveforms. The basis for relocation is that two earthquakes that occur in the same area should produce almost identical waveforms. By comparing waveforms to each other it is evidently possible to get much better relative arrival times. This leads to a significantly more reliable location of two earthquakes relative to each other. Hypocenters that have appeared to be scattered over a wide area have been collapsed into several discrete streaks and become more planar. This phenomenon has been more recently observed by Waldhauser and Ellsworth (2000) using double-differences on the Hayward Fault. The method makes use of the theory that two hypocenters which are relatively close have similar ray paths to a given station and so any difference between the travel times is due to the distance between sources. These results imply that the fractal dimensions of hypocenters may come into question since the distribution of the hypocenters is clearly affected by noise. While there does not appear to be any literature quantifying the variation for the fractal dimension of hypocenters due to noise, Gonzato et al. (2000) briefly examined how “noise” introduced by a scanner when digitizing maps will vary the calculation. They showed that different scanner resolutions will digitize fault traces with varying thicknesses, and straight lines which are at least 5 pixels in width may cause an overestimation of the dimension by more than 10 percent (Gonzato et al., 2000). Although there have been several papers that have calculated the fractal dimension of fault traces and hypocenters there has been little done trying to model such results. Other than Kagan (1982) one such paper is Turcotte (1986). However, Turcotte uses a fractal model to study deformation Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 22 and the model itself is made up of hexagonal strike-slip faults so it does not try to accurately simulate fault geometry. Additionally, D is assigned a value in the model rather than generated by the model. Turcotte’s paper though is one of several which correlate the fractal dimension of epicenter distributions to the b-value obtained in the Gutenberg-Richter relationship. He uses the equation that D=2*(b-value) to give yet another possible value of D = 1.9 for the San Andreas Fault although this is for epicenters not fault traces. Though the equation is frequently used, other authors say that it only is valid as a first approximation (Davy et al., 1992) or in fact that it does not hold true (Goltz, 1995). The one possible conclusion from all these works is that faulting does not appear planar, but the degree of self-similarity and range of values for D are not well determined at this time. Finding a model to simulate a nonplanar fault surface can therefore be useful as an improvement over models that assume perfectly planar faults, but it will be difficult to test with comparisons to current observations. Thus, this paper attempts on working towards such a model by building off the stochastic branching model introduced by Kagan (1982). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 23 3 Probability Distributions Probability distributions are vital to the model discussed in this thesis since statistical approaches are often useful for modeling complex situations. An introduction to various probability distributions is given in the following sections in order to understand and justify their use for different aspects of the model. After defining probability terminology, there is a brief discussion about discrete and continuous probability distributions. This will set the foundation to elaborate on four specific distributions relevant to the developed model. 3.1 Probability Terminology Random elements are going to be added to the model by using random number generators (RNGs). RNGs are algorithms which randomly select numbers as output based on some probability distribution. In order to properly discuss these probability distributions several terms associated with probability theory need to be defined. In general, probability distributions are associated with random variables usually denoted by X. Despite its name and notation, a random variable is a function not a variable. A random variable X assigns values Xk for every possible outcome k of an event being measured. For example, when working with discrete distributions any value Xk will equal P(k) (0 s P(k) s 1), the probability of k occurring. Since distributions only determine probabilities, not actual values, it is unknown a priori exactly what the actual distribution of numbers from an RNG will be. Therefore it is useful when discussing random number generators to define the expected value (or mean) E[X] of the random variable X to have some approximate idea of what numbers will be generated. By definition Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 24 E[X] = £ n*pn for the discrete case, where pn is defined as the probability n that the number n is generated. As an example, let a probability distribution be defined as po = 1/4 p i = 1/4, p2 = 1/2 , and pn = 0 for all other n. So E[X] = 0*1/4 + 1*1/4 + 2*1/2 = 5/4 in this case. Thus the expected mean of the numbers generated is 5/4, though it is important to note that since only integers are generated in this example, “5/4” itself will never appear as output. It is often tedious or impossible to explicitly write out an entire probability distribution. To resolve this problem there are two commonly used graphs to represent such information for the discrete case. They are the probability mass function (PMF) f(x) and the cumulative distribution function (CDF) F(x). As their names suggest, f(x) is a graph of px for every x and F(x) shows the cumulative sum of all the px up through x. Plots of f(x) and F(x) for the same probability distribution given in the previous example of expectation value are shown in Figures 2 and 3. Continuous distributions are used when a random variable is a continuous function. For this reason continuous distributions have some fundamental differences from discrete distributions. One difference is that the probability that a specific value occurs is always zero because there are an infinite number of possible values. Therefore it is more relevant to determine the probability that a value will occur within a particular range rather than being equal to a specific value. Continuous distributions have a probability density function (PDF) f(x) as opposed to a probability mass function. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 25 f(P x ) 3/4 . Po P i p 2 p 3 Figure 2: A graph o f the p ro b a b ility m ass function fo r {p 0 = 1/4, p, = 1/4, p2 = 1/2}. 3/4 1/4 , P Po Figure 3: A graph o f th e cum ulative d istrib utio n function fo r {p0 = 1 /4 , p, = 1/4, p2 = 1/2}. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 26 However, there is a CDF just as in the discrete case which is given by the integral of the PDF. As x approaches infinity the CDF approaches 1 since the probability of all of the possibilities occurring is always 1 by definition. In other words the probability that a generated value is less than infinity is 1. Most continuous functions do not integrate to one and therefore are not PDF’s, but many functions that converge to some positive value can be changed into a proper PDF by using a normalizing constant. A normalizing constant is a constant multiplied to an entire function so that it integrates to 1 and thus the normalized function can often represent a PDF. Note that a normalizing constant does not change the shape of the function only its amplitude. It is important to realize that the y-axis of a PMF only ranges from zero to one since it represents the frequency of an event. This is not the case for the y-axis of the PDF which can have any value. The values of the y-axis for a PDF taken alone have no meaning since the probability of any particular event is zero. However, the y-axis of the PDF is essential for finding the probability of an event occurring within a certain range of values. For example, the probability that an event occurs in a range [m.n] is the integral from m to n of the PDF. Similarly, the mean of a continuous distribution is not found with a summation as in the discrete case but rather by taking the integral of x*f(x)dx over the entire domain of f(x). 3.2 Uniform Distributions There are four relevant distributions in this work and the uniform distribution is the most straightforward of the four. It is a distribution that has no preferred outcome. Since any possible outcome is equally likely, it would Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 27 be a common choice for situations that deal with homogeneous spaces. The PDF of a continuous uniform distribution is a horizontal line and its CDF is a sloped line from zero to one. The other three distributions that will be needed are the power-law distribution, the Cauchy distribution, and the Poisson distribution which are more complex and are used for more specialized situations. 3.3 Power-Law Distributions Power laws are linked with fractals since self-similarity only occurs with power laws. Power laws are observed in the Gutenberg-Richter relationship and Omori’s Law as will be discussed in detail in the section on power-law earthquake statistics. A power-law distribution has a PDF which can be described by any power law such as x-2 or x~3. Although power law PDFs can have simple formulas they pose a problem in that they cannot be used in their entirety. The problem occurs with small positive values of x. When “x” is near zero, f(x) approaches infinity, or more importantly the integral does not converge and cannot be a probability distribution. This makes it impossible to use without proper modifications. A common solution is to impose a minimum cutoff when using power-law PDFs. This means that there is some minimum number below which values cannot fall. This must occur in the Gutenberg-Richter relationship and Omori’s law, and makes sense for power-law distributions which describe other natural phenomenon, since there must be a finite size and number of earthquakes (even if the minimum cutoff is as small as grain size). The minimum cutoff is an adjustable parameter and it can be difficult to justify what it should be for a particular situation. Another issue that occurs with a minimum cutoff is that it Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 28 also affects the value of the normalizing constant since it changes the lower bound of the integral. Even small changes in the cutoff value can have as significant an influence as the exponent in the power law on the distribution of values (Figure 4). The power-law distribution will be used to describe the variations in angles between parts of the synthetic fault surface that will be created. Unlike the case of a minimum earthquake size, a minimum rotation angle could easily be zero. However, using a power-law distribution without a minimum cutoff cannot be used since an infinite number of the rotations would be zero radians. One possible option would be to make the minimum value so small that it is essentially zero, although as previously discussed the distribution may change considerably depending on the minimum value. Another option is to use a probability distribution that closely resembles a power law for larger values of x (known as having a power-law tail), but is tapered to a finite value as x approaches zero. One such function that can fit these requirements for certain cases is the Cauchy distribution. 3.4 Cauchy Distributions The Cauchy distribution has a PDF of the form (a/n)/(x2+a2) where a is some positive value and the domain of x is all real numbers. Note that the integral equals 1 regardless of the value of a and therefore does not need a normalizing coefficient. A Cauchy distribution forms a bell shaped curve like a Gaussian distribution (also known as a normal distribution), however, there are significant differences between the two curves. They both can produce many numbers near zero when used for RNGs, but the tails of a Gaussian distribution decrease exponentially whereas a Cauchy distribution falls off as Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 29 PDRx) 3 2.5 2 1.5 1 0.5 x 0.25 0.5 0.75 1 1.25 1.5 1.75 2 CDRx) 0.6 0.4 0.2 x 0.25 0.5 0.75 1 1.25 1.5 1.75 2 Figure 4: The PDFs and CDFs for normalized x-2 distributions using three different minimum cutoffs, .001 (solid), .01 (long dash), and .1 (short dash). The minimum cutoff has a significant effect on the distribution. Note that the upper part of the PDFs are not shown. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 30 x ~ 2. This means that a Gaussian curve decreases much more rapidly than the Cauchy distribution. Another way to realize this is by comparing the mean of the curves centered around zero. The symmetry about zero of the bell shaped Gaussian curve causes the expectation value to be zero. Although the Cauchy curve is also symmetric about zero, the mean does not exist for the Cauchy distribution since the tails drop off so slowly (Rice, 1995). This can be verified by noting that the expected value, which equals the indefinite integral of (xAx)*(1/(1+x2))dx when a = 1, does not converge. Kagan (1982) prefers to substitute a Cauchy distribution for a power- law distribution in his model after determining that the best power-law distribution for the model is x ~ 2. The advantages of this is that a minimum cutoff parameter is not needed as in the power-law distribution since the curve tapers off as it approaches zero, and for large values of x the distributions are similar since it has a power-law tail. However, there are still several possible concerns in making such a substitution in distributions. The first is that there is still the decision of a value for the parameter a in the Cauchy distribution. Second, is that a Cauchy distribution has a range of [-co, oo] as opposed to the range of [minimum cutoff, °°l for the power-law distribution. Third, is that although the distributions may be similar for large values of x, they may not be similar for small values of x, which are where most values in the model reside. The first concern is resolved in the setup of the model. The other two issues were not addressed by Kagan (1982) and will be addressed here. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 31 Changing the range of the Cauchy distribution to make it similar to the range of a power-law distribution is not complicated. Due to the symmetry of the Cauchy curve, one can take the curve with x values only from [0 ,« j and use a normalizing coefficient of 2. Since the integral of the Cauchy distribution from [ - < » , 0] = [0,« ] = 1/2, the integral of the new curve will still be one (therefore a proper PDF) and retain the same general shape except for the amplitude. Since the minimum cutoff for the power law distribution is close to zero, the power-law and Cauchy distributions now span similar ranges. The final concern is whether it is possible for the two curves to have similar shapes. The two distributions appear very different if a= 1.0 for the Cauchy distribution as in many examples, and the power-law is x*2 with a minimum cutoff of .001 (the power law has been normalized). However, as the value of a is decreased the Cauchy distribution starts to have a sharper peak, and appears to resemble the power law x-2 (although with small enough values of a it is actually similar to curves with other values in the exponent also) (Figure 5). This is expected since the a in the numerator is just a scaling factor and the a in the denominator essentially becomes zero leaving only x' 2 times a constant. Fortunately, the model only uses small values for a since large output values are only wanted occasionally. It is still useful to check the similarity of the two curves. This can be done by examining the differences between the two CDF curves shown in Figure 6. There is actually a significant percentage (= 37.4%) of values generated by the Cauchy distribution that are less than .001, which is below the minimum cutoff of the power-law in this example. Noting that the two CDF curves never Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 32 PDF(x) 1.4 1.2 1 0.8 0.6 0.4 0.2 1.2 0.8 1.4 0.4 0.6 1 0.2 Figure 5: A comparison of a power-law PDF (solid) and two Cauchy PDFs (dashed). The Cauchy PDFs [0«J have a = 1.0 (short dash) and a = .005 (long dash) while the power-law PDF [.001,°°] is x-2 The curves approximately overlap when a = .0015. Note that only a portion of the PDFs are shown. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 33 CDRx) 1 0.8 0.6 0.4 0.2 - 0.002 0.004 0.006 0.008 0.01 0.012 0.014 Figure 6: A comparison of the cumulative distribution functions of a Cauchy (dashed line, a = .0015) and power-law (solid line, x-2 minimum cutoff of .001) distribution. Note that a significant percentage of numbers generated by the Cauchy distribution are less than .001, but by .0065 the distributions are quite similar. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 34 cross yet always get closer, the power-law distribution generates more values than the Cauchy distribution for any range above .001. When x = .0065 the two CDF curves are within .01 of each other. This means most of the difference between the two functions occur for values less than .0065 radians (~ .37 degrees) which is ideal since the small angles have little effect in the model and therefore the Cauchy approximation is valid. Another issue is that the Cauchy distribution with an unbounded domain is going to be used in the model to generate an amount of rotation which can only range uniquely from [0,2n]. Two questions that arise are how to transform an unbounded distribution into a bounded distribution, and more importantly, will that mapping keep the shape of the Cauchy curve. Although Kagan (1982) explains the transformation geometrically, it is algebraically equivalent to taking the arctangent of a Cauchy random variable to produce a new random variable with a distribution from [0,ti/2]. For small values the mapping will not change the values significantly since the arctan(6) ~ 6. Since most of the rotations will be small, most of the mapping will not alter the distribution. However, it is the large rotations that are the most important in distinguishing the model from a planar one, and those that will significantly be changed by the arctangent function. Kagan (1982) did not address this issue so it will be examined here. Just as with the comparison of Cauchy and power-law functions, it turns out that the value of a affects whether the transformation is reasonable or not. Taking a value of a = 1.0 leads to a distribution which generates many random numbers that will have a value above t t/2. When the arctan of the values is taken to change them into values less than n/2 the resulting Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 35 distribution is quite different in shape leading to a significantly different CDF than the original Cauchy distribution (Figure 7). However, the values chosen for a are much smaller than 1.0 in this model, leading to a very small proportion of randomly generated number to be above n/2 and most are very close to zero. This means that the arctangent will have a large effect on only a small proportion of the numbers. Taking the example of a = 0.1 leads to a much closer similarity between the CDFs (Figure 8). So the transformation should be acceptable for small values of a especially since a will be chosen to be much smaller than 0.1. This leads to the two conclusions that Kagan's Cauchy approximation of a power law and taking a transformation of it using the arctangent function are both acceptable because of the low values chosen for a . 3.5 Poisson Distributions Unlike the three previous distributions discussed, the Poisson distribution is a discrete distribution rather than continuous. It has the PMF (e~^)*(Xk/k!) where k can only be nonnegative integers. A convenient aspect of the Poisson random variable is that the only parameter X is also the expected value and variance of the distribution. Besides simply having useful mathematical properties, it is also a distribution that is used to model a variety of phenomenon ranging from scientific ones such as radioactive emission, to industrial ones such as telephone connection reliability (Leon- Garcia, 1994). These are modeled as Poisson processes because they are stochastic processes (since X is a rate) which satisfy three assumptions. The three assumptions are that the events rarely occur “simultaneously”, that the average rate of the events occurring does not change, and that the events Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 36 PDRx) 0.8 ; °-7 i 0.6 tV " °.S| \. 0 .4 1 ' 0.3 1 0.2 \ 0.1 ; CDRx) 0.8 I 0 .6 i 0.4; 0.2 ^ x X 1 2 3 4 5 6 7 8 8 Figure 7: A comparison of the PDFs and CDFs of a Cauchy RNG (dashed) versus arctangents of numbers generated by a Cauchy RNG (solid) when a s 1.0. Note that the two functions in both graphs are quite different and in fact the latter’s PDF is constant for the special case of «=1.0. The arctangent of a positive number has a range from [0,n/2]. PDRx) CDRx) 6 5 4 3 2 1 0.8 - 0.6 0.4; 0.2 [ ] x x 0.25 0.5 0.75 1 1.25 1.5 1.75 2 0.25 0.5 0.75 1 1.25 1.5 1.75 2 Figure 8: A comparison of the PDFs and CDFs of a Cauchy RNG (dashed) versus arctangents of random numbers generated by a Cauchy RNG (solid) when a=0.1. The functions are similar to each other unlike in Figure 8. Note that the scales are different from the previous figure. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 37 occur independently of each other. By defining “simultaneously” to be two events occurring in the same time interval, the first assumption can often be satisfied by making the time interval significantly small. The second assumption can usually be satisfied for events where the surrounding conditions of the event does not change, therefore keeping the average rate relatively constant over time. The final assumption of independence is a difficult condition to satisfy since it does not apply to all situations. This will be an important factor to consider in the model that is presented. 3.6 Power Law Earthquake Statistics The main reason for considering a nonplanar model is that faults do not appear planar. Just by visual observation, fault traces are not straight lines, and hypocenters do not plot on a single plane. This does not mean that the spatial distributions of earthquakes are fractal; all it suggests is that it may not be planar. It could be that the faulting occurs on several planes, making it look fractal. It can also be argued that faults actually are planar and it is just that surface traces are not representative of the fault at depth due to the free-surface boundary, and that locating hypocenter distributions are subject to noise which makes it appear nonplanar. However, the reason that the spatial distributions are often considered to have the potential to be fractal is because power-law distributions are observed in other types of earthquake statistics. As stated earlier fractals have an underlying power-law distribution, although power-law statistics are not solely generated by fractals. The two most common power laws for earthquakes are the Gutenberg- Richter relationship and Omori’s Law. The Gutenberg-Richter relationship Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 38 states that the log of the number of earthquakes greater than magnitude M decreases linearly as magnitude increases. This proportion can be rewritten in terms of seismic moment as N(M) = A * -5) where Mo is seismic moment, A is a constant, and b is the b-value. The b-value usually ranges from 2/3 to 1 though it appears that it can span anywhere from 1/2 to 2. It also varies between seismic catalogs and over different ranges of magnitudes (Frohlich and Davis, 1993). Omori’s Law shows that the rate of aftershocks decays proportionally to T P where T is time and p usually varies from 1.0 to 1.4 (Lay and Wallace, 1995). Just like the b-value, the p-value has also been claimed to be related to the fractal dimension of faults. Nanjo et al. (1998) found that Do for active fault traces was approximately equal to -2.1p+3.5. As this thesis will show, deducing such relationships from only limited fault information may be incorrect. However, these power-law observations suggests that size, and time distributions of earthquakes are fractal. What remains to be determined is whether the spatial distributions of earthquakes can be characterized as fractal. Results such as Kagan’s (1994) determination that hypocenter distributions are fractal suggest that this may in fact be the case. However, it is still inconclusive as to whether fractals truly characterize earthquake geometry due to problems such as noise errors as mentioned previously. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 39 4 Stochastic Branching Theory Stochastic branching theory is an area of mathematics which is used to study cascading processes that develop over time. The fundamental concept of the theory is that one starts with a single event which has the possibility of generating additional events according to a probability law. These new events (or offspring) may also have their own offspring and the process could continue indefinitely depending on the probability law governing this growth. One of the most common applications of branching theory involves examining the growth of family trees. For this reason branching is often represented as a mathematical tree having familial terms such as parent and generation to describe the relationships of vertices contained within the tree. However, this theory is widely applicable and has been used in many different sciences. For example, fission reactions, random electron emission, and the survival of mutant genes have all been studied with branching models (Karlin and Taylor, 1975). Part of the purpose of this paper is to determine the potential for applying branching theory to earthquake faulting. Although branching has not often been used in geophysics it may be ideal in this case due to the many analogous aspects between branching trees and fault geometry. For instance both may start with one particular event, continue to grow, branch into multiple directions during its growth, and eventually terminate. However, there has not been an extensive amount of research studying fault geometry with this approach. Two proponents of using branching in geophysics are Vere-Jones (1976,1977) and Kagan (1982). Although faulting is caused by many physical processes, branching is essentially a mathematical approach Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 40 to studying phenomenon that treats physical processes as having reached a complexity at which they can be thought of as a random process. Nevertheless, branching under certain probability distributions is still able to generate power-law statistics similar to those observed by earthquakes. It can also be incorporated into either a fractal or multiplanar Euclidean approach to studying earthquake phenomenon. 4.1 Branching Terminology The foundation of a branching model is a mathematical tree. It is thus useful at this point to define terms that are specific to trees. The following definitions are shown pictorially in Figure 9. A tree is a set of vertices (also referred to as points or individuals) which are connected by edges (lines) and contains no loops (edges which form an enclosed region). Each tree begins with an initial vertex known as the root The root vertex is referred to as a parent if it is adjacent to other vertices. All vertices connected to the root are denoted as offspring or children of the root. The offspring can simultaneously be parents by having their own offspring. This is the case if they are adjacent to vertices other than their parent. Following the trend of familial associations, the root is generation zero and each set of offspring is another generation (i.e., the root’s offspring is the first generation). The resulting tree creates a branching pattern in which every vertex is directly connected to at least one other vertex. 4.2 Fault Branching Model and Assumptions Mathematical branching trees are composed of a set of vertices connected by edges, but in order to apply it to a model for faulting the vertices and edges must represent physical quantities. Vere-Jones (1976) uses the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 41 Parent/Offspring Root / Vertices Generation Figure 9: A graphical representation of definitions associated with a mathematical tree. This example has 1 1 vertices, 10 edges, and 4 generations. All trees have only one root. The specific parent highlighted above has two offspring which are shown in a lighter shade of gray. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 42 edges to represent crack lengths while the vertices act as points of either branching or termination. The idea is that increasing the external stress will increase the probability that a given branch will propagate. The reason for using branching in the first place is because it uses a probabilistic view of crack propagation. This assumes that stresses which influence the orientation of the crack and the ability to continue rupturing will vary. However, Vere-Jones suggests that the edges may also be considered to represent an amount of energy released, due to increasing crack length causing a corresponding increase in energy released. In this paper though, as in Kagan (1982), the vertices are going to represent small fault patches and the edges will model various property relations between the patches that form the overall fault surface. There are several implicit and explicit assumptions when using a branching model. Although branching models are often given a time scale to represent the ongoing branching process, it is usually assumed that branching at a given vertex can only take place at one instance in time. Once branching or termination takes place at that point, no further activity will occur there. The time aspect is the stochastic part of this model and is an important factor in that it can be used to treat the extension of the crack as occurring over discrete time intervals rather than as a continuous process. Since all offspring of a given generation are created at the same time tit is possible to define the total number of vertices for any generation as the random variable Xf. The number of branches at each vertex is decided by a probability distribution function. The choice for this function is an important part of the model and will be explained in detail later in the paper. A single vertex Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 43 produces k (from 0 to °°) offspring with probability pk- The number of offspring for each vertex is assumed independent of all others. Note that m a m £ /n = 1 from the definition of probability and that the summation starts at 4=0 k = 0 because each vertex must have a nonnegative number of offspring. Also note that po represents the probability of a branch terminating because no offspring would be generated. Other important assumptions involve the choice of the probability distribution. Vere-Jones (1976) assumes that the number of branches and the lengths of the cracks at any point are statistically independent from any other point, in other words the propagation at one location does not affect the process at other locations. This means all the assigned probabilities are independent which considerably simplifies calculations. Another assumption is that the probability density function is constant throughout space. This does not necessarily imply that the space must be physically homogeneous, but it does mean that the model assumes a statistically homogeneous space. Therefore every vertex is given the same probability distribution function. The advantages of using this statistical branching model compared to a physical model are due to the apparent complexities of the underlying physics. While a crack is propagating there are many interacting processes occurring simultaneously. The orientation of the crack, the strength of the rock, the applied stress, and the heterogeneities within the rock may all influence crack propagation essentially causing elements of the propagation to occur as a random process (Vere-Jones, 1976). Thus studying the lengths and branching of individual segments together as a conglomerate rather than Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 44 individually determining all the physical influences occurring at each segment may be a proper representation that avoids complexities. Before looking at how branching has been used for studying cracks and faults, it is important to have some insight behind the mathematical foundation for branching specifically how one parameter, the average number of offspring, plays a critical role in the theory. 4.3 Criticality An important aspect of the branching model is the mean offspring parameter p (the average number of offspring per vertex).' Vere-Jones assumes this to be constant for every vertex in the tree. Depending on the value of this number, the tree may or may not continue indefinitely. Vere- Jones (1976) notes that since cracks do not propagate indefinitely, a tree which is meant to represent crack propagation should only continue for a finite amount of time. Therefore the mean offspring parameter should not be large enough to cause indefinite branching. There are three general cases to explore, p< 1, //= 1, andp > 1. The case of p < 1 is where each new generation on average would contain less vertices than the previous generation since each vertex on average generates less than one new vertex. Summarizing the proof presented in Lawler (1995) it can be shown that a tree where fi < 1 must terminate after a finite amount of time. To sketch the proof, if one assumes that there are k vertices in the nth-generation of the tree, and each vertex on average has n new offspring, then the next generation should on average 1 p. will be set equal to the value of X from the Poisson process for this paper, but a different symbol is used here for distinguishing the general case from that of a specific distribution. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 45 contain kp vertices. Following this reasoning the nth+2 generation contains on average k^ 2 vertices. Ultimately, if there are k vertices at the first generation, the n1 * 1 generation will on average have k/in_1 vertices. The expression k/tn~ 1 converges to zero when p < 1, and na-oo and therefore no offspring are expected after many generations and the tree terminates. Trees which have p < 1, are referred to as subcritical. In the opposite case ot p > ^ each generation is expected to be larger than the previous one. Therefore, there is a positive probability that these trees will never terminate. These cases are referred to as supercritical, but since they have the possibility of growing catastrophically, they are not used in branching models of fault propagation. A special and unintuitive situation occurs at p = 1. This is referred to as the critical state because it is the threshold between values that cause a tree to terminate and a tree to grow indefinitely. If every vertex had exactly one offspring, the number of vertices at every generation would be exactly the same and the tree would continue indefinitely which would seem to imply that p = 1 has the potential for unlimited growth. However, it is important to differentiate between the case of each vertex having one offspring versus averaging one offspring. The actual number of offspring per vertex was stated to be random, and therefore there is only one case out of infinite number of possibilities that every vertex has exactly one offspring. According to probability theory the chance of this particular case happening is 1A » which is zero. In all other cases the actual average number of offspring will Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 46 not sustain a value greater than or equal to 1 indefinitely.2 Any tree terminates the first time its actual average falls below one. This is because there will be an average of n/(n+1) offspring per vertex when the tree terminates since every vertex is an offspring except for the root. There is a mathematical proof in many sources on branching (e.g. Lawler, 1995) that the average of a random branching process with fi = 1 will at some time fall below one and therefore the tree will terminate. As the parameter // approaches chticality the distribution of total crack length is a power-law function with an exponent of -1/2 (Vere-Jones, 1976). Figure 10 illustrates this power law based on a numerical simulation. If the total length of the branching is equated to total energy as Vere-Jones (1976, 1977) suggests, the Gutenberg-Richter relationship of log(N) = k-bMs can be used to equate tree size to seismic moment to determine the b-value associated with branching.3 Substituting magnitude with seismic moment using the equation log(Mo) = 1 5*Ms+c, the Gutenberg-Richter relationship can be written as N = k*Mo‘ (b/1 -5) where N is now a function of seismic moment. Replacing N with the critical power law relationship lo r1/2 leads to the b-value for this branching model being .75 after setting the exponents equal to each other.4 Therefore the critical state may be ideal for modeling 2 The ‘actual average” is the mean calculated for a specific tree, as opposed to the expected average fi. 3 Ms is the magnitude, N is the number of earthquakes with a magnitude 2 Ms, b is the b-value and k is a constant (the constant k has a different value in each of the subsequent equations). 4 The relationship kx'1/2 jS shown in Figure 10 and is also the value derived by Vere-Jones. Interestingly, Fisher et al. (1997) derived an earthquake size distribution very similar to Vere- Jones (1977) that also led to a b-value of .75. They utilized a very different technique which involved using a 2-0 pl anar fault in a 3-D solid and setting the model’s parameters to critical values. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 47 Survival o f Trees as a Function o f Size 4 > N 5 5 4 > 4 > (_ (n 4 ) 4) L _ 4) . O E 10000 9000 -O .S 2 5 9 y = 1378 2x 8000 R2 = 0 .9 9 9 6 7000 6000 5000 4000 3000 2000 1000 10 20 30 40 Num ber o f V e rtic e s in T ree Figure 10: The survival function for approximately 10,000 trees (of size a 2) created with a random Poisson offspring distribution having a mean of 1.0. The survival function shows the number of trees that are greater than or equal to a particular size. The function above can be equated to the Gutenberg-Richter relationship which uses the number of earthquakes greater than or equal to a particular size. This graph only shows the survival through 50 vertices, but the power-law continues indefinitely. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 48 earthquake faulting. This is an aspect that Vere-Jones regards as promising for the usefulness of branching. However, there are several uncertainties which include the non-uniqueness of obtaining this value and more importantly the wide range of reported b-values as discussed earlier. 4.4 Branching Studies of Earthquake Fault Geometry One point worth noting is a subtle difference between Vere-Jones’ (1976) model and later branching models. Vere-Jones attempts to give a physical interpretation in addition to the geometrical one by assuming there are microcracks throughout the rock, and the crack propagation occurs by the extension of one crack coalescing and triggering other already existing microcracks. This though is just one of several possible mechanisms for fault propagation (Segall & Pollard, 1983; Sharon & Fineberg, 1996; Ben-Zion & Sammis, 2001). Both Vere-Jones’ 1977 work and Kagan’s 1982 work downplay a possible mechanism involved in the branching, but rather focus on the important strengths of a branching model, which are the potential geometrical and temporal similarities to crack propagation. They represent branching from the standard mathematical approach of starting at some initial crack which propagates into multiple cracks rather than coalescing with previously existing ones. The focus of Vere-Jones’ branching model is to examine the distribution of fault sizes generated whereas Kagan gives the model additional spatial characteristics to focus on the geometrical patterns it creates. This difference allows for the branching to also simulate jumps in fault surface traces. Additionally, Kagan uses branching for a fractal framework by using a model that branches thousands of times. However, a Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 49 similar model with significantly less branching could be used for a multiplanar Euclidean framework. It is difficult to assess the validity of applying branching to earthquake faulting. Some aspects in favor of using such a model are its geometric similarities to faulting at the surface, its apparent tendency towards power- law distributions as criticality is approached, and its ability to represent a complex process using probability. However, problems with branching involve the many degrees of freedom that come with changing it into a physical representation for faulting. Even with a justification for setting ^l-^, there are still several parameters such as the length of branches, the angle of branching, and the chosen distribution for the offspring, which need to be considered. This flexibility allows branching to be applied to a fractal framework or a Euclidean model where faulting occurs on several planes, but there is still the difficulty of justifying the values used for each of the parameters in a particular branching model. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 50 5 Methods The main focus of this thesis is to attempt to use probability distributions and stochastic branching theory to simulate nonplanar fault surfaces and test them using a fractal framework. I have created my own code using Fortran 90 for the computational aspects of the project and programmed in Mathematics for the graphical portion. I began by following the methods used by Kagan (1982) to write a program which would simulate nonplanar faults. The following paragraphs describe the process of creating the synthetic fault which is also shown as a series of diagrams in Figure 11. As in Kagan (1982) the basis for simulating nonplanar faults starts with a branching tree. This branching tree will be used to try and imitate the branching observed in fault traces. The tree for this model is created using a critical Poisson distribution. It starts with an initial root and the average number of offspring for each vertex is determined by a Poisson distribution using the single parameter X discussed earlier. Although the parameter X does not necessarily have to be an integer, all vertices will still have an integer number of offspring. Kagan (1982) used values of X slightly less than 1 to avoid the tree from continuing indefinitely. However, as explained previously, using a value of 1 will always create a tree with a finite number of vertices, so X is set to 1 in this paper.5 Kagan (1982) does not explain why he chose a Poisson distribution. However, as explained previously it has the convenient features of generating only nonnegative integers and it can also easily be made “critical" by setting the single parameter A . to 1 as done here. 5 While the majority of the code in this project is original, the code for the Poisson random number generator, along with other random number generators I used, are taken from Numerical Recipes (Press et al., 1996). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 51 a ) 7 7 7 First patch in x-y plane and centered at (0,0,0) b) m n Before Rotation of Offspring: Patches are in the same plane as parent C) PD F(0> » 23 After Rotation of Offspring: Patches are not in the same plane Figure 11: A diagram of the process that generates the synthetic fault, a ) The root (shaded gray) of the tree becom es a circular patch of unit radius in the x-y plane, b ) In this exam ple the root has two offspring (shaded dark gray) so two circular patches are random ly centered som ewhere on the perim eter of their parent's disk and are initially placed in the plane of their parents d is k c ) A probability distribution is used to generate two param eters 0 and < p which rotate the offspring. The process will be repeated until all th e vertices in the tree represent a patch of the fault surface. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 52 Note from the previous Figure 10 that although the number of offspring per vertex are Poisson distributed the total branch lengths have power-law frequency-size statistics. Once a Poisson distributed branching tree is established, the next step is to attempt to transform it into a model for earthquake faulting. Following Kagan (1982), each vertex represents a very small patch of a fault and the patches taken together construct a fault with an irregular surface determined by properties assigned to the edges. Note that a Poisson process is history- independent, whereas fault growth depends on the history. However, a Poisson process is a useful starting point for a statistical treatment of fault growth and may still provide a reasonable framework. My model as in Kagan (1982) has each vertex become a disk of uniform radius. He chooses a specific radius of 40 to 50 meters because the size of each patch is later converted into an amount of slip. Since this paper is focusing on the spatial distribution the actual value of the radius is not as important, thus the circular patches are each given a radius of one unit. However, it is important that the patches are small relative to the overall fault size to minimize discretization effects. Additionally, this helps justify the assumption that all of the patches are the same size and shape. The simulation continues by attaching one patch to another patch for every vertex that is connected to another vertex in the tree. They are connected by placing the center of the disk randomly along the perimeter of an adjoining disk with the placement determined by a uniform distribution. For example, the root of the tree corresponds to a circular patch initially placed in an x-y plane. If the root has two offspring then there are two more Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 53 patches added to the initial patch centered randomly along its perimeter and initially in the same plane as the parent patch. If the patches were never rotated then it would create a set of adjoining disks that all lie in a single plane. However in order to simulate the nonplanar aspect of faults, the disks are individually rotated according to a probability distribution. Kagan initially chose to work with power-law distributions since they are the foundation for fractals. Three power-law probability distributions, x_ 1 -5, x'2 0 an(j x-2.5 were tested in addition to a Cauchy distribution. So in this model, slightly differing from Kagan, the disk is placed in its parent’s plane and two random numbers are generated from the power or Cauchy distribution.6 The arctangent of these numbers are taken to give < t > (the amount of rotation about the disk’s x-axis) and 0 (the amount of rotation about the disk’s z-axis) which each range from -tt/2 (clockwise) to ti/2 (counterclockwise). The positive x-axis is defined here as towards the reader, the positive y-axis is to the right, and the positive z-axis is upward. The patch is then perturbed in both the 4 > and 6directions by the specified amounts. The values of < | > a n d ©are usually small due to the chosen distribution and the new patch typically resides in almost the same plane as its parent. The offspring however will always rotate by some amount if choosing a power distribution due to the minimum 0 and 0 values. Whether 4 >or 0 is rotated first does not change which plane the disk will finally reside in. The reason for the independence of the order of rotation can be understood by visualizing the disk as initially being embedded horizontally in a sphere with a pole that extends from the center of the disk to 6 Kagan (1982) uses a more complex quaternion method for rotations which is not necessary here since no slip is added to the synthetic fault surfaces. Further discussion of the advantages and disadvantages of his method are discussed in Kagan (1990). Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 54 the north pole of the sphere. The first rotation of phi will rotate the embedded disk counterclockwise around the x-axis so that the pole normal to the disk will be at a latitude phi degrees west of the north pole. The second rotation of theta will rotate the disk counterclockwise around the z-axis so that the pole normal to the disk will be at a longitude theta degrees counterclockwise of west. Changing the order of these two rotations will change the path the pole takes to its new position, but the final position of the disk will be the same. The process of choosing < | > a n d 0 is repeated if the offspring has its own offspring and will not terminate until the tree’s final generation. I created the Fortran 90 program Poipro18.f90 (which can be found in Appendix A) to generate a stochastic branching tree where each vertex represents a patch of the surface. Additionally, the program determines five values (x, y, z, 4 >,and 0) for every patch of the fault surface. The values x, y, and z are the three Cartesian coordinates which represent the center of the patch. The two other coordinates 4 >and 0 signify the latitude and longitude rotations, respectively, of each rotated disk. While the output generates five coordinates to describe the orientation and location of the disk, it is actually only three independent coordinates. The x, y, and z values can be uniquely determined by choosing an angle o>(from 0“ to 360°) along the perimeter of a parent’s disk, except for the case of the root disk which is automatically centered at the coordinate (0,0,0). The values for every disk in a tree are read into my Mathematica program Slices5.nb (which can be found in Appendix B) to graphically represent the fault surface. Although the program originally attempted to create every disk simultaneously, it was found to take a significantly long time Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 55 to run and take up a large amount of memory. However, it is not necessary for all the disks to be drawn together because only 2-0 slices of the 3-D fault surface are going to be examined at any given time and viewing it as such also makes it easier to visualize, understand, and interpret the fault surface. In the few cases where the fault will be represented as a 3-D object, the disks will be represented by a single point located at the center of the patch. Thus it is not necessary for Slices5.nb to consider the patches as disks since a 2-D slice of the fault surface, which is the intersection of disks with a plane, consists of a group of lines.7 If they do not intersect it is unnecessary to keep track of that patch for that particular slice which saves computing time and memory. The created surfaces are going to be examined using a fractal framework due to their complex properties. There are many possible ways one could attempt to calculate fractal dimensions of these simulated fault surfaces. For this paper the chosen method is usually to use a 2-D box- counting method on the slices taken of the 3-D fault surfaces. The box- counting method will be performed using a program obtained through anonymous FTP from the site ibogeo.df.unibo.it where it is freely distributed by Gonzato et al. (1998) who wrote the code. The program accepts monochrome bitmap (“.bmp”) files and performs box-counting on a graphic by checking for shaded pixels in a particular box. The program starts by using one large box that surrounds the entire region and continues to shrink the boxsize by factors of until the boxes are 1 pixel by 1 pixel. Although the program can be used to calculate the fractal dimension immediately from 7 It is extremely rare that a disk and plane are tangent, or for the disk to reside entirely within the intersecting plane. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 56 the data it generates, it is not an advisable way to do so since it uses every data point that it collects. Thus it would include points where saturation (as explained in a previous section) occurs and places where the box size is so small that the finite pixel resolution may alter the dimension. Additionally, it appears that the program does not try to minimize the number of boxes filled adding another source for error. In order to have some estimate for the reliability of the program I attempted to find the value of D for objects with known fractal dimension. Gonzato et al. (1998) provided graphical files for a Koch curve (similar to one side of a Koch snowflake with 0 = 1.262) and a Sierpinski carpet (D = 1.585). I have found that the program determines the fractal dimension of the Koch curve to have a D * 1.33 and the Sierpinski carpet to have a value of D * 1.60. Importing the data into a spreadsheet program and fitting lines to subsets of the data shows that a large fraction of the discrepancy comes from the endpoints. The largest box sizes have an inherent saturation problem while at the other end of the scale the smallest box sizes have a resolution problem since details finer than one pixel will not be resolved. The saturation of the large box sizes can occasionally be seen in the data when two different side lengths are tested yet the same number of boxes are needed to cover the data for each case. After testing various possibilities, removing the data points for the two largest and two smallest box sizes seems to be the best compromise and gives a value for the Sierpinski carpet to be D = 1.59 and the value for the Koch curve to be D «1.31 (many of the endpoints need to be removed for the value of the Koch curve to be within .01 of the actual value). Even after correcting for saturation, the value for D will still have some error Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. since there are a finite number of points used for the graphics though the number of pixels is not restricted by the size of the computer screen. However, in terms of time and memory the files can only be so large so they have been chosen to be 1028x1028 pixels and the finest details of the fault surface are lost. Testing the program with various other resolutions appears to cause some variation of the calculated dimension though the values are usually confined to a range of 0.1. Therefore the error of the values given for the following results that use this technique is considered to be approximately +/- .05. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 58 6 Results and Discussion It is hypothesized that the fractal dimension calculated for the 2-D slices and thus the fractal dimension for the entire surface will depend on many spatial and temporal factors, which include any or all of the following i) the initial input parameters, ii) the slice’s distance away from the initial patch at the coordinate (0,0,0), iii) whether the slice is parallel or perpendicular to the initial patch (which tends to govern the direction of the fault due to the high probability of small rotations), iv) the number of patches in the fault surface, and v) the percentage of the fault surface that has finished being generated. There are two input parameters when choosing a power-law distribution for the rotation. They are the exponent for the power-law distribution and the value for the minimum cutoff. The Cauchy distribution only has the single parameter a. Values that are kept constant are the mean offspring for the Poisson distribution and the radius of a patch which are each set to 1.0. Additionally there is a seed for the random number generators, which will cause the tree to be a different size depending on the seed. To begin testing whether the initial input parameters and the distance away from the root centered at (0,0,0) change the values for Do I began by choosing a seed which created a surface with 5,360 patches before terminating. I then had twelve different sets of parameters, nine of which were various power-law distributions and three which were Cauchy distributions. The chosen sets of values are shown in Table 1. Not surprisingly, the different values for the parameters create different values of D as shown in Table 2. This was in fact expected since Kagan’s (1982) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. F a u lt S u r fa c e # d-val u e m in , rotn,/« minx m axx miny m a x y m i n i m a x i x r a n g e ytan g e z rang e 59 C O If ) 5 © C M C O C O 5 © o © C M $ c - ■ M - © © o © © C O © ' ■ M - ’ C O o C O C M d © C O § 5 ) C M ID C O o C O © «• * . C O C O © © C M o C O C O 0 0 © © ' T— d © d d 0 0 © © o © C O © o © If ) © r> - o © r*» © © N . O C O C M C M © lr i C M in C M rs ; C M S i C M 8 c d C M C M 8 © C M © ' C M 8 C O C M C M r> - © 8 © 8 © TT C M C M C O © © £ 8 c *. ' co r- ^ C M o C M d © ' c m ' o © O ) « © o © ■ m ; o C M 8 C M C O © C M © h - o C M o © ' o C O d o C M C O 8 © 8 © © o * - 8 © 8 © If ) © ' if ) © © © © © © C O © ' S C M If ) h . 8 8 9 © 8 8 £ ■M - O © C M N © C M C M C M C M - o C M o o C M O o C O p 8 C M s © © 8 C M 1^ 3 O 2 © © C M C M C M C M 8 © T f C M C M C O C M S i N C M © ' C M C O C M C M © C O C M © co C O C M C M 1 ■M- • 1 • T • © ' t • • • 1 1 in © I ' . C O © C M C O C O % n < o < O < O < o < < o o < O < o < o < o < O < © © © to o ' o ' © o o o ia m in >»>»>. £ £ £ O O U 3 3 3 CO (0 <Q o o o * - CM C O ^ I f ) CO N C O O l O ^ ^ ^ ^ ^ ^ ^ ^ (y » - C M CM CM 0 ) 0 ) £ O ( 0 Q . 8 e o in c « c 8 CM CM C O a > 8 ■ c 3 C O ® c ' O . w ® + — ® 'T' -- £ lc « (0 X e ® B 5 £ Q-.te U rn s ® S T 3 3 <g x : 0 X O C T 5 3 M TO ® ® O - 2 r ; © S o £ © © h . c c o i= E u _ . B • CO _ c © £ o | S l i f f 3.2 ^ « 3 £ . a € r 1 - c © = 2 “ 8 2 o 2 § ® £ - . i ■ J S ® C 2 ~ E 812 g e (5 E © © A £ Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. instead. T h e follo w in g entries co n tain t h e m in im u m a n d m a x im u m values th at the center o f a p a tc h i s lo c a ted a t t o g iv e a general id e a o f t h e dimensions o f t h e fault. T h e final entries a r e t h e ran ges th at t h e disk-centers s p a n i n e a c h direction. A ll 12 surfaces u s e d p =1.0, h a d patches w it h a ra d iu s o f 1 unit, a n d u s e d s e e d number 3 3 3 w h ic h produced 5,360 d is k s f o r e a c h surface. Reproduced with permission o f th e copyright owner. Further reproduction prohibited without permission. Fractal D im ension Fault S urfa ce # x = 0,0 X = 4.0 X II O D © x = 12.0 x= 1 6 .0 A verage S tn d , D e v , 1 1 1 .2 5 1 .3 9 1 . 4 1 1 .3 2 1 .3 3 1 .3 4 0.06 12 1 .0 8 1 .2 4 1 .2 6 1.22 1 .1 9 1.20 0 .0 7 1 3 0.99 1.10 1 .2 5 1 .1 6 0.99 1.10 0.11 1 4 1 .3 6 1 .4 8 1 .4 3 1 .4 4 1 .4 5 1 .4 3 0 .0 4 1 5 1 .0 7 1 .2 9 1 .3 4 1 .2 9 1 .1 9 1 .2 4 0.11 1 6 0.95 1 .0 5 1 .1 8 1 .0 8 0.98 1 .0 5 0.0 9 1 7 1 .4 3 1 .5 9 1 .5 3 1 .5 0 1 .4 7 1 .5 0 0.0 6 1 8 1.20 1 .4 2 1 . 4 1 1 . 4 1 1 .3 3 1 .3 5 0.09 1 9 0.98 1 .0 9 1 .1 8 1 .1 1 1 .0 5 1 .0 8 0.07 20 1 .21 1 .3 7 1 .3 6 1 .3 7 1 .2 8 1 .3 2 0.0 7 2 1 1.00 1 .1 7 1 .2 5 1 .1 6 1 .0 6 1 .1 3 0.10 22 0.98 0.97 1 .1 3 0.97 0.99 1 .0 1 0.07 Table 2: Values for the fractal dimension of slices of the generated fault surfaces using a 2-D box counting technique. The slices are parallel to the y-z plane and span most of the x range of the fault surface. The estimated error for these value are +/- .05. The averages of the fractal dimensions calculated for slices of a particular surface are given along with the standard deviation to give a rough comparison of the results. These averages are ng| an estimate of the fractal dimension for the entire surface. o > o 61 plots suggest such a result, but additionally the probability distributions would naturally cause this to occur. When there is a higher probability of small angles of rotation the produced surface will stay relatively planar as patches stay in almost the same plane as their parent’s patch whereas higher chances of large rotations will cause more visible branching. An important aspect to note is that increasing the minimum rotation will cause large rotations whereas increasing the value of the exponent in the denominator of the power law will cause smaller rotations. By adjusting the two parameters accordingly similar fractal results can be obtained for different distributions as seen by the results in Table 2. This tends to show that Kagan’s choices for the best parameters are very non-unique and in fact several different power- law choices will produce similar effects. Thus using a Cauchy distribution is not as strong a choice as suggested by Kagan (1982) since it was specifically chosen due to Kagan’s preference for the results produced by the x*2 distribution.8 However, the Cauchy distribution can be similar to many different power laws with small enough values of a, so using it may still be best since it has the advantage of not needing a minimum cutoff. At this point it will be instructive to examine some actual plots that were used to determine the fractal dimensions. Three such examples are shown in Figure 12. Somewhat unexpectedly, most of the data sets were extremely linear, typically from 4 to 400 pixels length scales, which is desirable and works well for avoiding issues such as multiscaling. However, ‘ Kagan preferred the x*2 distribution because the results from using a four-point correlation technique (similar to the correlation dimension previously discussed) on surfaces created by the distribution were claimed to be the most comparable to hypocenter data from catalogs. Unfortunately, it is unclear if transforming the surface into a set of hypocenter points is a valid representation of the actual surface. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Ln(Boxes) vs. Ln(Side Length) Fault #11 at x = 0 8 7 6 ? s I 4 5 3 2 1 0 D - 1 .251 R" - 0 .9 9 7 3 0 2 4 6 8 L n (S td e L e n g th ) Ln(Boxes) vs. Ln(Side Length) Fault #21 at z - -1 4 3 2 1 0 0 . 0 .8 4 2 x I 'c ’ Ri - 0 .9 9 8 0 1 2 3 4 5 6 L n (S td e L e n g th ) Ln(Boxes) vs. Ln(Side Length) Fault #31 (3-0 boxcounting) 8 Tf^n996 I R * - 0 .9 9 9 4 ! 6 S 4 2 0 0 2 1 3 L n (S k te L e n g th ) Figure 12: Three fractal dimension plots using box counting methods. The top panel shows a linear plot for the box counting method of a fault surface slice which has a fractal dimension of 1.25. The middle panel shows a case of saturation for larger box sizes as seen by the plot leveling off at higher side length values (gray data points not included in best-fit line). The final linear fractal dimension plot is for a 3-D boxcounting technique discussed later, but included here for comparison. In general the closer RA 2 is to 1.0, the better the line fit Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 63 there were occasionally some plots which suggested some degree of saturation or multiscaling. Most of these cases arose when very few patches intersected the entire surface. Significant irregularities were usually noted in the tables when observed. Returning to the results from Table 2, it can then be summarized as follows. Any particular slice taken of the synthetic fault surface which is parallel to the y-z plane and thus perpendicular to the root patch, has a constant fractal dimension over a couple of orders of magnitude. The values of the slices though change throughout a particular surface, with larger values for the dimension tending to occur in the areas where the most patches occur. The lowest values are typically around 1.00 for the surfaces which have parameters that cause few deviations from a plane, such as fault surfaces number 13,16,19, and 22. At the other end of the spectrum the most erratic and unrealistic looking surfaces have values that approximately approach 1.60 such as surface number 17. Another approach to examining these surfaces through a 2-0 technique is to compare slices that are perpendicular to the x-axis to those perpendicular to the z-axis. To make this comparison I worked with the 12 surfaces used previously for Table 2. The fractal data for slices perpendicular to the z-axis are shown in Table 3. The first thing to be aware of is that the surfaces do not extend very far in the z-direction compared to how far they extended in the x-direction (Table 1). This is due to the small amounts of rotation keeping the faults usually rupturing in the x and y directions. This immediately limits how far the slices will be from the origin in the z-direction. Taking a first look at the results shows that any comparisons will be difficult for a variety of reasons. First, there is the problem that most of Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced w ith permission o f th e copyright owner. Further reproduction prohibited without permission. Fractal Dimension Fault Surface # Z = -2.0 z = -1.0 Z = 1.0 Z = 2.0 Average Stnd. Dev. 1 1 1.14 1.34 1.30 1.22** 1.25 0.09 12 0.94* 1.08* 1.11*' 1.06** 1.05 0.07 13 0.87* 0.99** 0.92 1.01** 0.95 0.06 14 1.28 1.50 1.28 1.22 1.32 0.12 15 X 0.98** 1.08** 1.02* 1.03 0.05 16 X X 1.00 X 1.00 NA 17 1.55 1.62 1.55 1.44 1.54 0.07 18 X 0.97 1.16 1.04* 1.06 0.10 19 X X X X NA NA 20 0.96 1.05* 1.18 1.27* 1.12 0.14 21 X 0.84* 0.98 1.00 0.94 0.09 22 X X X X NA NA X - The fractal dimension is not applicable since the surface does not pass through the location. * - Extra points were removed (usually those representing the larger box sizes) since significant saturation occurred for more than just the two points at either end of the scale. ** - Contains either a bend, significant curvature, or points that significantly deviate from the line, and thus it may not be best to fit a line to these data sets. Appears to be caused when only a small portion of the surface passes through the area. Table 3: Values for the fractal dimension of slices of the generated fault surfaces using a 2-D box counting technique. The slices are parallel to the x-y plane and span most of the z range of the fault surfaces. The estimated error for these value are +/- .05. The averages of the fractal dimensions calculated for slices of a particular surface are given along with the standard deviation to give a rough comparison of the results. These averages are ngi an estimate of the fractal dimension for the entire surface. o> a 65 the surfaces stay close to the z = 0 plane. Second, even when surfaces do stray at least a few units from the z = 0 plane, they are usually limited to only a small number of patches. This tends to cause plots with either significant saturation or data sets that are not necessarily linear as have been noted in Table 3. However, even with these difficulties there are at least two interesting observations that can still be made by this comparison. The first is a visual comparison as seen in Figure 13. The x-slices tended to create patterns that have many branches, but few gaps whereas the z-slices tend to show the opposite effect in that they usually show little branching, but quite a few gaps. The second comparison is that the values for the x-slices and z- slices seem to reach similar values for the same fault surface. However, it is important to note that the scales for each direction are quite different. Thus it appears that these surfaces, are best described as self-affine rather than self similar which is not surprising since fault surfaces tend to have preferred orientations. After examining the differences throughout the surface and those caused by parameter changes the issue of size changes was explored since based on the previous results larger faults may have larger values of D. Three different faults with sizes of 5,360,12,783, and 23,829 patches were examined with general attributes of the three shown in Tables 1, 4, and 5, respectively, and with comparisons of D values in Table 6. It turns out that trying to compare different fault sizes is problematic. Unlike studying different parameters or slice locations, different size trees use different seeds for the RNGs. This causes the creation of completely different trees with different branching structures and rotations. If values of D were the same throughout it Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 66 Fault Surface #11 5,360 patches At x = 0 D = 1.25 191 intersecting patches 2 . 5 i -10 -5 o 5 " At z = -2 D = 1.14 97 intersecting patches 0 10 20 x At x =16 D = 1.33 223 intersecting patches -10 -5 o 5 At z = 2 D = 1.22 165 intersecting patches 0 10 20 Figure 13: Four different slices of fault surface #11. Notice how the top two panels which are parallel to the x-plane shows significant branching, but few gaps whereas the bottom two panels which are parallel to the z-plane have the opposite behavior. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 67 9 § 2 N ® a c < a >. « a ' c < a N £ N C > X < f l E > c X X < a x c .g c p ® 4 3 a > •c 3 CO C O ao o P* . o> C O o> C D in © C D C O o o © co o t - i d C o ' C M C D o ir i C M d C O m C M Xp 3 * ( ■ © ao fx . px. C O C O 3 in fx . C M C O ir i C M 2 o C O 2 2 o C O in m C O o © © xp © C O C M o C M m C M C O o ir i C M o C O 8 co C M 8 o C O o C O n d C O 2 C M co (O C O C O 3 . C O o rx . in (O px C M d p*. ' o c m ' c m ' o © C O C M x r C O in C M f x» £ ao C O m C M © px. n : • C M • i n * ■ d • ■ m- O 1 • c m ' • p -i • d • C M s 3 m m 3 3 C M V r* . © o d C M rx l C D § ! G O ir i C M pxl © ' T f C O fx . 2 2 m © co px © C O co P > - cd I pxl • c m ' 1 C O ■ ■ C M * i C O 1 1 C M 1 3 © 2 X f C M © x r C M 3 2 3 ir i C M co C M in in C M C M 8 8 ir i C M P x . C M o> 2 © 2 CO C O o © © * in ' i 1 ir i 1 x r I X T 1 d • ir i 1 • in © r ^ op in C M C O < o < o « o < o < < o o < O < o < o in in in © o o in in in o o o co 32 33 34 35 36 37 38 39 © ^ © O C O C M C O * C M O *•10 0 a ® n v © o ' C M CO C O co o o C M * - O *■ 0 0 C O CO C O C M O O in in © in o ' o s s C M C M ' in co in o r s > ' a o co C M 3 a d in c m in co co c m ' c m ' o n S d i n i n C M C M C M *- © *- CM CO O in T t - • M ’ co ^ in • i * < < < o o o © 55^) C'O.OT ® + — © ~ 2 ®1< « — 1 1 £ 3 ^ 0 C O ^ ©B S Z — C O 0) © z © S H - O ■°-s > » 3 5 - c : o jT u eTa 2 © ® o -S © is o £ § C O O 3 CO a 8 CO O O ▼ “ CM ( 0 ® £ 0 CO a. s p ~ C M c a c 8 c o C M 1 rt < o § •c 3 c o © CO £ C O h E u . •2 ^ c/i _ c © .£ O 1 * 1 i l l 3 O 5 ® ® 3 £ o S ~ 1 * 2 © f 2 " S 2 o S § « £ -.§ ♦ ; © c S £ 6 © o ' © © o ^ 1 1 ? CD c co rt © J Q 2 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. instead. T h e fo llo w in g entries c o n tain t h e m in im u m a n d m a x im u m valu es th at the center o f a p a tc h i s lo c a te d a t t o g iv e a general id e a o f t h e dimensions o f t h e fault. T h e final entries a r e t h e ranges that t h e disk-centers s p a n i n e a c h direction. A ll 12 surfaces u s e d p =1.0, h a d patches w it h a ra d iu s o f 1 unit, a n d u s e d s e e d number 6 4 8 w h ic h produced 12,783 d is k s f o r e a c h surface. F a u lt S u r f a c e # d - v a l u e m in . rotn./a m in x maxx m i n y m a x y m ipz m a& z x r a n g e y r a n g e i ra ng e 68 I>x CO 0 5 3 i n CO co IX . m m CO 0 5 c o IX . CO IX . IX . CO 0 5 CO o ' e o o ' T f CO CO CM CO CM CM CM s in b CO c b o ' |X . fx * o o CM o p ? 0 5 CO IT ) c o IX . IX . 3 2 CO CM CM 0 0 C O CO CO CM x r 2 c b c v i n « CM CM o e o fx . CO CO 0 5 IX . 0 5 CM CM C O 8 CO CM m 5 i 0 5 O 2 CO CO e o CO CO CO 2 c b CO c b CM c b c o c b CO cm’ CO c b CO 0 5 in h - c o 2 s | x c p CO CO in |Xx 8 l>X m CM x f V o ' o ' CM s CM - |x l o ' c b x f o s CM o > 8 § CM V CM CM CO 8 8 CM m a i 1 a i • co • o ' • • 1 c b 1 0 0 • o ' 1 c o • • d i $ 8 8 in c o 8 CO CO i n 8 in c o CO i n CO r x i x f CO in ' c b c b in in i n 0 0 (X» X “ 0 5 0 5 CO 2 8 8 s 0 0 i n o ' CM • O ) i CO 1 CM I c b CM 1 c b CM 1 a i 1 o ' CO l |x l CM 1 c b CM • r x CM 1 c b CM I s CO CM CO CO 2 8 f f l 0 5 i n 0 5 § 8 e o i n CM 2 Si CM 2 CM 2 c b CM 2 c b CM 2 $ X “ in s 8 2 8 0 5 in 8 2 T “ C O m 2 | xJ © CM x f CM d c b cm ' cm ’ - CM CM in CO I X . CO T * 9 CM CO x f c o x f i n < o < o < o < o < < o o r - x * < • o < o < o < o < o < o 9 0 S O so 1.0 o o 1.5 1.5 1.5 C a uchy C a uchy C a uchy 5 CM in 8 2 55 8 57 8 59 09 r * CO 62 05 ® £ o ( D a. § G O m c v j 8 ( 0 C M C D I in os § £ 3 (0 £ ^.<2 ? ‘x « 1 JLS a s s © "O f i Q . . f a *1 C O (0 (I) S < D S “ O '°'S >* 3 © r : o j* O c T 3 3 © ® O A © iS'SS £ g fc h E u . . 5 2 o 0 ) o 3 si « ^ « 3 o ^ ,2 3 © o f i : ® £ ° i’ sg It'S cviio 2 ^ C O _ _ o 2 i i s --1 i n ill s i ? C D c © i n © n c o Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. instead. T h e fo llo w in g entries c o n ta in t h e m in im u m a n d m a x im u m values that the center o f a p a tc h i s lo c a te d a t t o g iv e a general id e a o f t h e dimensions o f t h e fault. T h e final entries a r e t h e ran ges th at t h e disk-centers s p a n i n e a c h direction. A ll 12 surfaces u s e d p =1.0, h a d patches w it h a ra d iu s o f 1 unit, a n d u s e d s e e d number 7 6 8 w h ic h produced 23,829 d is k s f o r e a c h surface. Reproduced with permission o f th e copyright owner. Further reproduction prohibited without permission. Fractal Dim ension Fault Surface # Patches x = 0,0 o c b n X 1 1 5,360 1.2 5 1 . 4 1 3 1 12,783 1.20 1.4 7 5 1 23,829 1.4 7 1.45 12 5,360 1.08 1.26 32 12,783 1.20 1.3 4 52 23,829 1.38 1.39 1 3 5,360 0.99 1.25 33 12,783 1.20 1.28 53 23,829 1.22 1.2 4 1 4 5,360 1.36 1.4 3 34 12,783 1.2 3 1 . 5 1 54 23,829 1.5 4 1.5 0 1 5 5,360 1.0 7 1.3 4 35 12,783 1.26 1.36 55 23,829 1.2 8 1.32 1 6 5,360 0.95 1.18 36 12,783 1.09 1.11 56 23,829 1.16 1.1 4 Fractal Dimension Fault Surface # Patches X I I O o o CO I I X 1 7 5,360 1 .4 3 1.5 3 37 12,783 1.3 3 1.57 57 23,829 1.6 0 1.4 3 1 8 5,360 1.20 1 . 4 1 38 12,783 1 . 3 1 1.42 5 8 23,829 1 .3 8 1.3 4 1 9 5,360 0.98 1.18 39 12,783 1.0 8 1.15 5 9 23,829 1.16 1.12 20 5,360 1.21 1.36 40 12,783 1.2 5 1.3 9 60 23,829 1.4 4 1.38 2 1 5,360 1.00 1.25 4 1 12,783 1.21 1.24 6 1 23,829 1.2 3 1.22 22 5,360 0.98 1.1 3 42 12,783 0.99 1.05 62 23,829 1.10 1.0 7 Table 6: Values for the fractal dimension of slices of different size fault surfaces using a 2-D box counting technique. The slices are parallel to the y-z plane taken at x = 0.0 and x = 8.0. The estimated error for these values are +/- .05. 70 would be simple to compare, but since they are not it is questionable if they are comparable at all. Making a comparison at x = 0.0 might be possible since the root disk causes all the faults to pass through x = 0.0. Looking at Table 6 shows that there may be a general trend towards an increasing D value with increasing size. This though is expected since there tends to be more intersecting disks. The comparison also shows that the variables causing the most planar surfaces tend to stay planar regardless of size whereas the more complex structures tend to become quite disordered as shown in Figure 14. Unfortunately the validity of this comparison is unknown since the nonrandom placement of the root may have too much influence on the geometry of the fault for the earliest created patches which are often located near the plane x = 0.0. For this reason an attempt was made to compare different size faults further away from the root patch. The plane x = 8.0 has been somewhat arbitrarily chosen since it is both away from the root patch and happens to be a location in which a significant portion of all the faults pass through. As seen in Table 6 there is no obvious correlation between tree size and fractal dimension. Though as stated earlier the comparison is misleading since the different size faults at a given location should be uncorrelated due to the randomness except for possibly near x = 0.0. It appears that the only way to get a good comparison between fractal values of faults of different sizes would be to get a 3-D fractal dimension of the entire surface. Such an approach will be discussed later. Whereas comparing different size faults may not work, analyzing the fractal dimension of a particular fault over time is possible, but has some difficulties that must be dealt with. The problem is that no real time scale has Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. - i o ! 71 Surface # 37 464 Patches D = 1.33 10 1 10 2 0 Surface # 57 2,699 Patches D = 1.60 10I 0 I -io r -20 -10 0 y Surface #39 i -------------------------------------------------- 398 Patches 2 _ [ j “ D = 1.08 L — -------------------5 --------------------- «- Surface #59 2,832 Patches D = 1.16 Figure 14: A comparison of the persistence (or lack thereof) of surfaces when using the more extreme parameter values. The top two panels have too much branching whereas the bottom two are essentially planar. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 72 been assigned to the branching trees at this point. Each generation happens at a later time than the generation that precedes it, but the actual amount of time that passes is currently arbitrary. There are several possible ways to resolve this in order to examine the fractal dimensions over time. One possibility would be to have 1 unit of time pass from every generation to the next. This though would not take into consideration the different number of patches created at each generation and would eliminate any randomness to the growth over time. A second possibility would be to not have actual values for the passing of time, but rather to consider the tree at different stages along its growth. For example, the surface could be examined after every additional 500 patches were added to the structure. This avoids the problem of attempting to designate a time to how long the process takes. However, this technique would be most effective if the entire surface was being examined in its entirety. Since the focus is on slices of the entire surface I preferred to deal with the passing of time by examining the slices after every n patches (for some number n) passed through the particular slice being examined. There is an additional consideration to be made though before actually plotting such data. There is the possibility of either examining the slices using “time windows” (e.g. examining patches 1-50, 51-100, etc.) or examining the entire surface up to that point (e.g. examining patches 1-50, 1- 100, etc.). Since these surfaces are supposed to represent faults, and to be occurring over a short enough period to neglect effects such as erosion, the latter approached was used. An example of how these surfaces grow over time is shown in Figure 15. Notice that there is a tendency for the surface to Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 73 2 | 1 After 50 Patches 2 01 — — D = 1.23 -Ij ^ ---------------o -----------------ro ----------- y After 200 Patches D = 1.31 2 0 * -2 - 4 -ttr Tor After 500 Patches D = 1.36 2 0 -2 - 4 After 1,000 Patches D = 1.38 2 0 -2 - 4 y After 1,494 Patches D = 1.36 2 0 2 4 iff y Figure .15: Time progression of a slice of fault surface #35 at x = 8.0. Of the 12,783 patches that make up the complete surface, 1,494 of them pass through x = 8.0. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 74 spread outwards along the x-y plane over time, but the surface also frequently backtracks onto itself. Some example plots which show the change in D versus the growth of the fault are presented in Figure 16. Somewhat unexpectedly the fractal dimension tends to level off in each case, and the values are generally approached asymptotically. This actually is somewhat similar to hypocenter data observed by Robertson et al. (1995). There are deviations from a perfect curve throughout, with the greatest deviations coming at the beginning. This is not surprising because there are only a few patches in a particular area at the beginning and they would tend to have a significant influence on the overall geometry in that area. It would be useful to have a curve fit which approximates these asymptotic approaches, but several different types of curves tend to approximate it equally, so no particular type of curve can be preferred. However, it seems that for several cases that were tried, the limiting fractal dimension is approximately reached by the time that half of the final number of patches in the area have been created. Unfortunately, this is a different number for each surface and it is not possible to know a priori when the halfway point is reached. Fault traces may be subject to similar conditions in that the fractal dimension may initially increase at the beginning and then stabilize in the long run. Once again this focus tends to lead to the question of whether these observations can be made about the entire surface simultaneously. Thus a 3-0 fractal approach is needed to see if it changes in the same asymptotic manner over time. Having a simulation which lets one see an entire surface has a significant advantage over field and experimental data which has limited Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Fractal Dimension vs. Number o f Patches for Fault #31 at x=8.0 1.6 S 1.5 o M C 1.4 3 1.3 (0 u 1.2 w 1.1 1 ♦ ♦ 500 1000 1500 2000 Number of patches that intersect x-8.0 2500 Fractal Dimension vs. Number of Patches for Fault #35 at x=8.0 1.45 1.4 S 1-35 S 1.3 1 1.25 ® 1.2 3 1.15 5 1 .1 1.05 1 500 1000 1500 Number of patches that intersect x=8.0 2000 Figure 16: Two plots showing the fractal dimension changing as time progresses. The bottom plot is for the same fault as that shown in Figure 15. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 76 exposure. Unfortunately, examining the entire synthetic surface as opposed to taken in slices is not as simple as it may seem. The problem which arises is that computers work with discrete space rather than continuous space. Although the slices represent the disks as line segments, the computer considers the segments to consist of a collection of pixels placed adjacent to each other. This is how the box-counting program is able to tell whether part of the segment falls within a particular box. Although the number of pixels used to represent a part of disk will depend on several considerations such as the size of the graphic, each line would be represented as a relatively small number of pixels. Representing the entire disk cannot be done with pixels, which are meant for 2-D objects. Points in 3-D space will work, but the number of points needed to accurately represent a disk is much larger than that needed to represent a line. It also has many boundary conditions to consider as opposed to only the two endpoints that a single segment needs. The tradeoff then is between having an accurate depiction of each patch by using more points versus the large amount of computer time and memory needed to process an entire complex structure at once. This issue has several interesting possibilities for future study such as how many points are needed to accurately depict the disk, and how those points should be arranged. It is not possible to fully explore the topic here, but since 2-D studies have been shown to have limitations a brief glimpse of the topic will be examined by looking at the simplest case. Each disk will be represented by a single point at its center and can possibly be though of as a “hypocenter,” for lack of a better term. This changes the geometry of the fault surface, but allows for a 3-D box-counting Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 77 dimension to be calculated. Note that the lack of points should lim it the lower bound of the fractal scaling range to be the radius of the disk since it will not recognize the disk itself. The box-counting program used for this part of the study has been provided to me and created by Shoshana Levin. I have made some small changes and additions to it so it can be used for this particular situation. Unlike the 2-D boxcounting program the boxes have sidelengths in units rather than in pixels, but spans less orders of magnitude with box sidelengths from about 1 to 15 units. The program was set to check 30 different box sizes and 64 positions at each size to minimize the number of filled boxes. Even just an initial look at the 3-D “hypocenter” fractal data shown in Table 7 provides some interesting results. The set of zero-dimensional hypocenters spread in 3-D space results in fractal dimensions both above and below 2.0. This immediately begins to show the problem of using epicenter data to approximate hypocenter data since epicenters cannot ever have a dimension above 2.0 since they all reside in the same plane. Additionally the results demonstrate how unintuitive the values of fractal dimensions are since the hypocenters do not fill 2-D space, but are probably near 2.0 since they enter 3-D space though do not adequately fill it. These values cannot accurately represent the final dimension of the entire fault surface as it is expected that such values will all be above 2.0 since the irregular surfaces take up more space than a plane. The 3-D fractal dimensions currently do not have a clear relationship to their 2-D counterparts. The fractal values increase as the faults become less planar, but they may be independent of fault size. How fast the Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Reproduced with permission o f th e copyright owner. Further reproduction prohibited without permission. Fault Surface # Fractal D Fault Surface# Fractal D Fault Surface # Fractal P 1 1 1.83 3 1 2.00 5 1 2.04 12 1.73 32 1.95 52 1.94 1 3 1.75 33 1.93 53 1.87 14 1.95 34 2.03 54 2.02 15 1.76 35 1.97 55 1.88 16 1.69 36 1.80 56 1.75 17 1.97 37 2.03 57 2.16 18 1.78 38 1.99 58 1.9 1 1 9 1.69 39 1.78 59 1.76 20 1.79 40 1.97 60 1.96 2 1 1.68 4 1 1.85 6 1 1.84 22 1.67 42 1.77 62 1.72 Table 7: Values for the fractal dimension of the patch centers of the generated fault surfaces using a 3-D box counting technique. T T ie estimated error for these values are +/- 0.01. 0 0 79 calculated values approach the actual values as more points are used to represent each patch, and how well another fractal method such as correlation dimension compares to the values in Table 7 are both questions worth examining, but will not be pursued here. Although the results presented have provided many insights into the fractal dimensions of nonplanar fault surfaces there are many concerns with the model that need to be addressed. One issue is that the surface is composed of many small circular patches. The problem with circular patches is that they cannot share sides to form a continuous surface. To resolve this, the patches are often overlapping and intersecting many times which weakens its ability to represent a physical object. Using alternate shapes for the patches such as a square would avoid this problem since they could be placed adjacent to each other and share an entire side. Squares unlike disks though create a bias since they do not have an unlimited number of axes of symmetry. The shape of the patch though may be irrelevant when the patches are small enough to avoid discretization effects. Clearly in the case with 5,000 patches and even with 23,000 there may still be discretization effects since the patches have a diameter of 2 units making it a significant percentage of the full length of the surface. Increasing the number of patches will decrease discretization effects, but there is still the question of the simplifying assumptions that have been made. For example X was set equal to 1.0 regardless of the state of the surface. However, it may make more physical sense to have X change as a function of time or length. Another possibility would be to not use a Poisson distribution and instead choose one that does not assume independence. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 80 The model presented here took an almost entirely stochastic approach to simulating reality. An alternate and opposite approach would be to have a purely physical, deterministic model. Reality probably falls somewhere in between the two and a better procedure may be to combine stochastic and deterministic fault models. Whether or not the model here is a dose enough depiction of reality, the complications that were shown when applying fractals to irregular surfaces should still apply. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 81 7 Conclusions This thesis began by recreating a model first developed by Kagan (1982) who in turn based the model off of work by Vere-Jones (1976,1977). At the time of Kagan’s work it was difficult to compare how well the model simulated earthquake fault geometry since it has only been within the last two decades that work has been done to quantitatively assess the irregularities of such phenomenon. This has been done using fractal dimensions. The original goal of this thesis was to see if Kagan’s model produced fractal dimensions comparable to lab and field data that has been published. Upon inspection of relevant papers, it was determined that too many difficulties occur in making such a conclusion. These issues include the various methods used to make the measurements, and the wide range of numerical results obtained. Thus, a second approach was taken which was to use this model to try and get a better understanding of fractal measurements of irregular synthetic fault surfaces. Two possible approaches could have been taken at this point. One would be to focus on a particular method for calculating fractal dimensions or on a particular aspect of the surface to measure (e.g. hypocenters or fault traces). However, this approach was not pursued due to the non-uniqueness of the model itself as seen by different sets of parameters creating similar surfaces. Instead the method taken was to look at many different aspects of the synthetic surfaces, sticking primarily with a 2-D box-counting technique. While this procedure prevented gathering large quantities of data in any one particular area, it did open up a wide range of concerns that should be explored in further depth. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 82 Although it cannot be determined whether or not the surfaces created here provide a realistic comparison to actual fault surfaces, it did not prevent important realizations from being made about how careful one needs to be when making fractal measurements. Fractals are a concept that originated with mathematical objects and are not as straightforward when working with naturally occurring objects. For example, a fault in the subsurface is unlikely to be easily related to its surface trace as it has been shown that 2-D slices of an irregular object are not easily related to the overall 3-D object. Thus it appears unrealistic to try and infer the 3-D fractal dimension by adding 1 to a value obtained from a 2-D slice or to determine the capacity dimension of epicenters or fault traces using the b-value or p-value, respectively. Additionally, altering earthquake geometry such as using epicenters or replacing surfaces with hypocenters should be avoided. There are two definite objectives for future work. They are to fully characterize the comparability of different fractal techniques and to determine how much information can really be gathered by using such techniques on irregular surfaces. The value of using fractal dimensions to characterize earthquake fault geometry is currently unclear. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 83 B ib lio g ra p h y Aviles, C.A., C.H. Scholz, and J. Boatwright, Fractal analysis applied to characteristic segments of the San Andreas fault, J. Geophys. Res., 92, 331-344, 1987. Avnir, D., O. Biham, D. Lidar. and O. Malcai, Is the geometry of nature fractal?, Science, 279, 39-40, 1998. Ben-Zion, Y. and C.G. Sammis, Characterization of fault zones, Pure Appl. Geophys., in press, 2001. Davy, P., A. Sornette, and D. Sornette, Experimental discovery of scaling laws relating fractal dimensions and the length distribution exponent of fault systems, Geophys. Res. Lett., 19, 361-363,1992. Eneva, M., Monofractal or multifractal: a case study of spatial distribution of mining induced seismic activity, Nonlinear Processes in Geophysics, 1, 182-190, 1994. Fisher, D.S., K. Dahmen, S. Ramanathan and Y. Ben-Zion, Statistics of Earthquakes in simple models of heterogeneous faults, Phys. Rev. Lett., 78, 4885-4888, 1997. Frohlich C., and S.D. Davis, Teleseismic b values; or, much ado about 1.0, J. Geophys. Res., 98, 631-644,1993. Goltz, C., Fractal and Chaotic Properties of Earthquakes, 178pp., Springer, Berlin, 1997. Gonzato G., F. Mulargia, and W. Marzocchi, Practical applications of fractal analysis: problems and solutions, Geophys. J. Int., 132, 275-282, 1998. Gonzato G., F. Mulargia, and M. Ciccotti, Measuring the fractal dimension of ideal and actual objects: implications for application in geology and geophysics, Geophys. J. Int., 142,108-116, 2000. Hirabayashi, T., K. Ito, and T. Yoshii, Multifractal analysis of earthquakes, Pure Appl. Geophys., 138, 591-610, 1992. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 84 Hirata, T., Fractal dimension of fault systems in Japan: Fractal structure in rock fracture geometry at various scales, Pure Appl. Geophys., 131, 157-170, 1989. Kagan, Y.Y., Stochastic model of earthquake fault geometry, Geophy. J.R. astr. Soc., 71, 659-691, 1982. Kagan, Y.Y., Random stress and earthquake statistics: spatial dependence, Geophys. J. Int., 102, 573-583,1990. Kagan Y.Y., Fractal dimension of brittle fracture, J. Nonlinear Sci., 1,1-16, 1991. Kagan, Y.Y., Observational evidence for earthquakes as a nonlinear dynamic process, Physica D, 77,160-192, 1994. Kagan, Y.Y., Are earthquakes predictable?, Geophys. J. Int., 131, 505-525, 1997. Karlin, S. and H.M. Taylor, A First Course in Stochastic Processes, 557pp., Academic Press, San Diego, 1975. Lawler, G.F., Introduction to Stochastic Processes, 176pp., Chapman & Hall/CRC, Boca Raton, 1995 (reprinted 2000). Lay, T., and T.C. Wallace, Modem Global Seismology, 521pp., Academic Press, San Diego, 1995. Lei, X., O. Nishizawa, and K. Kusunose, Band-limited heterogeneous fractal structure of earthquakes and acoustic-emission events, Geophys. J. Int., 115, 79-84,1993. Leon-Garcia, A., Probability and Random Processes for Electrical Engineering, 596pp., Addison-Wesley Publishing Company, Reading Massachusetts, 1994. Malcai, O., D.A. Lidar, O. Biham, and A. Avnir, Scaling range and cutoffs in empirical fractals, Phys. Rev. E, 56, 2817-2828,1997. Mandelbrot, B.B., The Fractal Geometry of Nature, 468pp., W.H. Freeman and Company, New York, 1983. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 85 Matsumoto, N., K. Yomogida, and S. Honda, Fractal analysis of fault systems in Japan and the Philippines, Geophys. Res. Lett., 19, 357-360,1992. Nanjo, K., H. Nagahama, and M. Satomura, Rates of aftershock decay and the fractal structure of active systems, Tectonophysics, 287,173-186, 1998. Okubo, P.G. and K. Aki., Fractal geometry in the San Andreas fault system, J. Geophys. Res., 92, 345-355,1987. Ouillon, G., C. Castaing, 0. Sornette, Hierarchical geometry of faulting, J. Geophys. Res., 101, 5,477-5,487,1996. Power, W.L., and T.E. Tullis, Euclidean and fractal models for the description of rock surface roughness, J. Geophys. Res., 96, 415-424,1991. Press, W.H., S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Numerical Recipes in Fortran 90, 552pp, Cambridge University Press, Cambridge, 1996. Rice, J.A., Mathematical Statistics and Data Analysis (second edition), 602pp, Duxbury Press, Belmont California, 1995. Robertson, M.C., C.G. Sammis, M. Sahimi, and A.J. Martin, Fractal analysis of three-dimensional spatial distributions of earthquakes with a percolation interpretation, J. Geophys. Res., 100, 609-620,1995. Rubin, A.M., D. Gillard, and J. Got, Streaks of microearthquakes along creeping faults. Nature, 400, 635-641,1999. Segall, P. and D.D. Pollard, Nucleation and growth of strike slip faults in granite, J. Geophys. Res., 88, 555-568,1983. Sharon, E., and J. Fineberg, Microbranching instability and the dynamic fracture of brittle materials, Phys. Rev. B, 54, 7128-7139,1996. Turcotte, D.L., A fractal model for crustal deformation, Tectonophysics, 132, 261-269 , 1986. Turcotte, D.L., Fractals and Chaos in Geology and Geophysics, 398pp, Cambridge University Press, New York, 1997. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 86 Vere-Jones D., A branching model for crack propagation, Pure Appl. Geophys., 114, 711-725, 1976. Vere-Jones D., Statistical theories of crack propagation, Math. Geol., 9, 455- 481, 1977. Waldhauser, F., and W.L. Ellsworth, A double-difference earthquake location algorithm: method and application to the northern Hayward fault, California, Bull. Seism. Soc. Am., 90, 1353-1368, 2000. Xu T., I.O. Moore, and J.C. Gallant, Fractals, fractal dimensions and landscapes - a review, Geomorphology, 8, 245-262,1993. Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 87 APPENDIX A: Program Polpro18.f90 ! P0IPR018. F90 i By E ric L ib ic k i l 6 /2 /0 1 f 1 This program is based on p a rt o f Kagan's 1982 paper. This program creates ! d isks which are in te rre la te d and combined to represent a nonplanar fa u lt su rface. I The disks are e s s e n tia lly v e rtic e s o f a tre e th a t are stru ctu red w ith ! a Poisson d is trib u tio n . Bach new d is k is assigned values to show it s ! p o s itio n in space and re la tio n s h ip to o th er d is k s . Bach d isk is ro tated ! based on power-law o r Cauchy d is trib u tio n s . i - Mean is used instead o f Kagan's parameter o f c r it ic a lit y . ] - d v a lu e makes th e exponent o f th e power-law l+d_value as in Kagan (1982). ! ! n rtyp e, n r, poidev, ra n i, and ran _state are used fo r random number generators t and taken from Numerical Recipes (Press e t a l . , 1996) M O D U L E points U SE nrtype USE n r, ONLY:poidev, ra n i USE ran _state IMPLICIT N O N E S A VE INTEGER, PARAMETER:: REAL:: meanx TYPE:: v ertex INTEGER:: id INTEGER:: gener INTEGER:: ch ild ren INTEGER:: parent REAL:: xcor ycor zcor ph i th e ta REAL: REAL: REAL: REAL: END TYPE TY P E (vertex), DIMENSION(N):: INTEGER:: maxidn INTEGER:: gin REAL:: radius REAL:: dvalue REAL:: m in ro tatio n N-100000 lAvoids tre e s w ith more than 100000 v e rtic e s [The mean value fo r th e Poisson d is trib u tio n tBach v e rte x has many numbers associated w ith i t (Each v e rte x gets it s own id number !G eneration number !Number o f c h ild re n th a t v e rte x has ! ID number o f parent ix_coordinate o f cen ter o f "disk" !Y_coordinate o f cen ter o f "disk'* lZ_coordinate o f cen ter o f "disk" [p h i a n g le (ra d ia n s )-ro ta tio n o f vector normal to disk lth e ta a n g le (ra d ia n s )-ro ta tio n o f vecto r normal to d isk v e rtic e s 1A vecto r o f N v e rte x [Maximum id number given o u t so fa r (G enertaion le v e l number [The radius o f th e c irc u la r d islo catio n s t Value fo r "d" in the exponent o f th e ro ta tio n d is trib u tio n [Minimum value fo r p h i in the ro ta tio n d is trib u tio n CHARACTER(len * 2 0 ):: o u tfile [The f i l e name fo r th e output END M O D U L E P R O G R A M poiprog Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 88 U S E points IMPLICIT N O N E INTEGER:: d is tr ib INTEGER:: i INTEGER:: j INTEGER:: k INTEGER:: re in itlo o p INTEGER:: z = l INTEGER:: seqno INTEGER:: o u tp u ts ty le iThe number produced by POIDBV(meanx) {Loop index I Loop index ILoop index 1 R e in itia liz in g loop index !A counter to avoid an in f in it e loop {Sequence number fo r random Poisson d is trib u tio n {Choice o f how th e output should be presented w r ite ( * , ' (A )1) 'E ric L ib ic k i’ WRITE ( * , ' (A /)■) 'POISSONX.OUTY' W R ITB (*,*) ' Please choose a mean value fo r th e Poisson d is trib u tio n ,' W R ITE (*,*) ’ don’ ' t fo rg e t the d ecim al.' READ(*,+) meanx WRITE( * , * ) ' You have chosen ' , meanx,' fo r th e mean o f the Poisson' WRITE( * , * ) ' d is tr ib u tio n .' WRITE(* ,* ) W RITE(*,*> ' Please choose an in teg er to in it ia liz e the Poisson' WRITE( * , * ) ' random number g e n e ra to r.' READ(*,*> seqno WRITE( * , * ) ' You have chosen ' , seqno,' fo r th e in it ia liz a t io n ' WRITE( * , * ) WRITE(* , * ) ' Please choose a radius (UNITLESS FO R N O W ) fo r the a l l o f th e ' WRITE( * , * ) ' c irc u la r d is lo c a tio n s , d o n ''t fo rg e t the d ecim al.' READ(*,* ) radius WRITE( * , * ) ' You have chosen ' , ra d iu s ,' fo r th e radius o f the d is lo c a tio n .' W R ITE(*,*) w r it e( * , * ) ' Please choose a value fo r "d", which is p a rt o f the exponent' WRITE(* , * ) ' in th e ro ta tio n a l d is trib u tio n o f th e d iscs, d o n ''t fo rg e t' W R ITE (*,*) ' the decim al. I f you want a Cauchy d is trib u tio n instead o f a ' WRITE( * , * ) ' power-law d is trib u tio n choose ' ' 9 . 9 ' ' . ' READ(*,*) dvalue WRITE( * , * ) ' You have chosen ' , d v a lu e ,' fo r "d ". ' WRITE(* ,* ) w ritE (* , * ) ' Please choose a value fo r the minimum ro ta tio n angle o f' WRITE( * , * ) ' the d is k s , don' ’ t fo rg e t the decim al. I f you have chosen' WRITE( * , * ) ' a Cauchy d is trib u tio n , type in th e parameter value you w a n t.' READ(*,*) m inro tation WRITE(* , * ) ' You have chosen ' , m in ro ta tio n ,' fo r the minimum ro ta tio n a n g le .' W RITE<*,*) WRITE( * , * ) ' Please choose th e form at fo r th e o u tp u t’ W RITE(*,») ' l ) P resen tation s ty le ' w r it e( » ,* ) ' 2) R o ta tio n a l output s ty le ' W R ITE (*,*) ' 3) Hypocenter output s ty le ' REM >(*,*) o u tp u ts ty le IF (o u tp u t_style = 1) WRITE( * , * ) ’ You chose presentation s ty le ' IF (o u tp u t_style “ 2) WRITE( * , * ) ' You chose ro ta tio n a l output s ty le ' IF (o u tp u t_style “ 3) WRITE( * , * ) ' You chose hypocenter output s ty le ' W R ITE(*,*) WRITE( * , * ) 'Please choose a name fo r th e o utput f i l e ' READ(*,») o u tfile W R ITE (*,*) 'You have chosen ' , o u t f ile ,' fo r th e f i l e name.' Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 89 W RTTE(*,*) O P E N (UNIT = 15, FILE = o u tfile , STATUS = 'new ', ACTION = 'w r ite ') U n itia liz in g Poisson random number generator c a ll R A N _SEEO (S EQ UEN CE = seqno) ! In it ia liz in g gener and id values fo r th e ro o t maxidn “ 1 v e rtic e s (l)% id * maxidn v e rtic e s (l)% g e n e r * 0 v e rtic e s (l)% p a re n t = 0 v e rtic e s (l)% x c o r = 0.0 v e rtic e s (l)% y c o r = 0.0 v e rtic e s ( l)% zcor * 0.0 v e rtic e s (l)% p h i • 0.0 v e rtic e s (1 )%theta » 0 .0 ! maxidn s ta rts a t 1 {ro o t's id is 1 (ro o t is generation zero {ro o t's parent is ’ ’ zero” {ro o t’ s cen ter x coordinate is {ro o t's cen ter y coordinate is (ro o t's cen ter z coordinate is (ro o t's p h i value is 0 .0 (ro o t's th e ta value is 0 .0 0.0 0.0 0.0 (The fo llo w in g makes sure th a t the ro o t has o ffs p rin g do d is tr ib * INT(poidev(m eanx)) IF (d is tr ib / - 0) TH E N v e rtic e s (1 )%children - d is tr ib EXIT ELSE z=z+l 1temporary check to make sure d it r ib d o esn 't always = 0 IF (Z “ 10) T H E N WRITE( * , * ) ' ***R o o t may not ever have o ffs p rin g ***' S TO P E N D IF B N D IF E N D D O ( In it ia liz in g o ffs p rin g D O j» l,N IF (v e rtic e s ( j ) % children /= 0) THEN CALL o ffs p rin g _ o f_ id n (v e rtic e s ( j)% id , v e rtic e s ( j)% c h ild re n , maxidn) END IF E N D D O (Updates values fo r the o ffs p rin g D O k=2,maxidn CALL ran d o m _circle(vertices(k)% id , vertices(k)% xco r, vertices(k)% yco r,« & v e rtic e s ( k)%zcor) CALL update_xy* (v e rtic e s (k ) % id , v e rtic e s ( k ) %xcor, vertices (k )% yco r,a f c vertices(k)% zco r) CALL p h i_ th e ta ( v e rtic e s ( k ) %id, dvalu e, m in ro ta tio n , v e rtic e s (k)% phi, a a v e rtic e s ( k ) tth e ta ) EN D D O Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 90 !Determines output IF (output s ty le = 1) TH E N WRITE (15,150) 150 FO R M A T ( IX , T1, ' iDnumber' , T10, ' Generatn1, T20, ' C hldrn •, T28, ■ ParntID ’ , 6 6T 40,‘x_cor’ , T 5 0 ,'y _ c o r',T 6 0 ,'z _ c o r' ,T 6 9 ,‘p h i’ ,T 7 4 ,'th e ta ') WRITE (15,155) 155 F O R M A T ( ' = = *r-z = = *-------------------------------- — i---------------------------- ■ ■ - - k * ) D O i = l , maxidn WRITE (15,160) v e rtic e s ( i ) % id ,v e rtic e s ( i ) % gen er,vertices( i ) % children, k k v e rtic e s (i)% p a re n t, v e rtic e s (i)% x c o r, v e rtic e s (i)% y c o r,6 k v e rtic e s (i)% z c o r, v e rtic e s (i)% p h i, v e rtic e s ( i ) %theta 160 F O R M A T (T3, 16, T13, 15, T22, 14,6 6T29, 16, T38, F 8 .4 , T47, F8.46 4T57, F 8 .4 , T66, F 6 .4 , T74, F 5.3) E N D D O ELSE IF (output js ty le “= 2) T H E N D O i=l,m axidn WRITE(15, 170) -1 .0 * v e rtic e s (i)% p h i, -1 .0 *v e rtic e s (i)% th e ta , 6 6 v e rtic e s (i)% x c o r, v e rtic e s (i)% y c o r, v e rtic e s (i)% zc o r 170 F O R M A T (T2, F 10.6, T16, F 10.6, T30, F 10.6, T44, F10.6, T58, F 10.6) E N D D O ELSE IF (output s ty le “ 3) THEN D O i = l , maxidn w r it e( 15, 180) v e rtic e s (i)% x c o r, v e rtic e s (i)% y c o r, v e rtic e s (i)% zc o r 180 F O R M A T (T5, F 1 0 .6 , T20, F 10.6, T35, F 10.6) E N D D O E N D IF ! R e in itia liz e s "v e rtic e s ” vecto r D O re in itlo o p = 1 , maxidn v e rtic e s (re in itlo o p )% id = 0 v e rtic e s ( reinitloop)% gener - 0 v e rtic e s ( re in itlo o p )% children “ 0 v e rtic e s ( re in itlo o p )%parent = 0 v e rtic e s ( rein itlo o p )% xco r ■ 0 .0 v e rtic e s ( re in itlo o p )%ycor = 0 .0 v e rtic e s ( re in itlo o p )%zcor = 0 .0 v e rtic e s ( re in itlo o p )% p h i =0.0 v e rtic e s ( re in itlo o p )% th e ta ” 0.0 E N D D O CLOSE (UNIT = 15) END P R O G R A M { * » » • • * * * * • • * • • * * * * • * * * • » * • • • * • • * • * * * • • • * • * * • • • • • • • • • • * • • • * • • • * • • • • » * • * • • • • • * « :Generates o ffsp rin g Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 91 SUBROUTINE o ffsp rin g _o f_id n (id n u m b er, num ber_of_children, naxnumber) U S E points IMPLICIT N O N E INTEGER, INTENT( IN ) : : id_number !Id number INTEGER, INTENT( IN ):: num ber_of_children !Number o f ch ild ren INTEGER, INTENT( INOUT): : maxnumber !Max idn given out so fa r INTEGER:: k ILoop index D O k * l, n is n fa e ro fc h ild re n maxnumber - maxnumber+1 IF (maxnumber > N) THEN WRITE( * , * ) ' stopping program, A rray going beyond ', N S T O P E N D IF v e rtic e s (maxnumber) %id « maxnumber v e rtic e s (maxnumber)*children * INT(poidev(meanx)) v e rtic e s (maxnumber )%parent = id_number v e rtic e s (maxnumber)%gener = & a v e rtic e s ( v e rtic e s (maxnumber) %parent) %gener+l I(g e n e ra tio n is one more than parent) E N D D O EN D SUBROUTINE lAssigns i n i t i a l x ,y , and z values randomly on the perim eter o f th e p aren t disk I * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * SUBROUTINE ra n d o m c ire le ( id_number, x coordinate, y c o o rd in a te , z c o o rd in a te ) U S E points IMPLICIT N O N E INTEGER, INTENT( IN ) :: id_number REAL, INTENT(OUT) : : x_coordinate REAL, nn'ENT(OVT):: yjco o rd in ate REAL, INTENT(OUT):: zc o o rd in a te REAL:: m !Output fo r Random Number Generator CALL RANl_S(m) xjco o rd in ate » radius*C O S(2.0*PI*m ) + v e rtic e s (vertices(id_num ber)% parent)% xcor yjco o rd in ate ” radiu s*S IN (2 .0*PI»m ) + v e rtic e s (v e rtic e s (id_num ber)tparent)% ycor zjco o rd in ate = v e rtic e s ( v e rtic e s (id_number)%parent)%zcor EN D SUBRO UTINE !updates x ,y , and z coordinates based on Phi and Theta o f p aren t. SUBROUTINE update_xyz( id_nixnber, new_x, newjy, newjz) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 92 U S E po ints IMPLICIT N O N E INTEGER, INTENT( IN ) : : id_nunber REAL, ItfTENT( INOUT) : : new x REAL, INTENT( INOUT):: new_y REAL, INTENT( INOUT) : : new z REAL:: P_phi (p h i value o f parent REAL:: p_theta (th e ta value o f the parent REAL:: P_* (x_center o f parent REAL:: p_y (y_center o f parent REAL:: P_z (z_center o f parent REAL:: te n p x (Takes on value o f new_x REAL:: temp_y (Takes on value o f new_y REAL:: temp_z (Takes on value o f new_z p_phi • v e rtic e s (vertices(id_nuraber)% parent)% phi p_theta “ v e rtic e s (vertices(id_num ber)% parent)% theta p_x « v e rtic e s (vertices(id_num ber)% parent)% xcor p_y = v e rtic e s (vertices(id_num ber)% parent)% ycor p_z = v e rtic e s ( vertices(id_num ber)% parent)% zcor tetqp_x = new_x tenp_y * new_y temp_z * new z new_x * p_x+(-p_x+tem p_x)*C 0S (p_theta)-(-p_y+tenp_y)*S IN (p_theta)*C 06(p_phi) + & * (-p _z+ tem p _z)*slN (p _th eta)*slN (p _p h i) new_y = p_y+ ( -p_x+tesrp_x) *SIN ( p _th eta)+ (-p_y+temp_y) *COS ( p _theta) *COS (p _p h i) - a « ( -p_z+tem p_z)*COS(p _ th e ta ) *S IN ( p_p hi) new z = p_z+ ( -p_z+temp_z) *COS ( p_phi)+ (-p_y+tesnp_y) *SIN ( p_phi) E N D SUBROUTINE !Assigns Phi and Theta values based on d and minimum ro ta tio n S U B R O U TI N E p h ith e ta ( idnum ber, d_value, m in r o t, phi_value, th eta_valu e) U S E po ints IMPLICIT N O N E INTEGER, INTENT( IN ):: id_number REAL, INTENT( IN ): : d_value REAL, INTEN T(IN): : m in ro t REAL, INTENT(OUT): s phi_value REAL, HRENT(OUT):: theta_value REAL:: pd (C u rre n tly a Power d is trib u tio n ( l/x ~ ( l+d _valu e)) REAL:: m (Output fo r Randan Number Generator REAL:: posneg (E ith e r +1.0 o r -1 .0 (cw/ccv) based on R N G REAL:: p_phi (p h i value o f parent REAL:: p _th eta (th e ta valu e o f th e parent REAL:: neg_d_value (negative o f th e d_value REAL:: re c ip d value (re c ip ro c a l o f th e d_value Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 93 ! Follow ing is a dependent d is trib u tio n th a t is determined by the t values o f "d_value" and "m in_rot" entered by the user. p_phi * v e rtic e s (v e rtic e s ( idnu m ber) %parent) %phi p th e ta « v e rtic e s (v e rtic e s ( id_number) %parent) %theta neg_d_value * -1 .0 *d_value r e c ip d v a lu e = 1.0/d _value IF (d_value = 9 .9 ) THEN !Cauchy d is trib u tio n was chosen CALL RANl_S(m) IF (m < .5 ) THEN posneg “ -1 .0 ELSE posneg = 1.0 E N D IF CALL RANl_S(m) pd * min_rot*TAN ( ( P I/2 ) *m) phi_value “ M O D U L O ( posneg-ATAN( pd)+p_phi, 2 .0 *P I) IKagan's O P norm alized to 1 CALL RANl_S(m) IF (m < .5 ) TH EN posneg = -1 .0 ELSE posneg « 1.0 E N D IF CALL RANl_S(m) pd = m in_rot*TA N ((P I/2)*m ) theta_value = M O D U L O (posneg*ATAN( pd)+p_theta, 2 .0 *P I) 1 Kagan' s O P norm alized to 1 ELSE CALL RANl_S(m) IF (m < .5 ) THEN posneg = -1 .0 ELSE posneg = 1.0 E N D IF CALL RANl_S(m) IF (d v a lu e “ 0 .5 ) TH EN m « (1 .0 - .29289321)*m t .29289321 ELSE IF (d_value = 1 .0 ) THEN m = .5*m + .5 i Hakes valu e o f ”m " be between .5 and 1 so x >= m inrot ELSE IF (d_value = 1 .5 ) THEN m * (1 .0 - .646447)*m + .646447 else WRITE( * , * ) "You have choosen a d_value th a t c a n 't c u rre n tly be processed” S T O P E N D IF pd = 1 . 0 / ( ( (m in_rot**negL_d_value)-( (m in_rot**neg_d_value)*m ) )**recip _d _v alu e)6 S-m in_rot phi_value * M O D U L O ( posneg*ATAN(pd)+p_phi, 2 .0 -P I) CALL RANl_S(m) IF (m < .5 ) THEN posneg = -1 .0 Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 94 ELSE posneg = 1.0 E N D IF CALL R A N I S (m) IF (d v a lu e «« 0 .5 ) TH E N m ■ (1 .0 - .29289321)*m + .29289321 ELSE IF (d v a lu e = 1 .0 ) TH E N m = .5*m + .5 i Makes value o f "m " be between .5 and 1 so x >= m inrot ELSE IF (d_value = 1 .5 ) TH E N m * (1 .0 - .646447)*m + .646447 ELSE W R I T E(*,*) "You have choosen a d_value th a t c a n 't c u rre n tly be processed" S T O P E N D IF pd = 1 .0 /( ( (m in_rot»*neq_d_value)-((m in ro t**n e q d value)*m ))*"recip_d_value)fc 4-m in_rot th eta_value » M O D U L O ( poBneg*ATAN(pd) +p_theta, 2 .0 *P I) E N D IF E N D SUBROUTINE Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. Appendix B: Program Sllces5.nb 95 SLICES5 .NB; By E ric L ib ic k i; 6/ 2/ 0 1; This is version 5 o f a Mathematics program which g ra p h ic a lly represents; 2 - D s lic e s o f a fa u lt surface created by the data from p o ip ro l8 .f9 0 .; NOTE: This version is s p e c ific a lly fo r s lic e s p a ra lle l to the y - z p la n e .; NOTE: The v a ria b le s "radius" and "disks" need to be changed by hand; General Comnents; - The comnand O ff [G en eral:: s p e ll] w ill surpress the e rro r th a t th ere are; s im ila r v a ria b le names; - The values o f "rad iu s ," "P lo tran g e," and "Viewpoint" may s ig n ific a n tly ; a ffe c t the tim ing o f th e program; - The comnand which ro ta te s th e disks M U S T be done before tra n s la tin g the; disks otherw ise i t w ill ro ta te th e tra n s la te d d isk around 0, 0, 0; - The w idth o f the xperpplane is c u rre n tly .0002 u n its ; - The fo llo w in g comnand is used to exp ort the graphic; to a bitmap f i l e so i t can be used w ith the boxcounting program; Tim ing[Export["pl817x0.bnp", p ic tu re , {ImageSize -> {1024, 1 0 2 4 }}]]; Timing[ FILENAME = " p o is l8 .o u tll”; ( * Nam e o f in p u t f i l e *) OpenRead[FILENAME] ; SetStreamPosition[FILENAME, 0 ]; radius * 1 .0 ; ( * The raid u s o f every d isk *) disks = 5360; ( * The number o f disks in th e f i l e used fo r the p lo t * ) xperpplane = 0 .0 ; ( * The s lic e taken perpendicular to the x - axis *) xminp = xperpplane - 0.0001; ( * The min x value in th e p lo t * ) xmaxp = xperpplane + 0.0001; ( * The max x value in th e p lo t * ) a n g le l “ 0 .0 ; ( * The ph i angle ro ta tio n o f th e d is k * ) angle2 >0.0; (* The th e ta angle ro ta tio n o f th e d is k * ) xcor = 0 .0 ; (* The x coordinate fo r th e c e n te r o f th e d isk * ) ycor = 0 .0 ; ( * The y coordinate fo r the c e n te r o f th e d isk *) zcor = 0 .0 ; ( • The z coordinate fo r the c e n te r o f th e d isk * ) (* In te rv a l value must go evenly in to 2 .0 so th e c ir c le ; s ta rtin g p o in t is repeated a t th e end; *) in te rv a l = .0 5 ; po ints = 2 /in te rv a l + 1; corA * corB = corC = corO * { 0 .0 , 0 .0 , 0 .0 }; corE > corF = {0 .0 , 0 .0 , 0 .0 }; Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 96 c o rl = cor2 cor4 = cor5 lis tA = {>; lis tB = {>; segmentAB; percentAB; segmentCD; percentCD; segmentEF; i ; j ; k; m ; cor3 =0.0; (* x, y , and z values o f corE *) cor6 = 0 .0 ; (* x , y, and z values o f corF *) * C reates a l i s t o f points which represents a d is k /c irc le *) * Keeps a l l the lin e graphics in a ta b le to p lo t to geth er *) * Segment o f the d is k /c irc le *) * used to fin d the p o in t on segmentAB th a t is on the 2 - D s lic e * ) * Segment o f th e d is k /c irc le *) * Used to fin d the p o in t on segmentCD th a t is on the 2 - D * Segment o f the in te rs e c tio n o f th e patch and 2 - D s lic e * ta b le v a ria b le fo r c irta b le * ) * loop v a ria b le *) * loop v a ria b le * ) * loop v a ria b le * ) s lic e * ) *) ( * Keeps tra c k o f max t min values o f in p u t and used fo r the p lo t range; in it ia liz e d to s a fe , opposite extreme values; *) maxxcor = -9 9 9 9 .9 ; minxcor = 9999.9; maxycor = -9 9 9 9 .9 ; minycor = 9999.9; maxzcor » -9 9 9 9 .9 ; minzcor = 9999.9; < * The fo llo w in g comnand represents a d isk by using many points along; it s perim eter; - I t is im portant to repeat the s ta rtin g p o in t o f the c ir c le ; - Need to use th e comnand Chop because otherw ise has 2*10~-16 ra th e r than zero; - Have to use .12 (o r o th er a lte rn a tiv e ) since i t does not take zero e x a c tly ; *) c irta b le = T ab le[ {C h o p [ra d iu s *S in [i*P i]] , Chop( radius*C os[ i* P i] ] , 0>, { i , 0 .1 2 , 2 .1 2 , in te rv a l} ]; For[m = l , m <= d is k s , mt+, ( ( * Reads in p h i, th e a ta , xcor, ycor, and zcor from filenam e *) a n g le l “ Read[FILENAME, R eal] ; ( * negative alread y taken by p o ip ro l8 .f9 0 * ) angle2 = Read[FILENAME, R eal] ; ( * negative alread y taken by p o ip ro l8 .f9 0 * ) xcor = Read[ FILENAME, R eal]; ycor = Read[FILENAME, R eal]; zcor “ Read[ FILENAME, R eal]; (* Keeps tra c k o f max and min xcor, yco r, and zco r, which represent th e; o u ter perim eter o f th e fa u lt surface; * ) maxxcor » Max[maxxcor, x co r]; minxcor - M in[m inxcor, xco r]; maxycor Max[maxycor, yco r]; minycor = Min[m inycor, yco r]; maxzcor = Max[maxzcor. zc o r]; minzcor = Min [m inzcor. zc o r]; Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. I f [ (xcor > (xper pplane - ra d iu s )) a (xcor < (xperpplane + ra d iu s )), { lis tA * {>; (* Creates a l i s t o f p oin ts to represent an in d iv id u a l d is k * ) F o r[j = 1, j < p o in ts , j+ + , { lis tA = Append [ L istA , Rotate3D[ P a rtf c ir ta b le , j ] , 0 , a n g le l, angle2 ] + {xco r, ycor, zcor} l ; > ]; ( * Finds the two points o f the d isk c lo s e s t to the in te rs e c tin g plane; The "1" is used because i t represents th e xcoordinate; k + 1 technique works because the f i r s t p o in t o f the c ir c le is the sane th e la s t po in t; Also works because k < p o in ts , not k s p o in ts; *) F o r(k • 1, k < p o in ts, k++, { I f [ (P a rt(P a rt(lis tA , k ], 1] < xperpplane) (P a rt[ P a rt( lis tA , k + 1 ), l ] > xperpplane), < corA * P a rt[lis tA , k ]; corB » P art [lis tA , k + 1 ]; ) f { <*ELSE IF * ) I f [ (P a rt(P a rt[lis tA , k ] , 1] > xper pplane) a t (P a rt[P a rt[lis tA , k + 1 ), 1] < xperpplane), ( core « P a r t[lis tA , k ]; corD “ P a rt [lis tA , k + 1 ]; ) I; } ] > ]; C le a r[ lis tA ]; The disks are a c tu a lly represented as a polygon w ith many sides such; as 40 and a c tu a lly o n ly by it s v e rtic e s , se^nentAB and segmentCD are; e s s e n tia lly th e two sides which in te rs e c t the plane th a t is ; being observed, percentAB fin d s where those lin e segments in te rs e c t; th e plane as corE and corF; Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 98 The avoids the {0 , 0, 0} th a t th e "cor" values are o rig in a lly s e t to ; * > I f [ corA a corB is core * corD, { segmentAB * {corA, corB}; percentAB « (Abe[ (xperpplane - P a rt [Part[segmentAB, 1 ], 1])])/ (A b s [P a rt[P a rt[ segmentAB, 2 ], 1] - Part[Part[segm entAB, 1 ], 1]]); c o rl = P a rt [P a rt[ segmentAB, 1 ], 1] + percentAB* (P a rt [P a rt [segmentAB, 2 ], 1] - P art [P a rt [segmentAB, 1 ], 1]); cor2 = P a rt [ P a rt[ segmentAB, 1 ], 2] + percentAB* (P a rt [P a rt [segmentAB, 2 ], 2] - P a rt [P a rt [segmentAB, 1 ], 2]); cor3 m Part[Part[segm entAB, 1 ], 3] + percentAB*( P a rt[ P a rt[ segmentAB, 2 ], 3] - Part[Part[segm entAB, 1 ], 3]); corE = {c o r l, co r2, cor3}; segmentCD = {co re, corD}; percentCD « (Abe[ (xperpplane - Part[Part[segm entCD, 1 ], 1])])/ (Ab8[Part[Part(segm entCD, 2 ], 1] - Part[Part[segm entCD, 1 ], 1]]); cor4 » Part[Part[segm entCD, 1 ], 1] + percentCD* ( P a rt [ P art [ segmentCD, 2 ], 1] - P art [P a rt [segmentCD, 1 ], 1]); corS = Part[Part[segm entCD, 1 ], 2] + percentCD*(Part[Part[segm entCD, 2 ], 2] - P a rt[ P a rt[ segmentCD, l ] , 2]); cor6 » Part[Part[segm entCD, 1 ], 3] + percentCD* (P a rt [P a rt (segmentCD, 2 ], 3] - P art [P a rt [segmentCD, 1 ], 3]); corF ” {c o r4 , co r5, co r6}; segments? = Graphics3D[L ine[ {corE , c o rF }] ] ; lis tB = Append[ lis tB , segmentEF]; Clear[segmentAB, segmentCD, segmentEF]; corA “ corB * core « corD • {0.0, 0.0, 0.0}; corE - corF - {0.0, 0.0, 0.0}; } ] ; > ]; (* E N D OF FIRST IF L O O P *) ) ]; (* E N D OF FIRST FOR LO O P *) ( * For the g rap h ical d is p la y on th e m onitor * ) Show[ lis tB , AspectR atio -> Automatic , Axes -> Automatic , AxesLabel -> {"*, y , z} , Boxed -> True , Ticks -> {None, Autom atic, Autom atic} , PlotRange -> {{xm inp, xmaxp}, {minycor - ra d iu s , maxycor + ra d iu s }, {m inzcor - ra d iu s , maxzcor + ■ ra d iu s }} , View point -> {1 , 0, 0} ]; (* For the g ra p h ic a l d is p la y needed fo r the 2 - D box cou nting *) Reproduced with permission of the copyright owner. Further reproduction prohibited without permission. 99 p ic tu re = Show[ lis tB , ]; C le a r[ lis t B ] ; C lose[ FILENAME1 ; ] AspectRatio -> Automatic , Axes -> None (* N O T AUTOMATIC * ) , Boxed -> False ( * NOT TRUE * ) , PlotRange -> {{xm inp, xmaxp}, {minycor - ra d iu s , maxycor + ra d iu s }, {m inzcor - ra d iu s , maxzcor + ra d iu s}} , Viewpoint -> {1 , 0, 0} Reproduced with permission of the copyright owner. Further reproduction prohibited without permission.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Management of large earthquake data using the Antelope Relational Database and seismicity analysis of the 1999 Turkey earthquake sequences

PDF

Modeling of continuous tiltmeter data from the 1984 rifting event at Krafla Volcano, Iceland

PDF

Systematic high -resolution imaging of fault zone structures

PDF

Numerical simulation of whole core squeezer radon pore water profiles: Methodological considerations and evaluation of benthic fluxes and rates of bio-irrigation and advection

PDF

A tectonic model for the formation of the gridded plains on Guinevere Planitia, Venus: Implications for the thickness of the elastic lithosphere

PDF

Using small events to measure the evolution of the strain and stress fields before the 1992 Landers, California earthquake

PDF

Kinematic history and tectonic implications of the Kokoweef-Slaughterhouse fault, eastern Mojave desert, California

PDF

Numerical modeling of long period events at the Piton de la Fournaise by the use of object-oriented scientific computing

PDF

Magma mingling/mixing in a heterogeneous, multi-pulse magmatic system: An example from the Jackass Lakes pluton, central Sierra Nevada batholith

PDF

Reinterpreting the tectono-metamorphic evolution of the Tonga Formation, North Cascades: A new perspective from multiple episodes of folding and metamorphism

PDF

Paleoenvironments and paleoecology of the disaster taxon Lingula in the aftermath of the end-Permian mass extinction: Evidence from the Dinwoody Formation (Griesbachian) of southwestern Montana...

PDF

Oxygen-related biofacies in slope sediment from the Western Gulf of California, Mexico

PDF

Silicic and germanic acids: Laboratory determination of their molecular diffusivities and field study of their benthic fluxes along the California margin

PDF

Three-dimensional reconstructions of fold growth using growth strata: An example from the Catalan Coastal Ranges, northeast Spain

PDF

Paleoecology and paleoenvironments of early Triassic mass extinction biotic recovery faunas, Sinbad Limestone Member, Moenkopi Formation, south-central Utah

PDF

Oxygen isotopic evidence for fluid infiltration in the Mount Stuart batholith, Washington

PDF

A proxy for reconstructing histories of carbon oxidation in the Northeast Pacific using the carbon isotopic composition of benthic foraminifera

PDF

Paleoseismology and geomorphology of the eastern Sierra Madre fault: Evidence for a >8ka age of the most recent event

PDF

A unified methodology for seismic waveform analysis and inversion

PDF

Structural evolution of the southwestern Daqing Shan, Yinshan Belt, Inner Mongolia, China

Asset Metadata

Creator Libicki, Eric Jason (author)

Core Title Simulating nonplanar fault surfaces using stochastic branching

Degree Master of Science

Degree Program Geological Sciences

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag geophysics,OAI-PMH Harvest

Language English

Contributor Digitized by ProQuest (provenance)

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c16-44100

Unique identifier UC11341104

Identifier 1409597.pdf (filename),usctheses-c16-44100 (legacy record id)

Legacy Identifier 1409597.pdf

Dmrecord 44100

Document Type Thesis

Rights Libicki, Eric Jason

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the au...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus, Los Angeles, California 90089, USA