S T U D Y OF T H E S C A L I N G A N D T E M P O R A L PROPERTIES OF A SIMPLIFIED E A R T H Q U A K E MODEL By Daniel Groleau B. Sc. (Physics) Universite du Quebec a Rimouski M . Sc. (Physics) Universite de Sherbrooke A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F THE REQUIREMENTS FOR T H E DEGREE OF D O C T O R OF PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES PHYSICS A N D A S T R O N O M Y We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH C O L U M B I A August 1997 © Daniel Groleau, 1997 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Physics and Astronomy The University of British Columbia 6224 Agricultural Road Vancouver, B.C., Canada V 6 T 1Z1 Date: - Abstract Certain driven systems consisting of a large number of elements evolve towards a critical state with no characteristic time or length scales. This class of phenomena is described as Self-Organized Criticality (SOC). SOC relies on the condition of a slow driving of the systems and the existence of fast burst-like responses of them and includes earthquakes. We employ a model proposed by X u et al. for ruptures in an elastic medium, subject to shear stress, and apply it to the study of earthquakes. In the model, the size of an earthquake is defined as the number of ruptures occurring sequentially on the basic units of discretization (squares) of the medium. A histogram of the earthquake sizes shows that the model is not completely scale invariant due to finite-size effects. To take them into account, we implement a finite-size-scaling analysis. The results of this analysis show that the model is scale invariant only when there is stress conservation. So, the model displays SOC in the conservative case only. We also study the dynamic quantities of the model, in particular the average stress in the system. The sets of average stress values are analyzed using two types of time series analysis. The nonlinear forecasting analysis investigates whether time series exhibit low-dimensional chaotic behavior as opposed to high-dimensional (or stochastic) behavior. We find that the above time series have a nonlinear structure, but with a substantial stochastic component, so SOC is inherently high-dimensional. The appearance of nonlinear structure is due to the fact that the system stops following linearly the external drive when it releases stress through earthquakes. The rescaled range analysis characterizes the time correlations (or memory effects) in the time series. We find strong positive time correlations in the above time series. Their presence is due to the nature of the driving ii in the model. These memory effects are destroyed as soon as a large earthquake resets the system. iii Table of Contents Abstract ii List of Tables vii List of Figures viii Acknowledgements 1 2 xvii Introduction 1 1.1 Self-Organized Criticality 1 1.2 Identification of Self-Organized Criticality 3 1.3 Geophysics as a Field of Application for Self-Organized Criticality . . . . 4 1.4 Interests of the Concept of Self-Organized Criticality for Physicists 5 1.5 Organization of the Thesis ... 6 Lattice M o d e l for the S t u d y of Earthquakes 8 2.1 Approaches to Fracture Modeling 8 2.2 Lattice Models 10 2.2.1 Lattice models considered by engineers 11 2.2.2 Lattice models considered by physicists 12 2.3 2.4 A New Lattice Model 15 2.3.1 A new discretization scheme 15 2.3.2 Advantages of that discretization scheme over others 20 Application to the Study of Earthquakes iv 21 2.5 3 4 Comparison of the Model of X u et al. to Other Models 27 2.5.1 Model of X u et al. versus other earthquake models 27 2.5.2 Model of X u et al. versus models of integrate-and-fire oscillators . 30 Study of the Scaling Properties of the M o d e l 32 3.1 Motivation for this Study 32 3.2 Finite-Size Scaling 33 3.3 Results of the Finite-Size-Scaling Analysis 35 3.3.1 Results in the conservative case (a = 1) 36 3.3.2 Results in the non-conservative case (a < 1) 46 3.4 Stress Ordering in One-Dimensional Lattices 50 3.5 Mean-Field Version of the Model of X u et al 54 3.6 Summary of the Important Results and Conclusions 57 T i m e Series Analyses 62 4.1 T i m e Series of Interest 62 4.2 Nonlinear Forecasting Analysis 69 4.2.1 Introductory discussion 69 4.2.2 State space reconstruction 76 4.2.3 Nonlinear function approximation 78 4.2.4 Nonlinear forecasting algorithm considered 80 4.2.5 Tests of the algorithm with time series having a known structure . 83 4.2.6 Results obtained with the time series of interest 89 4.3 Rescaled Range Analysis 103 4.3.1 Introductory discussion 103 4.3.2 Fractal time series analyses 112 4.3.3 Rescaled range analysis algorithm considered 113 v 4.4 5 4.3.4 Tests of the algorithm with time series having a known value of H 116 4.3.5 Results obtained with the time series of interest 117 Summary of the Important Results and Conclusions Conclusion 126 129 Bibliography 137 Appendices 142 A Lattice G r e e n Function Simplification and Properties 143 A.l 143 B Simplification of the Lattice Green Function A.2 Properties of the Lattice Green Function 147 C o a r s e - G r a i n i n g of the Lattice G r e e n Function 149 vi List of Tables Exponent r and characteristic earthquake size So for Lxl a , X = 0, 9i = 0.25 and three values of a th vii lattices, random List of Figures 2.1 One-layer planar region considered in the model of X u et al. T h e region is discretized into squares. It is embedded in an infinite medium (providing the necessary external drive) and the effect of any activity originating outside the region is neglected (open boundary conditions). In the model, the interaction between the external drive (represented by the arrows outside the region) and a square of the region is treated in a mean-field way, i.e. it is the same whatever the square location. T h e arrows at the corners of one of the square inside the region are a schematic representation of a double couple force redistribution following the shear rupture of that square. T h i s figure was modified from [114] 3.1 18 Log-log plot of the cumulative size-frequency distributions for L x 1 lattices, random ath, X = 0, 9i = 0.25 and a = 1. For comparison, we show the curve C ( S ) = S " 0 3.2 37 30 Log-linear scaling plots of the cumulative size-frequency distributions for I x 1 and L x 2 lattices, random a , th X = 0, Q\ = A l l curves use the fitting parameter r = 0.30. 0.25 and a — 1. For the left and middle curves, all ruptures are counted while for the right curves, only ruptures on distinct squares are counted. T h e left curves were obtained for L x 2 lattices using u = 1.38, the middle curves for L x 1 lattices using v = and the right curves for L x 1 lattices using v = 1.01 viii 1.38, 38 3.3 Log-log plot of the cumulative size-frequency distributions for L x 1 lattices, constant a , th X = 0, B\ = 0.25 and a = 1. T h e segment S < 50 is a crossover region and its origin is explained in section 3.4. Beyond S > 50, C(S, L) has a linear segment with a slope close to 0.30. T h i s linear segment gets longer as L increases but cuts off at large S. For comparison, we show a curve proportional to S 3.4 1-0 - 40 30 Log-linear scaling plots of the cumulative size-frequency distributions for L x 1 lattices, random a h, ®i — 0.25, a = 1 and two values of X. t Exponents used are r = 0.30 and v = 1.38 3.5 41 Log-linear scaling plot of the cumulative size-frequency distributions for a 1000 x 1 lattice, random a h, Q\ = 0.25, a = 1 and three values of t X. Exponents used are r = 0.30 and v = 1.38 3.6 42 Log-linear scaling plots of the cumulative size-frequency distributions for L x 1 lattices, random a h, X = 0, a = 1 and two values of B\. Exponents t used are r = 0.30 and v = 1.38 3.7 44 Log-linear scaling plot of the cumulative size-frequency distributions for a 1000 x 1 lattice, random cr^, X — 0, a = 1 and three values of #/. Exponents used are r = 0.30 and v = 1.38 3.8 45 Log-log plot of the cumulative size-frequency distributions for a 2000 x 1 lattice, random cr , X = 0, 9[ = 0.25 and three values of a. T h e dashed tfl lines are the fits to relation (3.3) 3.9 48 Log-log plot of the cumulative size-frequency distributions for L x 1 lattices, random a , X = 0, 0, = 0.25 and a = 0.94 49 th 3.10 Distribution of the stress differences Si (i = 1,2, ...,999) between the nearest-neighbor elements of a 1000 x 1 system with constant a h, X = 0, t 9 = 0.25 and a = 1 52 t ix 3.11 Histogram plots of the stress differences 8i between the nearest-neighbor elements of a 1000 x 1 system with constant o~ h, X = 0, 9i = 0.25 and t a = 1 at two distinct times. T h e solid line curve was obtained with data taken at the arbitrary time t 0 (data of figure 3.10) and the curve plotted using a dashed line and lozenges was obtained with data taken at time t\ > to- Fifty bins were used to get the histograms 3.12 Log-log plot of 1 — a ff e 53 values versus L for L x 1 lattices, random a h, t X = 0, 9i = 0.25 and a = 1 obtained with the original model of X u et al. T h e slope of the straight segment was found to be —1.02 58 3.13 Log-log plot of the cumulative size-frequency distributions for systems of L elements equally coupled with each other. T h e values of the parameters are X = 0 and the stress thresholds (i.e. the a h) are taken from a uniform t distribution with 9\ = 0.25. From left to right, the distributions refer to systems having L = 200, L — 400 and L = 1000. T h e values of a differ for all three distributions (see how they were adjusted in the text) 4.1 59 T y p i c a l time series of type 1: Size of an earthquake sequences (Si) as a function of the occurrence time (ti). T h e whole time series has 5 x 10 4 data but for clarity, only the first 2000 are plotted. T h e sequences were generated from a 100 x 100 system with random a ht T h e time at which the first sequence happens is arbitrarily set to t\ = 0 4.2 65 T y p i c a l time series of type 1 [variety (b)]: T i m e intervals ti — 2, 3 , 2 0 0 0 ) between consecutive earthquake sequences. (i = T h e sequences were generated from a 100 x 100 system with random o h (see figure 4.1). t x 66 4.3 T y p i c a l time series of type 2: Average stress (a) as a function of the long time (t). T h e data were generated from a 100 x 100 system with random a . th T h e whole time series has 5 x 10 data but for clarity, only the first 4 500 data (i.e. the end points of the vertical lines) are plotted. T h e time of the first datum is arbitrarily set to zero 4.4 68 Autocorrelation function for the time series of type 1 [variety (b)] plotted on figure 4.2. For convenience, the value r = 1 is not shown on the graph. 0 4.5 Autocorrelation function for N — 2000 data of the time series of type 2 plotted partially on figure 4.3 4.6 70 71 Normalized R M S forecasting error E (k) m as a function of k for time series having N = 1000 data generated by the logistic map (with fj, = 4.0) and the Henon map (with A — 1.4 and B = 0.3). T h e solid line is for the logistic map time series and the dashed line for the Henon map time series. T h e parameters of the nonlinear forecasting algorithm are Nf = 800 and m = 2 4.7 86 Normalized R M S forecasting error E (k) m as a function of k for a time series having N = 1000 data generated by the T A R model (with a = —0.4 and (5 = 0.5). T h e parameters of the nonlinear forecasting algorithm are Nf = 800 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line). 4.8 Normalized R M S forecasting error E (k) m 87 as a function of k for a time series having N = 1000 data generated by the A R ( 1 ) model (with 0 = 0.5). T h e parameters of the algorithm are Nf — 800 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line) 88 xi 4.9 Normalized R M S forecasting error E (k) m as a function of k for a time series of type 1 [variety (c)] for a 100 x 100 system with random o . ih The sampling time interval At is equal to 10 times the average time interval between the earthquake sequences. T h e parameters of the nonlinear forecasting algorithm are Nf = 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line) 91 4.10 Normalized R M S forecasting error E (k) m as a function of k for a time series of type 1 [variety (a)] for a 1000 x 1 system with constant a h- T h e data in t the time series (N = 331) are the time intervals between large earthquakes (S > 1000). T h e parameters of the nonlinear forecasting algorithm are Nf — 264 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line). 92 4.11 Histogram plot of the time intervals between large earthquake sequences (S > 1000) generated in a 1000 x 1 system with constant a .tn T h e data set has 331 data and 20 bins were used to get the histogram. T h e vertical axis is logarithmic in order to better illustrate that the distribution is exponential 93 4.12 Normalized R M S forecasting error E (k) m as a function of k for a time series of type 2 for a 100 x 100 system with random a h- T h e parameters t of the nonlinear forecasting algorithm are Nf = 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line) 4.13 Normalized R M S forecasting error E (k) m 95 as a function of k for a first- differenced time series of type 2 for a 100 x 100 system with random a ht T h e parameters of the nonlinear forecasting algorithm are Nf = 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line) xii 96 4.14 Normalized R M S forecasting error E (k) m as a function of k for a time series of type 2 for a 100 x 100 system with constant a h- T h e parameters t of the nonlinear forecasting algorithm are Nf = 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line) 4.15 Normalized R M S forecasting error E (k) m 97 as a function of k for a time series of type 2 for a 100 x 100 system with random a h when the parameter a t controlling the level of stress conservation in the system (see chapter 3) is set to 0.9. T h e parameters of the nonlinear forecasting algorithm are Nf = 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line). 98 4.16 Normalized R M S forecasting error E (k) m as a function of k for an evenly sampled time series of type 2 for a 100 x 100 system with random a h- T h e t sampling time interval A t is about equal to 1.5 times the average time interval between the earthquake sequences. T h e parameters of the nonlinear forecasting algorithm are Nf = 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line) 101 4.17 Normalized R M S forecasting error E (k) m as a function of A; for a time series which has the same power spectrum as the time series of figure 4.16, but randomized phases. T h e parameters of the nonlinear forecasting algorithm are Nf = 1600 and m — 2 (solid line), m = 3 (dashed line), m = 5 (dotted line) 102 4.18 Fractional Brownian motion (fBm) signals having 8193 data generated by the successive random addition method. From top to bottom, the curves have H = 0.1, H = 0.5 and H = 0.9 and have been shifted for clarity. . . xiii 107 4.19 R u n n i n g sum of an evenly sampled time series of type 2 for a 100 x 100 system with random o h having N t interval At — 16, 384 data. T h e sampling time is equal to 1.5 times the average time interval between the earthquake sequences. T h e trend and the mean have been removed from the time series before computing the running sum 110 4.20 Periodogram estimate of the power spectrum of the running sum of an evenly sampled time series of type 2 for a 100 x 100 system with random Oth- T h e sampling time interval of the time series A t is equal to 1.5 times the average time interval between the earthquake sequences, where the time series has 16, 384 data. We have removed the trend and the mean before calculating the power spectrum. T o reduce the fluctuations, we have computed an average power spectrum over 8 non-overlapping subsets (with 2048 data) of the time series resulting from the running sum operation. T h e dashed lines were obtained by doing linear regressions and have slope 0 = 2.94 ± 0.05 (left segment) and (3 = 3.46 ± 0.03 (right segment). 4.21 Log-log plot of (R/S)M N = 2 1 3 . . . Ill versus M for three synthetic time series having data generated by the successive random addition method. From top to bottom, the curves correspond to the cases H = 0.9, H = 0.5 and H = 0.1. T h e dashed line at the right of the figure has a slope of 0.1 and is there to illustrate that the curve for H = 0.1 will only slowly reach the slope value H RR = 0.1 as M increases. Note that each curve was obtained by averaging (R/S)M (for a given M) over 5 realizations corresponding to the same value of H. For each curve, we have excluded the point corresponding to M = 2 since it is trivial and always has (R/S)M=2 xiv = 1- 118 4.22 Log-log plot of (R/S)M versus M for an evenly sampled time series of type 2 for a 100 x 100 system with random a t h having A t equal to 1.5 times the average time interval between the earthquake sequences (N = 2 ) . 14 The curve corresponding to this time series is plotted using a solid line and lozenges. M T h e dotted line (H = 0.9) is a fit of the linear segment G [2 ,2 ] of the latter curve. T h e curve at the bottom was obtained 2 7 with 5 random shuffles of the time series: the curve of (R/S)M versus M was calculated for each realization and an average value of (R/S)M over the 5 realizations was computed for each M 4.23 Log-log plot of (R/S)M 121 versus M for evenly sampled time series of type 2 for L x L (L = 25,50,100) systems with random a . t h T h e sampling time interval A t is equal to 1.5 times the average time interval between the earthquake sequences and the number of data is N = 2 1 4 for all three time series. From top to bottom, the curves refer to the 100 x 100, 50 x 50 and 25 x 25 systems respectively 4.24 Log-log plot of (R/S)M 122 versus M for an evenly sampled time series of type 2 for a 100 x 100 system with constant a h having A t equal to 1.5 times t the average time interval between the earthquake sequences (N = 2 ) . 14 The curve corresponding to this time series is plotted using a solid line and lozenges. M T h e dotted line is a fit (H = 0.9) of the linear segment G [2 ,2 ] of the latter curve. T h e curve at the bottom was obtained 2 7 with 5 random shuffles of the time series: the curve of (R/S)M versus M was calculated for each realization and an average value of (R/S)M the 5 realizations was computed for each M xv O V E R 123 4.25 Log-log plot of (R/S)M versus M for an evenly sampled time series of type 2 for a 100 x 100 system with random a h- T h e sampling time interval A t t is equal to 1.5 times the average time interval between the earthquake sequences and the number of data is N — 2 . T h e curve plotted using a 14 solid line and lozenges was obtained with the original time series whereas the one plotted using a dashed line and pluses was obtained with 5 surrogate time series having the same power spectrum as the original time series but randomized phases. For the latter curve, (R/S)M was calculated for each realization, for a given M, and an average value of (R/S)M 5 realizations was computed for this value of M xvi over the 124 Acknowledgements I would first like to thank my supervisor, Dr. Birger Bergersen, for his help and financial support over the past five years. His enthusiasm has helped me finding the energy to finish my Ph.D. project, especially in the difficult times. I would also like to thank Dr. Huang-Jian X u who initiated this work by introducing the earthquake model. I am particularly grateful to him for having helped me in the early stages. A special thank to Dr. Lawrence Ward of the Department of Psychology at U B C who allowed me to use the code of his nonlinear forecasting algorithm. In addition, I benefited a lot from my discussions with Dr. Ward. I also thank him for having accepted to sit as a member of my Ph.D. committee in the last year. Thanks also to my other Ph.D. committee members, Dr. Rosemary Knight of the Department of Earth and Ocean Sciences and Drs. Robert Parsons and Nathan Weiss of the Department of Physics and Astronomy, who helped me in different ways. I gratefully acknowledge the financial support, through graduate fellowships, of Le Fonds FCAR of the government of Quebec as well as of the National Science and Engineering Research Council of Canada (NSERC). Merci a mes parents, Louisette et Georges, pour m'avoir toujours supporte moralement durant mes annees d'etudes. Finalement, sans le support de ma compagne de vie, Renee, je n'aurais probablement jamais termine ce doctorat. Je la remercie d'avoir toujours ete la lorsque j'avais besoin de reconfort et d'encouragement. Sa grande comprehension m'a ete indispensable. Je ne pourrais terminer ces remerciements sans mentionner ma fille, Carmen, qui est nee durant la derniere annee de mon doctorat, plus precisement en pleine redaction de cette these. Son arrivee dans la vie m'a donne l'energie necessaire pour terminer mon doctorat. xvii Chapter 1 Introduction Once upon a time, science was not divided into disciplines as it is today. There was no distinction between social and natural sciences, and certainly not the fragmentation of the latter into physics, chemistry, biology, geology, etc. However, when one is confronted with a new problem or phenomenon, it is usually impossible to neatly assign the problem to a particular discipline. One needs to recall and fuse together the knowledge one learned before and, if that is insufficient, to consult scientists in other areas [56]. In recent years, under the banner of complex systems (or complex phenomena to emphasize the fact that the interest is in problems from the real world), some scientists have ventured out of their trained disciplines. Complex systems research is a synthetic approach in the investigation of systems consisting of a large number of interacting, simple elements, irrespective of their origins (see [56] for a general outline of the complex systems research). It has led recently to the introduction of the concept of Self-Organized Criticality ( S O C ) by Bak, Tang and Wiesenfeld [3]. T h i s concept has been introduced in order to understand the origin of the power laws observed in the study of various phenomena [81, 37, 57, 70, 61, 28]. 1.1 Self-Organized Criticality According to S O C , certain dynamically driven systems consisting of a large number of elements evolve spontaneously towards a critical globally stationary dynamical state with no characteristic time or length scales. T h i s critical state is reached without fine tuning 1 Chapter 1. Introduction 2 of an external control parameter, unlike for phase transitions in equilibrium statistical physics. Those systems are made critical by the choice of a threshold (nonlinear) dynamics that forbids them to follow adiabatically the external drive. T o illustrate the basic ideas of S O C , Bak, T a n g and Wiesenfeld [3] used a cellular automaton inspired by the creation of avalanches in a pile of sand. In this model, sand is added grain by grain at a random location in the lattice (random driving) until somewhere a local slope becomes unstable (i.e. the local slope is larger than a threshold value), and an avalanche is initiated. T h e pile eventually reaches a stationary critical state, characterized by a critical slope, in which additional grains will fall off the pile via avalanches of all sizes, distributed in size and lifetime according to a power law. It is of interest to note that the relative number of large and small avalanches in the model does not depend on the microscopic mechanisms, i.e. the global behavior of the sandpile can be synthesized from a simple model such as the one above. Over the years, three facts have been established about S O C [94]. First, some systems qualify as self-organized critical if their large scale evolution obeys a diffusion equation which satisfies a global conservation law [49, 41]. Second, some systems, which exhibit diffusion-like response but do not obey a global conservation law, nevertheless exhibit S O C [77, 26, 1]. In these systems, S O C is believed to be related to the problem of synchronization of coupled threshold oscillators of relaxation [67, 7, 19]. T h i r d , a feedback mechanism must operate which describes the action of the order parameter onto the control parameter [92] and attracts the dynamics to a critical state. Recently, a unifying formalism has been presented by G i l and Sornette [36] to describe the general conditions under which a system exhibits S O C . T h e y have introduced a continuous formulation in terms of the Landau-Ginzburg theory of phase transitions which takes into account the nonperturbative nature of the threshold dynamics. S O C was identified as the regime where diffusive relaxation is faster than the instability growth Chapter 1. Introduction 3 rate whereas synchronization of threshold oscillators was observed in the limit of slow diffusion. Now that the S O C concept has been described, we can talk about the properties that are necessary for a model to display S O C . 1.2 Identification of Self-Organized C r i t i c a l i t y It is often asserted that the hallmark of S O C is the simultaneous existence of two properties [93]: (1) power law distribution of events, and (2) spatial and temporal correlation functions decaying algebraically. According to [93], property (1) is more characteris- tic of the class of phenomena described by S O C , since it relies on the condition of a very slow driving of the system and the existence of fast burst-like responses. A num- ber of self-organized critical phenomena can be described as ordinary phase transitions made self-organized [94]. Apart from the absence of tuning and the property of selforganization, the latter are similar to ordinary phase transitions. We will not consider this class of phenomena in this thesis. For a system to be self-organized critical, it must exhibit a power-law distribution of event (response) sizes. B u t the opposite is not true: the observation of a power-law distribution of event sizes does not necessary imply that the system is self-organized critical [93]. A s an illustration of this, it is of interest to note that a number of numerical and experimental works claiming the observation of S O C in fact rely on a different physical mechanism, called sweeping of an instability [93]. T h i s mechanism involves the slow sweeping of a control parameter on the approach to and possibly beyond an instability at the critical value of that control parameter. In other words, the power law arises if a control parameter crosses over from one side of the critical point to the other side. A s a Chapter 1. Introduction 4 result, the system can not operate persistently since the control parameter does not remain constant. M a n y examples of how that mechanism operates are provided in [93] and, among them, foreshocks (shocks preceeding a large earthquake) and acoustic emissions during the progressive deformation of materials submitted to increasing stresses. T h e S O C concept applies to many different phenomena occurring in nature and in particular, to phenomena in the field of geophysics. 1.3 Geophysics as a F i e l d of A p p l i c a t i o n for Self-Organized C r i t i c a l i t y Geophysics is a field in which many phenomena are well suited to the concept of S O C . In particular, the earthquake phenomenology was argued by several authors [96, 2] to bear the signature of S O C . T o illustrate this, we summarize a part of the discussion in the paper of Sornette et al. [94]. Let us consider a model elastic tectonic plate [20] scaled down in the laboratory to perform a mechanical experiment. A shear force (F) is, in the simplest case, imposed on two opposite borders of the plate, the other two borders being free. A s F increases, the plate starts to deform. For low F, after some transient during which the plate deforms and adjusts itself to the applied force, the plate becomes static: the strain rate (r) is zero (if any ductile behavior is neglected). A s F increases, the transient becomes longer and larger deformations develop within the plate. A plasticity threshold (F ) will c be eventually reached and then the plate will start to flow with a non-zero strain rate (r > 0) under fixed F. A s F increases above F , r increases. In laboratory experiments, c such a transition from a brittle to a ductile behavior may be critical. Parameters F and r qualify respectively as the control parameter and the order parameter (r is non-zero only when F > F ). c Now, instead of controlling the shear force exerted on the plate, we imagine that Chapter 1. Introduction a constant 5 (very small) shear rate of deformation (r —» 0) is applied at the borders of the plate. T h e n , the natural boundary condition for earthquakes and plate tectonic deformations is recovered. Such a situation has been studied in many works showing the existence of S O C , in particular in the shape of a power-law earthquake size distribution (the well-known Gutenberg-Richter power law [43]). T h e above boundary condition is nothing but driving the plate at the critical point (F —¥ F ) c by controlling the order parameter (r) to a very small value (r —>• 0), thus ensuring the critical properties of the plate. T h e fact that r is non-zero even if it is very small proves that the plate is an out-of-equilibrium system. T h e above example illustrates the relevance of the concept of S O C to the earthquake phenomenology. We will focus on the earthquake phenomenon throughout the present work. T h e S O C concept is interesting not only from a practical point of view (e.g. understanding earthquakes), but also from a theoretical point of view. 1.4 Interests of the Concept of Self-Organized C r i t i c a l i t y for Physicists T h e interests of the physics community in the study of S O C is twofold. Firstly, S O C , as was shown by Sornette et al. [94], can be mapped onto ordinary criticality in many cases. T h e n , rather than searching for the underlying mechanism for S O C , physicists can study unstable critical points (i.e. critical points which are not self-organized). T h i s task can be performed by using the toolbox developed in the last decades in this field. However, setting the order parameter of an unstable critical point to an infinitesimal number does not allow for a standard renormalization group procedure, because there is no control parameter and the exponents related to the distance from the critical point do not exist. There have been some attempts [76, 49, 83, 24] at finding a general renormalization group Chapter 1. Introduction 6 theory for systems right at their critical point, but much work along this line remains to be done though. Such a theory could provide a new class of tools for critical phenomena as suggested in [94]. Secondly, S O C applies to phenomena actually occurring in nature and so may be useful to their understanding. A m o n g these phenomena are, apart from earthquakes, formation of river networks [100], formation of clouds [72], evolution of populations [1], neuron dynamics [7, 19], etc. T h e models used to study such phenomena are often idealized in the sense that they capture only the essential features of real systems. If these simple models would turn out to be robust with respect to various modifications, then extrapolation to real situations could be possible. (It is of interest to note that such an approach has been successful in equilibrium statistical mechanics, where universal features in real systems can be understood from the study of simple models). Also, if models used to study a priori not closely related phenomena were to have similar features, then this could prove that a unique universal dynamics is behind those phenomena. Indeed, why could we not learn about the dynamics of neurons of the human body by studying earthquakes? Our interest in the concept of S O C is oriented towards the understanding of the earthquake phenomenon in the current work. 1.5 Organization of the Thesis T h e present work is devoted to the detailed study of a model introduced by X u et al. [113, 114] which was intended to be applied to earthquakes. T h i s model was shown to be self-organized critical [114]. However, we believe that a more rigorous proof is needed before to state whether or not the model displays S O C . In addition, the temporal aspects of the model have almost entirely been left out in [113, 114] and so, we want to study Chapter 1. Introduction 7 them thoroughly in this work. T h e organization of the present work goes like this: in chapter 2, we describe the approach used to model fractures in a planar region, the occurrence of earthquakes being closely related to the appearance of fractures in nature. We use a particular lattice model, which has been introduced recently by X u et al., to do so. In chapter 3, we study the scaling properties of the model by means of a finite-size-scaling analysis, a technique which has been extensively used in the study of equilibrium phase transitions. T h i s will enable us to clearly state whether the model displays S O C . In chapter 4, we implement two types of time series analysis, namely a nonlinear forecasting analysis and a rescaled range analysis. T h e time series are generated from the model. Finally, in chapter 5, we summarize our results, draw the relevant conclusions and make suggestions for future work. Chapter 2 Lattice M o d e l for the S t u d y of Earthquakes In nature, the occurrence of earthquakes is closely related to the appearance of fractures. T h e modeling of fractures can be done using different approaches (see chapter 5 of [45] for instance). These approaches are discussed in section 2.1. T h e approach presented in section 2.2 follows an intermediate route and uses lattice models as the tool to study fractures in the presence of disorder. T h e lattice model that we use in this work is presented in section 2.3 and applied to the study of earthquakes in section 2.4. Finally, it is compared to other earthquake models and to models of integrate-and-fire oscillators in section 2.5. 2.1 Approaches to Fracture M o d e l i n g T h e description of fracture of samples of materials in the cases where disorder is important is the main goal of the material scientist interested in the fracture phenomenon. One possible approach he (or she) can use is to set up some microscopic rules which determine the mechanical behavior of the medium as well as how rupture will take place. If the medium is homogeneous and the rupture process deterministic, one can very simply know what the macroscopic behavior will be. However, in the presence of disorder, the solution to the problem is not that easy. A m o n g the various ways of handling disorder directly, three different approaches can be distinguished [45]: mean-field, local models and renormalization methods. discuss each of them. 8 Next, we Chapter 2. Lattice Model for the Study of Earthquakes 9 The mean-field approach treats the surroundings of a volume element as a homogeneous embedding medium. In this way, each element interacts equally with all the others. However, because fracture is sensitive to disorder, this approach gives poor results. The approach using local models tries to reproduce as closely as possible the real microscopic behavior starting from first principles. This is the case, for instance, for molecular dynamics that we now briefly describe. Suppose one chooses an interatomic potential, for example the Lennard-Jones (6-12) potential, to decribe the interactions between the particles of the system under study. The particles are placed on a regular lattice. On the outer boundary of the system, an external load is applied by assigning a constant force vector to the particles lying on this boundary. To obtain the motion of the system, time is discretized into slices of duration At and one calculates the position that each particle has after the time At. In the simplest case, one calculates, using the equations of motion, the velocity of each particle. When all velocities v are known, all particles are simultaneously displaced by vAt. Then, the step is repeated for the next time interval starting from the new positions, and so on. The molecular dynamic framework provides naturally the fracture properties without having to insert by hand any rheological behavior. However, the study of a system with several cracks in a disordered medium is, with the present means of computation, not possible. Much work remains to be done along this line though. The approach using renormalization methods views the fracture process as a critical phenomenon. It expresses the invariance of some physical quantities under a change of length scale. However, the few situations where this technique has been used contain some severe approximations, or use unphysical and simplistic geometries [91, 74]. The above approaches are not adequate to study in an effective way the fracture phenomenon. This pessimistic picture seems to close all doors towards a solution. Nevertheless, the situation is not so hopeless. The approach used in [45] follows an intermediate Chapter 2. Lattice Model for the Study of Earthquakes 10 route: it describes very simply the mechanical behavior of mesoscopic (not microscopic) elements and considers the collective behavior of these elements in the presence of disorder. We now summarize this approach. T h e interested reader is referred to [45] for the details. 2.2 Lattice M o d e l s A numerically tractable lattice formulation of fracture, alternative to molecular dynamics, is given by discretizations of continuum equations. In this case, the medium is reduced to a set of points embedded into a grid. O n l y local laws (e.g. the balance of force) are considered and their implementation involves for each point only a few neighbors. T h e calculation of collective properties like the equilibrium displacement field is reduced to solving a set of coupled linear equations. Lattice models have their validity at much larger length scales (than the atomistic length scale) where the medium can be described by continuous fields. T h e breaking of the lattice is not a natural consequence of the simulation but has to be put into the model by hand as an additional rule of the model. T h i s rule can lean on experimental data or on phenomenological laws, but the fact that it does not follow from first principles shows that lattice models are less fundamental than molecular dynamic models. One advantage of lattice models is that they allow for the easy introduction of disorder. A n even bigger advantage is that they are more rigid than a crystal in molecular dynamic models because the number of neighbors is fixed and the breaking of a bond is an irreversible yes-no decision. Lattice models are therefore good candidates to.overcome the numerical difficulties encountered in molecular dynamic models. Lattice models have been studied by engineers and physicists, but the methods used in the two science communities are different. Indeed, engineers use mostly Finite Element Chapter 2. Lattice Model for the Study of Earthquakes 11 Models ( F E M ) whereas physicists focus on less formal models. We now discuss in some details lattice models as considered by both engineers and physicists. 2.2.1 Lattice models considered b y engineers In engineering practice, the most commonly used lattice models are F E M . We want here to give a brief outline of the finite element method as applied to fracture problems. T h e reader is referred to one of the many pedagogical textbooks on the subject for more details (see for instance [105]). Basically, the F E M are ways of discretizing the general continuum rheological equations on a lattice (usually triangular in two dimensions). T h e displacement field is ap- proximated by a best fit within a subspace of trial functions. T h e best fit is found by minimizing the total potential energy over all admissible fields. In the F E M , the description of a real elastic solid will be accurate only if the displacement field varies slowly over the size of the elements used. In particular, if one wants to introduce a crack in a medium, the displacement field around the crack will be accurate only if the elements are much smaller than the size of the crack itself. T h i s is a serious drawback for the modeling of a dense array of microcracks in a medium. Similar drawbacks, although less dramatic, will appear if one wants to model a heterogeneous solid. In fact, in such a case, the size of the heterogeneities should be large compared to the mesh size. Also, fluctuations of elastic moduli from element to element should be small. FEM are popular because of their simplicity (although elaborate codes often use all sorts of sophisticated optimizations) and flexibility in dealing with very different problems. These codes are tools and used like black boxes rather than being considered to have a physical interest on their own. T h e lattice models used by physicists are less sophisticated than F E M . However, they have the advantage of being built in such a way Chapter 2. Lattice Model for the Study of Earthquakes 12 as to have all the steps directly accessible and easily interpreted physically. 2.2.2 Lattice models considered by physicists Models considered by physicists tend to have a certain common setting which can be summarized in the following way [45]: the medium is discretized such that all spatial sites are equivalent (contrary to F E M , the grid is not made finer in regions of higher stress) and have the same number of neighbors (regular lattice). T h e variables that characterize the medium are placed on the sites of the lattice. T h e equations that describe the medium are discretized so that, for each site, one has one equation per variable which only involves variables on the neighbor sites. Also, a boundary condition on the outer perimeter is implemented. T h e boundary condition on the internal boundary (crack surface) reflects the response of the bonds that constitute the crack. T h e n , the set of equations has a unique solution. A s mentioned in [45], the simulation of a rupture process must be done in an iterative way: the equations must be solved in order to determine which bonds should be broken, but once the bonds have been broken, the internal boundary condition gives rise to a change in the internal stress. Consequently, the equations must be solved again if one wants to know which bond to break next, and so on. T h e process is computationally demanding because the set of linear equations must be solved as many times as bonds are broken (unless several bonds are broken simultaneously). T h e algorithm performed in one iteration can be decomposed into four steps [45]: (1) solve the set of equations, (2) determine the set of bonds eligible to be broken, (3) select from the eligible bonds one which will be broken, Chapter 2. Lattice Model for the Study of Earthquakes 13 (4) break the bond, i.e. change its elastic properties. Each of these steps allows for a large variety of options that describe many possible physical situations. Step (1) describes the nature of the medium and the externally applied constraints, step (2) the connectivity of the crack, step (3) the breaking rule and disorder and step (4) distinguishes for instance between a breakdown and a fuse problem (see [45] for the details). We now briefly discuss each of these steps. To do step (1), the type of model and the type of external constraints need to be specified. For the elasticity, scalar, vectorial and tensorial models are found in the literature. T o our knowledge, there exists only one scalar model dealing with elasticity and it was introduced in [68]. T h i s model simulates antiplane shear deformation along the axis O z of a thin plate placed in the x-y plane. A constant slow antiplane velocity along Oz is applied at the boundaries of the plate. Vectorial and tensorial models allow for a larger variety of externally imposed constraints than scalar models. Vectorial models of common use are the central-force model [32], the bond-bending model [52] and the beam model [86]. T h e latter vary in the discretization scheme used to solve the equations of elasticity. T h e central-force model, the simplest of the three, is a network of Hookean springs which can freely rotate around the sites of the lattice. So, in this model, force can act only in the direction of the bond. T h e other two vectorial models fall in the same class in the sense that a bond can not rotate around a site without a cost in energy. Tensorial models have been introduced only recently by Chen et al. [15] and X u et al. [113, 114] and have some -advantages over vectorial models, in particular they are computationally less demanding. T h e discretization scheme used in [113, 114] will be presented in the next section and compared to the one used in [15] in section 2.5. T o study crack propagation [step (2)], one can pursue two philosophies: understand how one crack grows or see how in a stressed medium cracks appear, coalesce and break 14 Chapter 2. Lattice Model for the Study of Earthquakes the system. T h e first philosophy is artificial in most experimental situations since the surface of a crack in a lattice is not uniquely defined. T h e second philosophy corresponds to the majority of experimental situations because all not yet broken bonds are eligible to be broken. T h e breaking rule and disorder are two major ingredients of the model [step (3)]. T h e breaking rule must be put by hand into the model to reflect the phenomenological nature of lattice models. It is usually chosen as simple as possible. Disorder has a major impact on the breaking rule as well as on the breaking process. If disorder is unaffected by the breaking process, then it is quenched and otherwise it is annealed (interplay of disorder and fracture has to be taken into account). O n the length scales that lattice models are supposed to describe, disorder can be reduced to a spatial dependence of the elastic modulus a n d / o r strength, often chosen to be a random distribution. T h e bond represents the system on a mesoscopic scale and is characterized by a constitutive law (e.g. force F linearly dependent on the displacement 5). T h e relation between F and 5 can be chosen to remain linear up to a breaking threshold 5 . T h e C strength of the bond is given by the value of F at 6 = 5 and the elastic modulus by C the slope of the straight line. Disorder can be implemented by randomly choosing the strengths a n d / o r the elastic moduli of the bonds of the lattice. In quenched systems, the random variables (strength and elastic moduli) are fixed at the beginning of the simulation and the following rupture process is deterministic: to break a bond, the equations are solved and the bond for which S/S is the largest is chosen to be broken. c O n the other hand, annealed disorder leads to a probabilistic rupture process. In this case, the bond to be broken may be chosen according to a probability p oc d~ , where n is n a given exponent [66]. Finally, regarding step (4), once the bond to be broken has been chosen, its actual removal or replacement has to be implemented. In the case of fracture, the elastic Chapter 2. Lattice Model for the Study of Earthquakes 15 modulus of that bond can be set to zero or reduced by a certain amount. 2.3 A N e w Lattice M o d e l T h e central-force, bond-bending and beam models (vectorial models) are physical in the sense that the basic elements (Hookean springs or elastic beams) have a real physical meaning. However, these models have some drawbacks [45]. One of them, and surely the most severe one, is the computer time involved. It is enormous due to the fact that, at each time step, one has to solve a huge array of linear equations. T h i s is a major obstacle to progress in this area. T h e (tensorial) lattice model that we use was introduced by X u et al. [113, 114] and it overcomes some of the difficulties encountered in the above vectorial lattice models. It proposes a new method of discretizing the continuum elastic equations and has some advantages over the latter. 2.3.1 A new discretization scheme X u et al. have introduced recently a new scheme to discretize the continuum elastic equations. Because it is not as well known as the discretization schemes presented in [45], we summarize here the most important steps of this new scheme. T h e interested reader is referred to [114] for the details. T h e first step is to construct discrete derivatives. We do this only for a two-dimensional system for simplicity. (Generalization to a three-dimensional system is straightforward in principle). These derivatives are of the form (say for a function 9\ r _ x + y\ + —T" (_ +9 [r +— — x — x — y\ x + y\ + -^r- (_ ~ 9 \r H \-9\r _ x + y\ — g(f)) (_ x — y\ (_ — ) + g[r x — -g[r y — and 9\ r y\ — (_ - g\r x + y —- Chapter 2. Lattice Model for the Study of Earthquakes 16 where x and y are unit vectors oriented along the X-axis and Y-axis respectively, and the lattice constant (a) is taken to be unity. Note that if g sits on the nodes of a square lattice, then D g (a = x, y) sits on the centers of the squares formed by the nodes of the a lattice, and vice versa. T h e above derivatives ensure that the basic unit of discretization is a square instead of a bond. Using these discrete derivatives, the discretization of elastic theory is straightforward. First, a displacement vector u is defined on each node (corner of a square). T h e distortion of a square is characterized by the strain tensor (defined at the center of each square): Initially, all deformations are elastic and the stress tensor is related to the strain tensor through the special Hooke's law: o-ij{r) = X'SijUuif) + 2p,Uij(f), where / J is the shear modulus, and A' a Lame constant modified to take into account the plane-stress geometry of the problem (see [59] for instance). indices Summation over repeated = x,y) is implied. T h e discretized elastic equations are given by the force balance condition = DjOij = 0 which must be satisfied on each square. We further restrict our consideration to the shear mode of rupture (most ruptures in earthquakes are shear mode fractures!). In this case, only two stress components are needed to determine where a rupture will occur: o\ = o xy ° 2 = {o-yy — a )/2 xx (rupture along the X-axis or Y-axis) and (rupture at 4 5 ° or 1 3 5 ° to the X-axis). In this thesis, we consider only the stress component o\ = a xy for simplicity. We do not think that to have two stress components would have a major influence on the results presented in this thesis as well as in [42]. Now, a rupture is assumed to occur at the square centered at f*o in an infinite medium. T h i s happens when a (f ) xy 0 > a , where a h is the stress threshold of that square. T h e th t Chapter 2. Lattice Model for the Study of Earthquakes 17 additional shear stresses (cr|-) caused by the rupture are separated into an elastic part (ofj) and a non-elastic part (at r 0 only) [to express the violation of Hooke's law at this location]. T h e new stress tensor components (cr™™) can be expressed as a function of the stress tensor components before the rupture (cr°j ) in the following way: d ^(f)=ag (f) + d ay(f), where = - (2-1) f ^ r o <4(*) = °tif) (2-2) °' y(f) = °S WV At equilibrium, o™ w (-) 2 f must satisfy the force balance condition on each square and as a result, DjO-'ij = 0 since Djoff (2.4) = 0. Using (2.1), (2.2) and (2.3) in (2.4), it can be shown [114] that L>g + Ff = 0, where and -T„, = ;= 0 \Pi T , T , x+v — 2 0 r i+y > 0--p r — 0 i-v r,ro-^ -\- 0 , x—v r,r +-2 J I 0 Therefore, afj can be viewed as generated by the external source F d = (F ,F ), d d which is only non-zero at the corners of the ruptured square, and is represented on figure 2.1. F d satisfies the conditions that its net force and net torque are zero, and thus is a double couple (a term used in seismology to describe the earthquake source having double couple 3 Chapter 2. Lattice Model for the Study of Earthquakes - 18 f Figure 2.1: One-layer planar region considered in the model of X u et al. The region is discretized into squares. It is embedded in an infinite medium (providing the necessary external drive) and the effect of any activity originating outside the region is neglected (open boundary conditions). In the model, the interaction between the external drive (represented by the arrows outside the region) and a square of the region is treated in a mean-field way, i.e. it is the same whatever the square location. The arrows at the corners of one of the square inside the region are a schematic representation of a double couple force redistribution following the shear rupture of that square. This figure was modified from [114]. 19 Chapter 2. Lattice Model for the Study of Earthquakes force / ) . The discretization scheme of X u et al. is in agreement with the fact that a shear rupture can be modeled by a double couple force redistribution [11]. Next, ofjif) can be determined by using the above earthquake source and the method of Fourier transformation. Since o is the only relevant component of the stress tensor, xy we drop the subscript (from now on, stress simply refers to the shear stress component a ). It can be shown [114] that the redistributed field is given by xy o-'(r-)=a' (r) = -fG(r-f ), xy (2.5) 0 where / = fV2(X' + p)/{X' + 2p) and G is a lattice Green function expressed as dk f dky sin k sin k — r d 2TV J-TT 2TT 2-K (1 — cos k kx cos kky) y)' w G(f) = f V J-ir d x k 2 S i n 2 k y x S xi f l 2 k 2 y t y 2 f (2 6) x In appendix A, we illustrate how to simplify the calculation of G by reducing the number of integrals in (2.6) from two to one. This has a positive impact on the accuracy of the computation of G(f— f*o) especially at large value of f— r . Also, in appendix B , we do a 0 coarse-graining procedure to overcome an unphysical property of G. Our coarse-graining procedure is different from the one introduced by Chen et al. [16] in a study of the model of X u et al. in the sense that it is not implemented in the derivation of the equations of the model, but instead applied directly to the Green function (2.6) [see appendix B]. In the model, the unit of discretization (square) may be sufficiently large that the whole square will not break in an individual rupture event. Instead, it will retain some shear stress. Assuming that a (r )=Xa (r ), new (2.7) old 0 0 where X is the fraction of the original stress which is not released, (2.5) can be rewritten as a (f) = - ( 1 - X)a (f )G(rold 0 f )/G(0). 0 (2.8) Chapter 2. Lattice Model for the Study of Earthquakes 20 After the original rupture at r , which is followed by stress redistribution, it can 0 happen that the stress exceeds the stress threshold on several squares. We will then let them all break independently and add the respective stress redistribution contributions. This procedure was used as well in [113, 114]. Other rules for the order in which the squares break could have also been considered. For instance, we could have picked only the square whose stress exceeds its stress threshold most. We could, as well, have considered a self-consistent scheme in which the redistributed stresses from different broken squares are coupled as in [6]. We employ the independent rupture procedure because it is the simplest. 2.3.2 Advantages of that discretization scheme over others The discretization scheme introduced by X u et al. is rather simple and easy to implement. Before applying it to the study of earthquakes, we mention the advantages it has over discretization schemes discussed in [45]: (i) Each square is characterized by a strength. The length scale is set by the size of the square. This makes it possible to study engineering type problems (e.g. crack growth in ceramics) as well as geological scale problems (e.g. plate tectonics and earthquakes). (ii) The stress tensor is dealt with directly, so the simulation process is explicit in its physical interpellation. (iii) When open boundary conditions are used, one focuses on a finite region embedded in an infinite medium and the effect of any activity originating outside the region is neglected. In this case, the method of X u et al. becomes straightforward since the calculation of the lattice Green function can be done as for an infinite region. Chapter 2. Lattice Model for the Study of Earthquakes 2.4 21 A p p l i c a t i o n to the Study of Earthquakes The phenomenology of earthquake occurrence is complex (for a review, see [55]). Many statistical relations have been deduced over the last century from the observation of earthquakes, which can be viewed as sequences of ruptures releasing the elastic stress built up slowly by the movement of tectonic plates. Among them are the well-known Gutenberg-Richter [43] and Omiri [78] power laws. The former relation expresses the frequency of occurrence of earthquakes of a given magnitude whereas the latter expresses the aftershock occurrence decay with time from the mainshock. The global dynamics of earthquakes is rich and complex. There have been many attempts at earthquake modeling, yet there are still many unsolved problems due to the fact that one does not know how to model the physics well [55]. Knopoff suggested to introduce at least five constraints when modeling the complex phenomenology of earthquake occurrence [55]: (i) plate tectonics to restore the energy dissipated in earthquakes, (ii) a wide range of earthquake sizes, (iii) a brittle-ductile rheology to provide for the time delays between clustering events (aftershocks in particular) and the abrupt occurrence of ruptures, (iv) stress redistribution upon rupture to provide the mechanical coupling between earthquake events, (v) a complex geometry, to take into account the distribution of earthquakes as well as the finite size of the fracture surfaces. The model of X u et al. is in principle able to incorporate those five constraints as well as the tensor character of the stresses (the latter is also believed to be important for Chapter 2. Lattice Model for the Study of Earthquakes 22 earthquake modeling [55]). However, in practice, aftershocks will appear in this model only if a static fatigue law is introduced by hand into the model [113, 114] (static fatigue is believed to be the mechanism responsible for aftershocks [89]). Also, the model is twodimensional and only a one-layer system is considered (see figure 2.1), so the complex geometry of fracture surfaces is approximately dealt with. In the version of the model of X u et al. used in this work and in [42], aftershocks are absent since no static fatigue law is considered. However, even with this simplification, we will see that this version of the model still has many complex features. To apply the model of X u et al. to the study of earthquakes, we implement the following algorithm, divided into twelve steps: 1. Consider a finite (planar) region embedded in an infinite medium, neglecting the effect of any activity originating outside the region (open boundary conditions). 2. Divide the region into L\ x L squares each having a size equal to unity. 2 3. Set the stresses on all squares to zero at the beginning of the simulation. 4. At the beginning of the simulation, introduce disorder by assigning random numbers between Q[ and 9 to the stress thresholds (cr ) of the squares. The random U t/l distribution of stress thresholds is taken as uniform in this work. 5. Drive the system, i.e. increase the stresses uniformly until on a square the stress exceeds the stress threshold (this is called a rupture). The amount of stress added to every square is ACT. A n earthquake sequence then begins and the driving process stops. 6. Increment the long time scale (defined by variable t throughout this work) by the amount ACT. 23 Chapter 2. Lattice Model for the Study of Earthquakes 7. Reset the stress threshold of the ruptured square to a random number between 9i and 9 (rule I) or to the constant value 1.0 (rule II). U 8. Calculate the stress redistribution [using (2.8)]. 9. Check if other ruptures occur in the region afterwards. If so, then step 7 is repeated for all the ruptured squares. 10. Repeat step 8 for all the ruptured squares (independently) and do step 9 again. The earthquake sequence stops when the stress is lower than the stress threshold on all squares. 11. Count the number of units (S) that have ruptured during the sequence. a measure of the magnitude of the earthquake. S is Note that if a square ruptures more than once during the sequence, every occurrence is counted (unless otherwise specified). 12. Go to step 5 to generate the next earthquake sequence. We describe next most of the above steps in some details. Step 1 describes the region and the boundary conditions considered. Physically, the planar (two-dimensional) region represents a thin piece of the earth's crust, or at least a region for which the largest horizontal dimension is much larger than its thickness (vertical dimension) [see figure 2.1]. Therefore, it means that as soon as earthquakes occur in this one-layer region, the faults generated break the region over its entire thickness. The activity happening outside this region is neglected. However, the medium surrounding the region of interest (see figure 2.1) is essential since it provides the necessary external drive (see step 5 for the details). The boundary conditions that represent best such a situation are the open boundary conditions. The latter greatly facilitate the calculation of the lattice Green function (2.6) since it can be done as for an infinite region. Chapter 2. Lattice Model for the Study of Earthquakes 24 In most of the models studying S O C , the continuous space is discretized into elements (step 2). It is the existence of such distinct elements that is the basis of the complex response of the model to an external drive. T h i s is due to the fact that each element acts as an oscillator of relaxation. T h e coupling between many such oscillators with different characteristics is known to generate complicated behaviors. It could be argued that the results presented in this work are specific to the discrete nature of the model of X u et al. and are thus not relevant to real geology. However, as suggested by Sornette et al. [95], it is clear that the presence of heterogeneities in the earth's crust generates a distribution of characteristic scales. Somehow, models such as the one of X u et al. provide simple ways to represent the heterogeneities which are rescaled from the infinitesimal scale up to the smallest scale considered in the models, i.e. the lattice constant (a). In step 4, we take into account the fact that the region modeled is heterogeneous. However, it is of interest to note that disorder is incorporated solely in the stress thresholds and not in the elastic constants. T h e latter form of disorder would destroy the translational invariance of the lattice Green function as mentioned in [114]. T h i s represents one limitation of the model. T h e driving of the region (step 5) is done slowly and uniformly. T h e slow driving (this defines the long time scale) is not an artificial constraint imposed but, as we have seen in chapter 1, is fundamental to get S O C in a model. T h i s constraint makes sense for plate tectonic problems. T h e uniform driving condition is a mean-field of what happens (or happened) in nature. approximation Indeed, in the case of the Asian continent deformed by the indentation of India for instance, the stress increase was distributed non-uniformly over large distances. In such a case, the model of X u et al. would treat the interaction between the external drive (India in the above example) and the region (the Asian continent in the above example) as being independent of the position in the region. T h i s is another limitation of the model. Chapter 2. Lattice Model for the Study of Earthquakes 25 In step 6, it is implied that the rate of stress increase (p) is set to the value 1. T h e choice of such a value is motivated by practical considerations: because no aftershocks or foreshocks can happen in the current version of the model of X u et al., it is not necessary to use a small value of p (e.g. p = 10~ ) 9 since no rupture can occur between two earthquake sequences. A s a result, the amount ACT is added as a whole rather than being added little by little to all the squares of the lattice, and this makes the computation much faster. It is of interest to note that the value of p determines the order of magnitude of the long time scale (defined by variable t): t ~ p~ l begun, t is frozen. T h i s reflects the fact that t/t seq = 1. Also, once a sequence has 3> 1, where scale (or time scale for the duration of an earthquake sequence). t seq is the short time A s a result, the two time scales in the model are well separated, a fact which is important to obtain S O C in a model as discussed in chapter 1. T h e resetting rules of the stress threshold of a ruptured square are defined in step 7. Rules I and II do not have a physical basis contrary to the rules used by X u et al. in [113, 114], for which slow annealing of the fracture surfaces was allowed. In particular, rules I and II avoid the problem of having to deal with the whole range of time scales implicit in the annealing process of a fracture surface. T h e big advantage of using rules I and II is that they allow for the easy investigation of the effect of disorder on the scaling properties of the model (see chapter 3). Step 10 defines an independent redistribution procedure when more than one square break simultaneously. W h e n such a procedure is used, the model of X u et al. is said to be a series model (according to the terminology in [35]). T h i s means that the time for the stress on a ruptured square to be released is shorter than the time needed to redistribute (transmit) that stress to the other squares. Therefore, we have a series of stress releases, followed by a series of stress transmissions, followed by a series of stress releases and so on. We will see in chapter 4 that this independent redistribution procedure gives rise to Chapter 2. Lattice Model for the Study of Earthquakes 26 earthquake sequences that are temporally not correlated. T h e model of X u et al. is also said to be quasistatic because it does not describe the details of the dynamical process of rupture, but only the stress distribution before and after an earthquake sequence (treated as an instantaneous event on the long time scale). In other words, the inertial forces on the long term, which are very small compared to other forces, are neglected and thus, the region is always in mechanical equilibrium. Therefore, the model of X u et al. describes the avalanche analogue of a fracture, or a moving dislocation [97]. According to Sornette et al. [95], a series (or dislocation) model is more adapted to the study of large earthquakes, i.e. the ones which break the entire thickness of the earth's crust. T h e latter situation is precisely the one which we consider in this work (see the above paragraph related to step 1). In this work, we estimate the magnitude of an earthquake in a simplified way by counting the number of squares that have broken during an earthquake sequence (step 11). T h e use of such an estimator of the earthquake magnitude (i.e. S) is convenient since most of the work that has been done on earthquake models considered a similar estimator. A s a result, the comparison of our results to the ones obtained by other people will be easier. Obviously, other estimators of the earthquake magnitude could have been used, such as a^j — o^aft, where &b f and a / e a t are the average stress in the system before the earthquake and the average stress in the system after the earthquake respectively. However, we do not believe that the scaling properties of the model of X u et al. chapter 3) will be qualitatively (see dependent on the estimator of the earthquake magnitude used, and so, we consider only the above definition in the present work. T h e algorithm outlined above describes how to apply the model of X u et al. to the study of earthquakes. In its current version (and also in [42]), this model conserves the stress if the small fraction of it lost at the boundaries of the region due to the long-range stress redistribution is ignored. (This fraction goes to zero as the size of the region goes Chapter 2. Lattice Model for the Study of Earthquakes 27 to infinity). In addition, stress dissipation can be incorporated into the bulk of the region as well by multiplying a (f) at all rexcept r*o (rupture location) by a factor a between 0.0 and 1.0 [see (2.8)]. T h e reason for introducing parameter a is to investigate the effect of the degree of stress nonconservation on the scaling properties of the model, a task which will be performed in chapter 3. In all cases, the system organizes itself into a stationary state where on average the stress added to the region balances the stress removed from it. T h e stationary state of the model can be obtained by initially discarding a large number of earthquake sequences. Afterwards, by generating a sufficient number of subsequent sequences, it can be checked whether the cumulative size-frequency distribution obeys the Gutenberg-Richter power law [43] C ( S ) oc S~ , T where C(S) (2.9) is the fraction of sequences during which S or more squares have ruptured, and r is the exponent of the power law. 2.5 C o m p a r i s o n of the M o d e l of X u et al. to O t h e r M o d e l s Before we study the scaling properties of the model of X u et al. (see chapter 3), we compare the latter to other earthquake models and to models of integrate-and-fire oscillators. 2.5.1 M o d e l of X u et al. versus other earthquake models Friction and fracture are two equally important phenomena occurring during earthquakes. T o our knowledge, only the model of M o r a and Place [71] is able to incorporate both friction and fracture in a unique model. Most of the earthquake models deal with either friction or fracture. T h e models dealing with friction only are stick-slip ones dealing with fracture only are lattice models and the models (see section 2.2). T h e model of X u et Chapter 2. Lattice Model for the Study of Earthquakes al. 28 (presented in section 2.3) is a lattice model. We now describe in some details the model of M o r a and Place, stick-slip models and other lattice models applied to the study of earthquakes. T h e model of M o r a and Place is based on molecular dynamics. While it is infeasible to directly simulate the number of particles required for geophysical scales, this model can be used to simulate larger units such as pieces of rock. T h e interactions between the particles are specified through potentials which are compatible with geophysical observations. U p to now, the number of particles considered in the model of M o r a and Place has not been sufficient to simulate in a realistic way the behavior of a two-dimensional tectonic plate (see [25]). Stick-slip models focus on regions with a pre-existing fault. One of the first model of this type was the one-dimensional model of Burridge and Knopoff [12]. It consists in a chain of blocks lying on a rough stationary surface and connected by springs. Also, the blocks are all connected by flat springs to a top surface moving with a certain speed relatively to the bottom surface (system driving). T h e dynamics of the model is defined as follows: as long as all the Fj's (Fi is the total force on block i) are smaller than a constant threshold value (F ), th then all the blocks are stuck. A s the system is being driven, eventually one block i will have Fi > F h and the excess force will be relaxed. t T h e excess force will be transferred equally to the two neighboring blocks, which can in turn slip if they overcome the static friction force with the bottom surface. A chain reaction (sequence) might thus happen, which can involve an arbitrary number of blocks (N). T h e frequency distribution of sequences involving N slipping blocks was shown not to be a power law by Burridge and Knopoff [12]. T h e model of Burridge and Knopoff has inspired a great deal of work in the last decade. For instance, Carlson and Langer [13] added a nonlinear friction law to the original model and got the power law that Burridge and Knopoff were unable to obtain. Chapter 2. Lattice Model for the Study of Earthquakes 29 Nakanishi [73] also got a power law by studying a coupled map lattice version of the original model of Burridge and Knopoff. Despite these findings, the model of Burridge and Knopoff as well as its many versions were shown to N O T be self-organized critical [93]. O n the other hand, the model of O l a m i et al. [77] was shown to display S O C [17]. T h e latter model was also inspired from the model of Burridge and Knopoff for its derivation but, apart from this, it can be regarded as a two-dimensional generic representation of a nonconservative model, since the final equations of the model of O l a m i et al. do not contain the physical constants of the original model. In this sense, we do not consider the model of O l a m i et al. as a stick-slip model. Stick-slip is a label which will be reserved for the models displaying explicitly the physical constants of the model of Burridge and Knopoff. T h e stick-slip models are in principle able to solve the fully dynamic problem where, during an earthquake, a rupture propagates with a speed close to that of shear waves. However, these models are too simplified from a geometrical point of view: the region modeled is that of a pre-existing fault with many asperities coupled only to their nearest neighbors and so, nothing is said about how the fault was generated and about the effect of the surrounding medium. T h e creation of faults is related to the fracture phenomenon and can be modeled using lattice models such as the ones discussed in sections 2.2 and 2.3. These models use a continuum (macroscopic) description of matter as opposed to the discontinuous (microscopic) description used in the model of M o r a and Place [71]. T h e continuum description of matter allows for the easy prediction of where the rupture will occur but poorly describes the evolution of discontinuities such as faults [27, 106]. A p a r t from the models discussed in [45], other models using a continuum description of matter include the ones of Chen et al. [15] (see also a modified version of it in [6]) and X u et al. [113, 114] (see section 2.3). T h e latter have been applied to the study of earthquakes. Chapter 2. Lattice Model for the Study of Earthquakes 30 T h e y differ in the discretization scheme used to solve the continuum elastic equations. In [15] (and [6]), the basic unit of discretization is a bond whereas in [113, 114], it is a square. A s a result, the earthquake source is a dipole in [15] (and [6]) and a quadrupole (double couple) in [113, 114]. A rupture is followed by long-range stress redistribution in both models however. T h i s should be contrasted with stick-slip models and the model of O l a m i et al. [77] in which the stress redistribution is to the nearest neighbors only. 2.5.2 M o d e l of X u et a l . versus models of integrate-and-fire oscillators A class of models, which are mainly applied to the study of biological systems, has appeared in the last decade. These models focus on large assemblies of oscillators which can spontaneously evolve to a state of large scale organization. Collective synchronization is the best known phenomenon of this kind. T h i s effect has attracted much interest in the study of fireflies [10], pacemaker cells of the heart [82], cells of the pancreas [90], etc. Most of the works on synchronization have used models in which the interactions between the oscillators are smooth and continuous in time. Comparatively few works have been done with models in which the interactions are episodic and pulselike, although these models are relevant to the study of certain species of fireflies for instance. Models of this type are often called models of Integrate-and-Fire Oscillators (IFO). Models of I F O are closely related to the models of X u et al. [113, 114] and Olami et al. [77]. T h i s can be seen by identifying the state variable of earthquake models (usually a shear component of the stress tensor) with the voltage magnitude of models of I F O . It is of interest to note that the models of I F O study the level of mutual entrainment (i.e. the synchronization) between the elements whereas the earthquake models focus on the size distribution of the sequences generated (if it is a power law, then the model might display S O C ) . Considering the fact that both of these classes of models contain the same basic ingredients then, from this point of view, S O C and synchronization can Chapter 2. Lattice Model for the Study of Earthquakes 31 be considered as two sides of the same coin. For completeness, we mention that, contrary to earthquake models, the models of I F O could allow for more general driving processes than slow uniform driving (e.g. random driving of individual elements). Models of I F O exist in at least two varieties. In some models [69, 7], the oscillators are equally coupled with each other, so the firing of one oscillator leads to an equal redistribution of voltage to all the other oscillators in the system. In other models (see [19] for example), the oscillators are coupled only to their nearest neighbors, so the firing of one oscillator leads to a voltage redistribution to the nearest neighbors only. Therefore, as in earthquake models, the firing of one element leads to the redistribution of the state variable either to all the other elements or to the nearest neighbors of that element only. Chapter 3 Study of the Scaling Properties of the Model The scaling properties of the model of X u et al. [113, 114] are studied by means of a finite-size-scaling analysis. In section 3.1, we give the motivation for such a study and in section 3.2, we present how to implement finite-size scaling in the context of SOC in an earthquake model. The results of the analysis, obtained by considering one-dimensional lattices, are presented in section 3.3 in both the conservative and non-conservative cases. The stress ordering in one-dimensional lattices is investigated in section 3.4. This investigation is useful in explaining some of the results obtained in section 3.3. We introduce in section 3.5 a mean-field version of the model. Finally, in section 3.6, we summarize the important results obtained in this chapter and draw our conclusions. 3.1 M o t i v a t i o n for this Study The Gutenberg-Richter power law (2.9) is obtained only in the limit of an infinite system. Obviously, the power law may also appear for a finite system, but it will prevail only over a limited range. The system is never completely scale invariant. It is then important to sort out effects due to the finite size of the system and those due to the finite size of the basic unit of discretization. Such a problem arises because systems of interest are heterogeneous on scales much smaller than the discretization unit, so that the continuum limit is not meaningful. Quantities relating to the stress threshold (a h), the degree of stress retained t by a ruptured square (X) [see section 2.3] and the degree of stress conservation (a) [see section 2.4] may therefore depend upon the size of the discretization unit and the size 32 Chapter 3. Study of the Scaling Properties of the system in a non-trivial way. of the 33 Model T h i s type of dependence was realized by Weibull [110], but has since often been forgotten. T o illustrate more specifically these ideas, we discuss the effect of parameter X introduced in section 2.3 [see (2.7)]. T h e case X = 0 corresponds to a rupture opening the whole element, as in [6]. T h e interpretation of a rupture event and the scaling properties are somewhat different if X parameter X ^ 0. We expect (as well as the stress threshold distribution) to depend on the size of the discretization unit. In fact, two discretizations of the same system can be considered, one with a coarser mesh than the other. T h e n , a coarser mesh element will rupture more often than a finer mesh element, but the total number of rupture events in the whole system will be the same. These ideas will be illustrated quantitatively in section 3.3. Disentangling the effects due to the finite size of the system from those due to the finite size of the basic unit of discretization is the motivation of the study in this chapter. To achieve this, we need to vary the various parameters of the model of X u et al. F r o m what we said above, it is clear that varying the parameters of the model is not merely a mathematical procedure: it is a physically motivated one. T h e results obtained from the study in this chapter will ultimately serve to confirm whether the model of X u et al. really displays S O C . Before we present our results, we briefly describe finite-size scaling. 3.2 Finite-Size Scaling It is now well known that a finite system can not have a true singularity in its thermodynamic properties at a non-zero temperature. For instance, this is easily seen in the case of the Ising model in any dimension. Experimental systems, even if very large, are always finite. For all practical purposes, true phase transitions do occur in these finite systems and this fact must be reconciled with the rigorous proof that no phase transition can occur in any finite system. T h i s is one of the roles of finite-size scaling. A s an illustration Chapter 3. Study of the Scaling Properties of the 34 Model of this, we consider a L x L x L system. Such a system falls into the category mentioned above of a system without a true thermodynamic equilibrium phase transition. In this situation, the susceptibility and specific heat remain finite at all temperatures and the power law divergences of the infinite system are replaced by rounded peaks. Finite-size scaling describes the dependence of the peak height and the range in temperature over which the power-law behavior is observed on the system size L. Finite-size scaling has been extensively used in the study of equilibrium phase transitions (see for instance [84]). Its use in the present work is justified by the fact that, as argued by Sornette et al. [94], S O C can often be mapped onto ordinary criticality. The way to implement mathematically finite-size scaling, in the context of S O C in an earthquake model, is to replace the Gutenberg-Richter power law (2.9) (valid in the limit of an infinite system) by C(S,L) = S- [{S-1)L- ], T where L is the largest dimension in a L\ x L exponent that expresses how the (3.1) V 9 finite-size 2 lattice, g a scaling function and v an effects scale with the size of the system. W h e n comparing systems with different parameter values, we also attempt C(S, L) = S- g[(S - \)L- /f(X, T v 9 a)], h (3.2) where we assume that parameter 9 , the upper bound for the stress threshold distribution U (see step 4 of section 2.4), is set to a fixed value. T h e function / is included to test our hypothesis that, if the mesh size for a given system is changed, then some of the parameter values must also be changed. T h e term (S — 1) in the argument of the scaling function ensures that g(0) = 1, i.e. when scaling, we map the size interval [l,oo) on [0, oo). It is of interest to note that the function g acts as a correction factor for the finite size of the system and the finite size of the discretization unit. Chapter 3. Study of the Scaling Properties We are now ready to do our of the finite-size-scaling 35 Model analysis. O u r hope is that finite-size scaling will turn out to be as useful for S O C as for ordinary criticality. 3.3 Results of the Finite-Size-Scaling Analysis We consider lattices with L x L units, where L 2 2 = 1 (unless specified otherwise). T h e lattice configuration considered is therefore a strip (i.e. a one-dimensional lattice). We believe that the scaling properties of the model of X u et al. will be sensitive to this particular lattice configuration though, in the sense that with two-dimensional lattices, we would probably find different values for the exponents r and v [see (3.1) and (3.2)}. However, the use of one-dimensional lattices allows one to go to large lattice sizes, and check if the scaling laws are valid over the whole range of lattice sizes. Now, regarding the parameters of the model, 9 is fixed to the value 1 while 9i, X U (see subsection 3.3.1) and a (see subsection 3.3.2) are allowed to vary. In particular, parameter a can be set to 1 (conservative case) and to values smaller than 1 (nonconservative case). Also, the degree of disorder in the distribution of stress thresholds (ath £ [9i,9 = 1.0]) can be changed by using one or the other resetting rules (see step 7 in u section 2.4). T h e use of rule I maintains the degree of disorder of the initial distribution of stress thresholds, whereas the use of rule II makes disorder disappear after a small number of iterations. Below, the latter two situations will be respectively identified as random a h and constant a ht t We let the model reach a stationary state by discarding the first 10 earthquake 5 sequences and include 10 sequences to do the finite-size-scaling analysis. 6 Chapter 3.3.1 3. Study of the Scaling Properties of the Model 36 Results in the conservative case (a = 1) We show in figure 3.1 the cumulative size-frequency distribution C(S,L) and various L values. for random a th T h e main features are the power-law behavior for small events, and the L-dependent fall-off of the frequency of large events. T h e power-law behavior extends over two or three decades depending on L. In figure 3.2, we note that the data sets of figure 3.1 obey quite well the scaling behavior predicted by (3.1) [see the middle curves]. T h e exponent r = 0.30 was found and cannot be changed by more than 0.01 before there is a noticeable deterioration in the data collapse. A s can be seen from figure 3.1, this value of r is smaller than the one obtained from a linear fit to the straight segment of the log-log plot (especially for L = 200 and L = 400). However, the exponent r is the same if we count the number of ruptures on distinct units or if we count all ruptures. T h e exponent v, on the other hand, is different in these two cases, as is the scaling function g. W h e n all ruptures are counted, the same exponents are found for the L x 1 and L x 2 lattices, but the scaling functions are different. We attempted a data collapse of the data sets for the L x 1 and L x 2 lattices (when all ruptures are counted) by rescaling the horizontal axis for one of the sets, but found it to be significantly worse than that of figure 3.2. We now turn to the effect of having constant a h- From figure 3.3, we see that large t events become less common. We were unable to find a scaling function which fits both small and large events for L x 1 lattices, no matter if all ruptures or ruptures on distinct units (results not shown) are counted. However, from figure 3.3, it is clear that the shoulder broadens as L increases, so C(S, L) might scale with L if the small events (i.e. those with S < 50) are not considered in the data collapse. We found that, for r = 0.30 and v = 1.38 (values obtained above with L x l systems having random ath), the scaling function (3.1) works well with systems having constant a h when only the events with t Chapter 3. Figure 3.1: Study of the Scaling of the Model 37 Log-log plot of the cumulative size-frequency distributions for L x 1 lat- tices, random a h, X t C(S) Properties = 5-°- . 3 0 = 0, 9i = 0.25 and a = 1. For comparison, we show the curve Chapter 3. Figure 3.2: L x Study of the Scaling Properties of the Model 38 Log-linear scaling plots of the cumulative size-frequency distributions for 1 and L x 2 lattices, random a h, X = 0, 9i = 0.25 and a = 1. A l l curves use the t fitting parameter r = 0.30. For the left and middle curves, all ruptures are counted while for the right curves, only ruptures on distinct squares are counted. T h e left curves were obtained for L x 2 lattices using v = 1.38, the middle curves for L x 1 lattices using v = 1.38, and the right curves for L x 1 lattices using v = 1.01. Chapter Study 3. of the Scaling Properties of the Model 39 S > 50 are included in the fit (results not shown). T h i s shows in particular that exponent r = 0.30 is universal in the sense that it is independent of the degree of disorder in the model. O n the other hand, the crossover region (S < 50) on figure 3.3 is not a universal feature since it appears only for systems with constant o ht We have checked that the crossover region is a robust feature which is not sensitive to the initial configuration of the model. Therefore, even if an ensemble average of initial configurations is carried out (see [45]), the crossover region remains. T h e presence of the crossover region will be explained in the next section by studying the stress ordering in one-strip systems with constant o ht We next investigate the effect of the parameter X controlling the fraction of stress retained by a ruptured square. T h e data collapses on figure 3.4 are seen to be good for systems with random o~ h- Also, for a fixed L, we note that an increase in X extends t the power law to larger S. T h i s makes sense because a larger X means that a ruptured element retains more stress, so it will be easier to have it involved in the near future in an earthquake sequence (because stress remains close to stress threshold). A s a result, the largest size for an earthquake sequence increases. In figure 3.5, we show, for a 1000 x 1 lattice and random a h, a fit to the scaling form t C(S,L) = S- g[(S-l)L-»/f(X)}. T T h e fit, although not perfect, implies that a system which has been discretized by a fine mesh and a small value of X (i.e. a ruptured square almost opens) can be approximated by a coarser mesh and a larger value of X (i.e a larger fraction of the stress is retained by the ruptured square). T h e effect of changing the lower bound 0; for the uniform stress threshold distribution is shown on figure 3.6. T h e data collapses are very good for systems with random a ht For a fixed L, we note that a decrease in Q\ extends the power law to larger S. Comparing Chapter 3. Study of the Scaling Properties of the Model 40 Figure 3.3: Log-log plot of the cumulative size-frequency distributions for L x 1 lattices, constant a h, X = 0, 0/ = 0.25 and a = 1. The segment S < 50 is a crossover region and its origin is explained in section 3.4. Beyond 5 > 50, C(S,L) has a linear segment with a slope close to 0.30. This linear segment gets longer as L increases but cuts off at large S. For comparison, we show a curve proportional to 5 . t - 0 3 0 Chapter 3. Figure 3.4: Study of the Scaling Properties of the Model Log-linear scaling plots of the cumulative size-frequency distributions for L x 1 lattices, random cr , 9i — 0.25, a = 1 and two values of X. th T = 41 0.30 and u = 1.38. Exponents used are Chapter 3. Study of the Scaling O.O Properties of the Model 0.5 ( S - 1 ) L _ V 42 l.O / f ( X ) Figure 3.5: Log-linear scaling plot of the cumulative size-frequency distributions for a 1000 x 1 lattice, random a h, 9i = 0.25, a = 1 and three values of X. t are r = 0.30 and v = 1.38. Exponents used Chapter 3. Study of the Scaling Properties of the 43 Model figures 3.4 and 3.6, we see that, for a given L, an increase in X has much the same effect on the power law as a decrease in We have also attempted a fit to the scaling form C{S,L) = S-rgftS -l)L-»/M)]. T h e quality of the fit (see figure 3.7) is comparable to the one on figure 3.5. T h i s result implies that a system which has been discretized by a fine mesh and a large value of 9i, can be approximated by a coarser mesh using a smaller value of B{. We expect that this near equivalence can be made even better if we simultaneously adjust both X and #/. We have also investigated the effects of changing parameters X and 9i in systems with constant a ht Again, we fail to find a satisfactory scaling form which is valid for both small and large events since a crossover region at small S is still present. A s we did above, we attempted data collapses for various X and 9 by focusing only on the events t outside the crossover region and found that with r = 0.30 and v = 1.38, the scaling function (3.1) works well. T h i s shows that r m 0.30 is a universal value. Before we consider the case of non-conservative systems (a < 1), we comment on the value obtained for exponent r . We have obtained through a finite-size-scaling analysis r = 0.30 ± 0.01, which seems to be robust. T h i s value can be compared with the one obtained in [68] (0.31 ± 0.04), which was obtained with a long-range scalar model with quenched disorder (stress thresholds are kept equal to their initial random values throughout a simulation). O u r value of r is in agreement with the latter, but is not in agreement with the value (0.4) obtained by X u et al. [114] and Chen et al. [15] for a two-dimensional system. A s we have checked for a 100 x L 2 system (results not shown), the discrepancy between our value of r and theirs is significant since as Li increases, the cumulative size-frequency distribution settles onto a curve with a linear segment having a slope close to 0.4. T h i s was observed for systems with both random a h and constant a ht t (In some way, this discrepancy shows that if we would have studied the scaling properties Chapter 3. Figure 3.6: Study of the Scaling of the Model 44 Log-linear scaling plots of the cumulative size-frequency distributions for L x 1 lattices, random cr h, X t T = Properties 0.30 and v = 1.38. = 0, a = 1 and two values of 0/. Exponents used are Chapter 3. Study of the Scaling Properties of the Model 45 Figure 3.7: Log-linear scaling plot of the cumulative size-frequency distributions for a 1000 x 1 lattice, random a n, X = 0, a = 1 and three values of f?;. Exponents used are t T = 0.30 and v = 1.38. Chapter 3. Study of the Scaling Properties of the Model 46 of isotropic lattices instead of one-strip lattices, then we would have obtained a different value of r , but apart from that, qualitatively similar results). Finally, Christensen and O l a m i [17] obtained a different value (0.22 ± 0 . 0 5 ) through a finite-size-scaling analysis of the short-range, two-dimensional model with constant o h in [77]. Therefore, the models t in [77] and [114] seem to belong to different universality classes. Some light is shed on this question in section 3.5. 3.3.2 Results in the non-conservative case (a < 1) In this subsection, we examine two effects for L x 1 systems with random a h t a n d constant G h when X is fixed to 0 and 6i to 0.25. T h e first one is the effect of parameter a for a t fixed L, and the second is the effect of L for a fixed a < 1. In figure 3.8, we present the cumulative size-frequency distributions for different values of the parameter a controlling the level of stress conservation for a 2000 x 1 system with random o~ h- A decrease in a makes the power law cut-off happen at increasingly smaller t S. Also, as a decreases, the magnitude of the slope of the straight segment of the log-log plot increases, an effect which was also noticed in [17]. for a = 0.94, random a h and different lattice sizes (L). t In figure 3.9, we show C(S,L) A s expected, as L increases, the straight segment of the log-log plot extends to larger S. However, a saturation is noticed for L > 2000. We have fitted the cumulative size-frequency distribution for fixed L and a (see the dashed lines on figure 3.8) to C(S) = S~ T exp[(l - S)/S + bS 2 0 + cS + 3 ...}. (3.3) For large L, we found that only the first term in the argument of the exponential is significant, where So is a characteristic earthquake size. In this case, (3.3) reduces to a form similar to the one fitted in [58]. In table 3.1, we give the values of r and So for various values of a and L. For a fixed a, SQ{L) increases with L and appears to saturate. It is Chapter 3. Study of the Scaling Properties of the a L T 0.94 0.94 0.94 0.94 0.94 0.97 0.97 0.97 0.99 0.99 0.99 4000 2000 1000 400 200 2000 1000 400 2000 1000 400 0.482 0.477 0.490 0.419 0.354 0.428 0.405 0.361 0.386 0.367 0.277 47 Model s 0 52.2 50.0 48.2 37.2 27.4 115 98.5 70.6 407 307 149 Table 3.1: Exponent r and characteristic earthquake size So for L x 1 lattices, random o~th, X = 0, 6i = 0.25 and three values of a. possible with our data that SQ(L) diverges, but slower than L (SQ{L) ~ IP with /? <C 1) as suggested by Janosi and Kertesz [50]. The discussion in section 3.5 suggests, however, that this is unlikely. From table 3.1, we also note that, for a fixed L , So increases with a and r increases as a decreases. We therefore corroborate quantitatively conclusions obtained qualitatively from the curves on figure 3.8. We have also studied the above effects for systems with constant a h (results not t shown). Again, as in the conservative case, a crossover region is observed. Then, we found that a straight line can be fitted, for fixed L and a, to the log-log plot of C(S, L) versus S for values of S outside the crossover region up to a cut-off value much smaller than in the conservative case. We noted that the slope of the straight segment increases and the cut-off value decreases as a decreases. Also, for a fixed cv, we found that the cut-off value saturates as L increases. Because the cut-off value (for a fixed L) is smaller for a < 1 than for a = 1, we can again associate a characteristic earthquake size So to Figure 3.8: Log-log plot of the cumulative size-frequency distributions for a 2000 x 1 lattice, random a h, X = 0, B\ = 0.25 and three values of a. T h e dashed lines are the fits t to relation (3.3). Chapter 3. Study of the Scaling Properties of the Model 49 Figure 3.9: Log-log plot of the cumulative size-frequency distributions for L x 1 lattices, random cr th: X — 0, 9 — 0.25 and a = 0.94. t Chapter 3. Study of the Scaling Properties of the Model 50 C(S,L). Finally, it was observed that simulations running for longer times are providing data which are not altering our results. T h i s is in contrast to a study of the short-range model in [77], which found that the time for total invasion of the interior of large lattices by S O C can be very large, especially for values of a close to 0.0 [67]. We believe that because the model of X u et al. is long-range, the invasion is much faster than in the model in [77]. 3.4 Stress O r d e r i n g i n O n e - D i m e n s i o n a l Lattices In the previous section, we found that for one-strip systems with constant a h, the cumut lative size-frequency distribution C(S, L) exhibits a crossover region for S < 50, which seems to be robust. We also discovered that for S > 50, up to a cut-off value scaling with L, the log-log plot of C(S,L) versus S displays a straight segment (see figure 3.3). T h e goal of this section is to explain the presence of that crossover region. T o achieve this, we study the stress ordering in one-strip systems, and in particular in one-strip systems with constant a . th T o our knowledge, a study of the stress ordering in a system has only been performed by Grassberger [38] in the case of the short-range model in [77]. We begin by describing the stress ordering study and afterwards, attempt an explanation for the presence of the crossover region for systems with constant o ht The stress ordering in a system can be either local or global [38]. T o explain the presence of the crossover region, it is sufficient to study the local stress ordering. In fact, in the case of a one-strip system with constant a h, the study of both kinds of t stress ordering is redundant since if some local stress ordering is found in a system, then automatically global stress ordering also prevails in the system. Local stress ordering is studied by considering the stress differences between nearest-neighbor elements in a system, i.e. 5{ = | a i i+ — Oi \ (i = 1,2,...,L — 1 for a L x 1 system). In comparison, Chapter 3. Study of the Scaling Properties of the 51 Model global stress ordering can be studied by focusing on the fitness, defined as / ; = o\ h where i = 1, 2 , L for a L x 1 system, o\ h — o~i, being the stress threshold of the iih element. A histogram plot of the <5j's can be obtained and used to conclude whether local stress ordering exists in a system at any given time. T h e presence of peak(s) in the histogram ultimately proves that local stress ordering prevails in the system. We begin our exploration by considering a system with constant a h- O n figure 3.10, t we plotted the Si's for a 1000 x 1 system. A series of patches of various sizes are observed at this particular time (arbitrarily defined as to, where here time refers to the long time scale) for small <5;. In particular, the two biggest patches on figure 3.10 have sizes 61 (at % ~ 800) and 46 (at i « 600), which means that the number of elements at about the same stress level is 62 and 47 in these two cases respectively. Assemblies of elements next to one another, which are at about the same stress level, will be referred to as of elements below. GROUPS A s long as one element in a group has its stress over o~th, then the other elements in the group will likely rupture in the same sequence. We have checked that the groups of elements are not static features, i.e. elements at later times. they do not involve the same A histogram of the data of figure 3.10 was obtained and the result is plotted on figure 3.11 (solid line). T h e other curve on figure 3.11 was obtained with data taken at time t\ > t . 0 T h e histograms at both times display a peak at Si ~ 0, showing that local stress ordering prevails in the one-strip system with constant a ht We have also checked that local stress ordering is maintained at any time by varying parameters X and 9i (results not shown). X T h e quantitative difference a larger value of (or a smaller value of 9{) makes, as compared to X = 0 (or 9i = 0.25), is to restrict the Si's to be smaller numbers. For instance, we have checked that for X = 0.8, the Si's cut off approximately at 0.2 (as compared to about 1.0 for X = 0 on figure 3.11). A t this point, we attempt an explanation for the presence of the crossover region observed for systems with constant a th (see figure 3.3): during an earthquake sequence, 52 Chapter 3. Study of the Scaling Properties of the Model Figure 3.10: Distribution of the stress differences 6i (i = 1,2, ...,999) between the nearest-neighbor elements of a 1000 x 1 system with constant o h, X = 0, 0i = 0.25 and a = 1. t Chapter 3. Study of the Scaling Properties of the Model 53 300- 250- 200- <0 . , 0 , 0.2 ! „ 0.4 , 0.6 , 0.8 , 1 Si Figure 3.11: Histogram plots of the stress differences <5j between the nearest-neighbor elements of a 1000 x 1 system with constant o , X = 0, 6i = 0.25 and a = 1 at two distinct times. The solid line curve was obtained with data taken at the arbitrary time to (data of figure 3.10) and the curve plotted using a dashed line and lozenges was obtained with data taken at time t\ > to. Fifty bins were used to get the histograms. th Chapter 3. Study of the Scaling Properties of the 54 Model when the ruptured elements are found mostly in a given group of elements or in a couple > of small groups of elements away from other groups in the system, then the size of the sequence (S) is limited by the number of elements in the group(s), so S tends to be small (i.e. in the crossover region). O n the contrary, if the ruptured elements are found mostly in groups of appreciable sizes close to one another, then S tends to be large. In other words, the localized (versus delocalized) nature of the activity in a sequence determines whether S is small or large. We believe that this explanation is plausible considering the fact that for a 1000 x 2 system with constant o~ h, no crossover region is observed (results t not shown) since the presence of a second strip of elements allows the derealization of the activity to take place. A l o n g the same line, we have noticed that systems with random a h never display patches (and in a more general sense, local stress ordering) such t as the ones on figure 3.10 (results not shown), so the localization of the activity never happens in these systems. T h i s explains why systems with random a h always display t simple cumulative size-frequency distributions (as on figure 3.1 for instance). We now turn to a simplified version of the original model of X u et al. by considering a mean-field version of it. 3.5 M e a n - F i e l d V e r s i o n o f t h e M o d e l o f X u et a l . We present in this section a mean-field version of the model of X u et al. in which all the elements are coupled equally with each other. A s a result, it is not necessary to use a Green function and to refer to a lattice structure anymore. T h e additional stresses (a ) 1 caused by the rupture of a given square are now all expressed as • = ^ f , (3.4) where L is the number of squares in the system and the other symbols have the same meaning as in the case of relation (2.8). Therefore, the rupture of a square is followed by Chapter 3. Study of the Scaling a stress decrease [(1 — X)o ] Properties of the on this element and a stress increase old 55 Model [(1 — X)a° /(L ld — 1)] on all the other L — 1 elements, which leads to the same total stress value before and after the rupture. Such a situation leads to a total stress growing without limit in time. To prevent this from happening, we can multiply a' in (3.4) by a constant a between 0.0 and 1.0. We use the same symbol as in subsection 3.3.2, but it is important to keep in mind that in the mean-field model, nonconservation of stress is essential for the model to describe earthquakes. In the mean-field model, a should be adjusted to values close to 1 in order to not have earthquakes releasing artificially too much stress from the system. We will say more about this below. In section 3.3, we presented results obtained with a version of the original model of X u et al. Systems with random a h and constant o~ h were both considered. In the case of the t t mean-field model, systems with constant a h lead, in the stationary state, to trivial sizet frequency distributions, since groups of elements which ruptured simultaneously in the past, will rupture simultaneously in the future as well. Indeed, for a given such system, depending on the number of elements (N ) e in each group, the resulting size-frequency distribution consists of a series of peaks, one for each distinct value of N . e Therefore, the mean-field model is N O T self-organized critical when disorder in the stress threshold distribution is absent, since no power law can be obtained for the size-frequency distribution. In systems like this, it remains the possibility of studying the synchronization in the system, but we will not do such a study here (see however [69, 7] for a discussion of mean-field models of integrate-and-fire oscillators). T h e study of systems with random a h with the mean-field model is non-trivial. t The study of systems with random a h with the mean-field model can be done by t implementing the algorithm of section 2.4, except for steps 1 and 2. O n l y resetting rule I (see step 7) of the algorithm is meaningful and so is the only one we use. T o do step 8 of the algorithm of section 2.4, i.e. to calculate the stress redistribution after the rupture Chapter 3. Study of the Scaling Properties of the 56 Model of an element, relation (3.4) is used rather than (2.8). Again, we let the model reach a stationary state by discarding the first 10 earthquake sequences generated and include 5 10 6 sequences to obtain the cumulative size-frequency distribution C(S,L) for a given number of elements (L) in the system. Now, if we would set a to a given value and plot C(S, L) for a few values of L, then we would obtain a graph similar to the one on figure 3.9, with no data collapse of the distributions for different values of L. Because we want to attempt a data collapse of the distributions for different L, we use an alternate route: for a given L, a will be set equal to a / / , the effective value of a for a L x 1 system e (when a — 1) studied with the original model of X u et al. We plotted on figure 3.12 log(l — a ff) versus log(L) for L 6 [100,1000], where cv // was obtained by averaging e e over 2000 earthquake sequences generated by the original model of X u et al. with X = 0, 0i = 0.25 and a = 1. From figure 3.12, we see that log(l — a jj) decreases about linearly e with log(L). T h e slope of the straight line was found to be —1.02. T h e fact that a ff ^ a e even if these two quantities are close together results from the long-range nature of the model of X u et al.: there is always a fraction of the redistributed stress from a ruptured element which is lost at the boundaries of the system. T h i s fraction goes to zero as the size of the system goes to infinity (i.e. L —» oo). Now, using for a in the meanfield model the value of a ff e on figure 3.12 for the appropriate L, we d i d simulations of the mean-field model for L = 200, L = 400 and L — 1000. A fit of the cumulative size-frequency distributions to the finite-size-scaling function (3.1) was attempted with the values r = 0.3 and v = 1.38, but the data collapse was found to be bad (results not shown). For this reason, we simply plotted on figure 3.13 C(S, L) versus S without finite-size scaling. It is obvious from the figure that even if C(S, L) extends to larger S as L increases, it is not possible to get a statisfying data collapse for both small and large S. A s a result, the exponent r seems to be sensitive to the details of the model, i.e. in this case the fact that the elements are equally coupled with each other and not Chapter 3. Study of the Scaling Properties of the Model 57 arranged in a lattice as in the original model of X u et al. It is also possible that the trick we have used above to fix a for a given L in the mean-field model does not allow to get optimal results, which could explain why no satisfying data collapse was obtained. It is relevant to mention that a /f e obtained from the original model of X u et al. is not unique for a given L in the mean-field model since it depends on how the squares are joined together in the original model. For instance, a ff e obtained for L = 400 is not the same if the lattice configuration is a strip (400 x 1) rather than a two-dimensional 20 x 20 arrangement. T h i s study of the mean-field model was just a first attempt at trying to investigate the effect of changing the rules of the model. M u c h work along this line remains to be done. 3.6 S u m m a r y of the Important Results and Conclusions In this chapter, we have studied the scaling properties of the original model of X u et al. as well as a mean-field version of it by means of a finite-size-scaling analysis. In the mean-field model, the elements are equally coupled with each other. In the original model, one-dimensional (L x 1) lattices were considered for convenience. We summarize the results obtained in this chapter and draw our conclusions. For the original model, the finite-size-scaling function (3.1) works well over the whole range of S for systems with random a h and for S > 50 for systems with constant a ht For systems with constant a h, t t there is a crossover region for S < 50 and it appears whatever the parameter values in the model. T h e existence of such a crossover region is related to the localized nature of the activity on a group of elements or on a couple of small groups of elements at about the same stress level (and therefore likely to rupture together as long as one of them will do so) isolated from any other group in the system. Chapter 3. Study of the Scaling Properties of the Model 58 Figure 3.12: Log-log plot of 1 — a / / values versus L for L x 1 lattices, random a^, X — 0, 6i = 0.25 and a = 1 obtained with the original model of X u et al. The slope of the straight segment was found to be -1.02. e Chapter 3. Study of the Scaling Properties of the Model 59 0.1- 0.01- 0.001- 0.0001- 1e-05- 1e-0610 100 S 1000 10000 100000 Figure 3.13: Log-log plot of the cumulative size-frequency distributions for systems of L elements equally coupled with each other. The values of the parameters are X = 0 and the stress thresholds (i.e. the a h) are taken from a uniform distribution with 9i = 0.25. From left to right, the distributions refer to systems having L = 200, L = 400 and L = 1000. The values of a differ for all three distributions (see how they were adjusted in the text). t Chapter 3. Study of the Scaling Properties of the Model 60 T h e crossover region disappears for two-dimensional systems with constant a h t a s well as for systems with random o ht The effects of parameters X (controlling the degree of stress retained by a ruptured square) and 9i (controlling the width of the uniform distribution of stress thresholds) were also investigated and we found a fair agreement with the finite-size-scaling function (3.2) over the whole range of S for systems with random a h and for S larger than about t 50 (outside the crossover region) for systems with constant a h- T h i s result implies that t a system which has been discretized by a fine mesh and a small value of X (or a large value of 6i) can be approximated by a coarser mesh using a larger value of X (or a smaller value of 0i). T h e exponents of the finite-size-scaling functions (3.1) and (3.2) were found to be T = 0 . 3 0 ± 0 . 0 1 and u « 1.38. T h i s value of r is robust and in particular, it is independent of the degree of disorder in the system (i.e. it is the same for systems with random o h t and constant a h)- However, we have checked that for systems with more than one strip, t there is a crossover of the exponent r to r % 0.4. Therefore, the results in this chapter are quantitatively specific to one-strip systems. Introduction of stress nonconservation (controlled by parameter a) results in the appearance of a characteristic earthquake size (So), which is smaller than L for large L. For a fixed L, So increases with a and the exponent r increases as a decreases. The latter observation is in agreement with the findings of Christensen and O l a m i [17]. These results were found to be independent of the degree of disorder in the system. For the mean-field model, only the study of systems with random a th is non-trivial. In addition, this model necessitates stress non-conservation in order to apply to earthquakes. We have used a trick to fix parameter a for a given number of elements (L) in the meanfield model. W h e n such a trick is used, we get a rather poor fit to the finite-size-scaling function (3.1) when r = 0.30 and v — 1.38 are used. We conjectured that this could be Chapter 3. Study of the Scaling Properties of the Model 61 due to the sensitivity of the exponent r to the details of the model. However, it is also possible that the trick that we have used to fix a does not work well. The general conclusions that we can draw from our study of the scaling properties are that the model of X u et al. displays S O C in the conservative case but not in the non-conservative case. In the conservative case, this was proven by the fact that C(S, L) is a power law which cuts off due to finite-size effects. O n the other hand, in the nonconservative case, a characteristic earthquake size smaller than the system size appears. In such a case, C(S,L) is best fitted to a function of the type (3.3) N O T self-organized critical. and the model is Chapter 4 T i m e Series Analyses In this chapter, we implement two types of analysis of time series generated by the original model of X u et al. These time series are presented in section 4.1. In section 4.2, we perform a nonlinear forecasting analysis. T h i s analysis investigates whether the data in a time series generated from a dynamical system exhibit low-dimensional chaotic behavior, as opposed to high-dimensional (or stochastic) behavior. T h e second time series analysis consists in doing a rescaled range analysis, which allows computing the Hurst exponent (see section 4.3). T h i s exponent enables one to learn about the time correlations (or memory effects) in a system. Finally, in section 4.4, we summarize the results obtained in this chapter and draw our conclusions. 4.1 T i m e Series of Interest A t least two types of time series can be generated by the model of X u et al. T h e first type involves time series consisting in the size of the earthquake sequences (S) at different values of the long time (t): {Si} (i = 1, 2 , N ) . A typical time series of this type is plotted on figure 4.1 for a 100 x 100 system with random a ht A time series looking qualitatively similar is obtained for a two-dimensional system with constant a h as well t as for one-dimensional systems (see chapter 3). (Most of the work in this chapter will be performed on time series for two-dimensional systems, but for completeness, we will sometimes analyze time series for one-dimensional systems with constant a h since we t saw in chapter 3 that stress ordering prevails in these systems. 62 Also, unless otherwise Chapter 4. Time Series 63 Analyses specified, the values of the parameters of the model of X u et al. will always be assumed to be a = 1, X = 0 and 9i = 0.25). From figure 4.1, we note that Si fluctuates greatly over time and that the sequences are instantaneous events on the long time scale. Also, even if it is not completely clear from the figure, the time intervals between consecutive earthquake sequences vary over several orders of magnitude. Therefore, a time series such as the one plotted on figure 4.1 is highly variable in both its dependent variable (Si) and independent variable (ti). A t this point, it is desirable to fix one of these variables and allow the other one to vary. T h i s can be done according to one of the three following procedures: (a) focus only on the earthquake sequences having S = So ° r S > So- T h e n , the time intervals between these consecutive sequences constitute a time series of interest, (b) focus on all the earthquake sequences whatever S. T h e n , the time intervals between these consecutive sequences constitute another time series of interest, (c) sample the original time series every At time units. Now, because there is in general not a sequence happening at a particular sampling time, the sequence which is closest in time to the latter will be treated as having happened at that sampling time. T h i s sequence will have its size retained as a datum for a newly generated time series consisting in the size of earthquake sequences about equally sampled in time. (For completeness, we mention that another related procedure that takes advantage of the fact that the earthquake sequences are point events could be used: sum the size of all the earthquakes occurring in a time window of a given size. T h e time corresponding to the midpoint of the window, for instance, can then have this sum of earthquake sizes assigned to it. T h e resulting time series is also evenly sampled. T h i s procedure is appealing but we will not consider it in this thesis). Chapter 4. Time Series 64 Analyses Consequently, the time series of type 1 exist in three varieties, defined as varieties (a), (b) and (c). A n example of a time series of variety (b) is plotted on figure 4.2. It can be obtained from figure 4.1 by calculating the time intervals between the consecutive earthquake sequences. A s mentioned above, it is obvious from figure 4.2 that these time intervals vary over several orders of magnitude. Now, regarding the time series of variety (a), their structure depends on the value of the parameter So. W h e n So is small, the time intervals between the consecutive sequences are rather small whereas for So large, the time intervals are larger (results not shown). T o this effect, we will try to confirm by a different technique the results obtained by Chen et al. [16] with the model of X u et al.: the large earthquakes are temporally correlated but the small earthquakes are not. Finally, for the time series of variety (c), we noted that parameter At has little effect on their structure (time series not shown). T h e second type of time series (hereafter called type 2) which can be obtained from the model of X u et al. involves the average stress in the system (a) as a function of the long time (t), where a is defined as N s being the number of elements in the system and o\ the stress on element i. Here, we have chosen not to take the absolute value of o~i in calculating a. A s we have checked, this choice has no effect on the time series for one-dimensional systems since the cr^'s are all positive. O n the other hand, for two-dimensional systems, <7; is sometimes negative, but we have checked that this occurs for about 1% of the elements only and so, the effect on a is minor. ah t A time series of this type for a 100 x 100 system with random is plotted on figure 4.3. A time series looking qualitatively similar is obtained for a two-dimensional system with constant o . as well as for one-dimensional systems. tn From figure 4.3, we note that a well-defined pattern appears: diagonal lines are connected to Chapter 4. Time Series Analyses 65 450-r 400- 350- 300- 250- •50 J i 0 i 0.2 r 0.4 ' i 0.6 i 0.8 i 1 • i 1.2 i 1.4 i 1.6 U Figure 4.1: Typical time series of type 1: Size of an earthquake sequences (Si) as a function of the occurrence time (tj). The whole time series has 5 x 10 data but for clarity, only the first 2000 are plotted. The sequences were generated from a 100 x 100 system with random o h- The time at which the first sequence happens is arbitrarily set to«i=-0. 4 t Chapter 4. Time Series Analyses 0.005- 0.004- - f o * 0.003- + 0.002- 0.001 - 0- 1 1 1 1 1 1 1 1 1 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 Figure 4.2: Typical time series of type 1 [variety (b)]: Time intervals £j — (i = 2 , 3 , 2 0 0 0 ) between consecutive earthquake sequences. The sequences were gen erated from a 100 x 100 system with random c h (see figure 4.1). t Chapter 4. Time Series 67 Analyses vertical lines, which are themselves connected to diagonal lines, and so on. T h e diagonal lines symbolize average stress increases due to system driving. T h e slope of the diagonal lines (i.e. the rate of stress increase) is symbolized by p and was set to 1 for convenience (see section 2.4). T h e vertical lines represent average stress decreases due to earthquake occurrences and give a hint about the magnitude of the earthquakes. A time series such as the one on figure 4.3 is unevenly sampled. We will also consider below evenly sampled versions of the time series of type 2 by sampling every At time units the original time series of type 2. T h i s is possible since, as explained in section 2.4, we assume that a increases linearly with time between two earthquake sequences. As a preliminary study, we can analyze the time series plotted partially on figures 4.2 and 4.3 by considering the linear correlation coefficient [85] ^2(xi - x){x i+i - x') where / is called the lag and x and x' are the mean of the time series {XJ} and (i = 1,2,...,N) respectively. A s suggested in [54], r is reliable only for I < N/4 and thus, we do not calculate 77 for I > N/4. ; T h e set of values r (I = 1, 2 , N / 4 ) and the ( plot of 77 against I is referred to as the Autocorrelation Function ( A C F ) . T h e A C F for the time series plotted on figure 4.2 and the time series plotted partially on figure 4.3 are displayed on figures 4.4 and 4.5 respectively. T h e A C F for the time series of type 1 (N = 1999) exhibits oscillations around 0, and 77 remains small for all / except / = 0 (ro = 1 is not shown on figure 4.4). T h i s seems to reveal that the data in this time series are not correlated. O n the other hand, the A C F for the time series of type 2 (N = 2000) decreases with / (for I < 100) and fluctuates around 0 at larger It is of interest to note that the bump of the A C F between / « 160 and I « 220 (see figure 4.5) is a feature reflecting the presence of non-trivial correlations in the time series of type 2. However, Chapter 4. Time Series Analyses 68 0.35 0.325- 0.32 — i 1 0.05 J 0 1 0.1 t Figure 4.3: Typical time series of type 2: Average time (t). The data were generated from a 100 x 100 time series has 5 x 10 data but for clarity, only the the vertical lines) are plotted. The time of the first 4 1 1 0.15 0.2 stress (a) as a function of the long system with random a h- The whole first 500 data (i.e. the end points of datum is arbitrarily set to zero. t Chapter 4. Time Series Analyses 69 the A C F is not powerful enough to reveal exactly what type of correlations are present in the latter time series. More sophisticated analyses are required to do so. T w o of them are considered in this chapter and are outlined below. In summary, we consider in this chapter two types of time series. Note that only the time series of type 2 are explicitly linked to the state variable (i.e. the stress) of the model. T i m e series consisting in the value of the stress at a given location in the system as a function of the long time could have also been considered but, as we have checked, they exhibit large fluctuations and depend on the location chosen to measure the stress, so they are not adequate for a time series analysis. 4.2 Nonlinear Forecasting Analysis In this section, we implement the first of two types of time series analysis, namely a nonlinear forecasting analysis. We begin with a general introductory discussion in subsection 4.2.1. Afterwards, we describe the two important ingredients entering into a nonlinear forecasting algorithm, namely the state space reconstruction (subsection 4.2.2) and the nonlinear function approximation (subsection 4.2.3). T h e forecasting algorithm that we use is presented in subsection 4.2.4 and tested by considering time series having a known structure in subsection 4.2.5. In subsection 4.2.6, we apply this algorithm to the analysis of the two types of time series generated by the model of X u et al. and discuss the results obtained. 4.2.1 Introductory discussion Until recently, it was assumed that randomness was caused by extreme complication, i.e. the presence of many irreducible Degrees O f Freedom ( D O F ) . T h i s led to Kolmogorov's theory of random processes, which he defined in terms of the joint probability distribution 70 Chapter 4. Time Series Analyses I Figure 4.4: Autocorrelation function for the time series of type 1 [variety (b)] plotted on figure 4.2. For convenience, the value r = 1 is not shown on the graph. 0 Chapter 4. Time Series Analyses 71 Chapter P. 4. Time Series For a time series 72 Analyses (i = 1 , 2 , N ) , the d ih '-P(fi. •••>&) = Probability^ order probability distribution is < • < T h e theory of random processes copes with situations in which inadequate information is available and makes no statements about the causes of the randomness. A s a result, with more data or more accurate observations, a phenomenon that previously seemed random might become more predictable, and hence less random. M a n y of the classic examples of randomness are not complicated. T h i s is the case for a flipping coin for example. The dynamics of a flipping coin involves only a few D O F . Its randomness comes from a sensitive dependence on the initial conditions: a small perturbation causes a much larger effect at a later time, making prediction difficult. W h e n sensitive dependence on the initial conditions occurs in a sustained way, it is called chaos (see [21] for an outline at the introductory level). Since chaos is defined in the context of deterministic dynamics, in some strict sense chaos is not random. B u t it is relevant to note that chaos amplifies noise exponentially, so any uncertainty present initially is amplified to macroscopic proportions in a finite time. A s a result, short-term determinism becomes long-term randomness. The flipping coin being a chaotic system, it is inherently nonlinear. A nonlinear system is not necessarily chaotic though! T h e flipping coin is an autonomous system since it receives no time dependent input from the outside. Non-autonomous (or driven) systems are rather common: they are found in many fields of science, for example in engineering, chemistry, biology, physical geography and economics. A n external agent provides the necessary drive, which can take on many different forms: periodic and random are the most common forms. W h e n the driving process is periodic, a Poincare section can reduce the system to an autonomous mapping problem. O n the other hand, when the driving process is random, it is not possible to simplify the problem in such a way. T h i s is the case in vibration tests performed on mechanical systems. A generic representation of Chapter 4. Time Series 73 Analyses the kind of mechanical systems that are studied in engineering is provided by the driven Duffing oscillator. It is described by the equation (see [47] for instance) x"~x " + 2(uj (x'-x ')+u (x-x ) + aLo (x-x ) 2 0 n 0 n 2 0 n 2 0 + (3u 2 n (x - x^) \ X-X Q | = 0 , (4.1) where u> is the natural frequency in the absence of nonlinearity, £ the damping coefficient n and f3 (respectively a) the coefficient of the antisymmetric (respectively symmetric) nonlinear restoring force. T h e oscillator is driven through the acceleration x ", from which 0 XQ and XQ can be computed. Substantial high-order harmonic responses alert to the presence of a nonlinearity in the system. For given values of a and (3, the behavior of equation (4.1) is chaotic for appreciable driving levels and nearly linear for low driving levels [47]. T h e driven Duffing oscillator is useful to illustrate what is going on in the model of chapter 2, with the important difference that the former has only a few D O F . A s seen in section 2.4, the systems studied with the model of X u et al. are slowly driven and thus the amount of stress added to the system is just enough to cause one rupture in it. A s a result, the systems studied with the model of X u et al. will not be treated as input-output systems: we will consider instead that only a single time series is available for the forecasting analysis and use the techniques appropriate for such situations (see below). Note that the important role played by the external drive in the model of X u et al. should not be minimized though. B o t h chaotic and random time series have broadband spectra, so that distinguishing them requires some other statistics than the autocorrelation function and the Fourier power spectrum. There are a number of statistics which are sensitive to nonlinear behavior, namely the short-term prediction error, the Lyapunov exponents [112] and the correlation dimension [39]. T h e short-term prediction error is small for a chaotic system and large for a random process. T h e Lyapunov exponents quantify the rate at which Chapter 4. Time Series Analyses 74 nearby trajectories in the state space of a system diverge. For a system to be consid- ered as chaotic, at least one of its Lyapunov exponents must be positive (necessary but not sufficient condition). T h e correlation dimension (v) estimates the dimension of the subset of the state space of a system (often called the attractor of the system) that the trajectories approach. For example, fluid flows have an infinite-dimensional state space, but can have low-dimensional attractors [60]. T h e correlation dimension can be obtained by applying the algorithm of Grassberger and Procaccia [39]: v n where C n limlogC„(e)/log(e), is the correlation integral for the n-histories of the time series »W = (iv(FrT)) g ^ " II * " C 9(x) = '"' being the Heaviside function (9(x) = 1 for x > 0 and 9(x) = 0 otherwise) and the norm used is the maximum norm. T h e n-histories of the time series are constructed by using the vectors •^i — ip^i—n+l j %i—n+2•> • • • > %i)• If the time series comes from a chaotic system, then v n does not increase with n and converges to the value v. A n alternative way of characterizing the attractor of a system is through the fractal (or Hausdorf) dimension (D). T h i s number can be obtained by using the box-counting algorithm for instance (see [40]). In most cases, if D < 2, then v and D are close together, but it is important to keep in mind that D has only to do with the geometrical structure of the attractor whereas v is sensitive to the dynamical process of coverage of the attractor [39]. A s a result, u tends to be a more relevant measure of the attractor than D. T h e correlation dimension (y) is useful to distinguish between chaos and randomness. However, the use of the correlation dimension as a statistic to distinguish between chaos and randomness has its limitations (see for instance [79]). Chapter 4. Time Series 75 Analyses Along the same line, three economists (Brock, Dechert and Scheinkman) have recently introduced a test based on the correlation dimension: the B D S test [9]. It is a test against IID (Independent Identically Distributed) noise (null hypothesis). The test can be modified to single out the nonlinear behavior by subtracting out the linear component of the data and applying the B D S test to the residuals. A procedure advocated in [102] is instead to combine the B D S test with a technique called bootstrapping data sets. with surrogate The result is a powerful tool to detect any nonlinearity in a time series. In this work, we consider only the short-term prediction error. This statistic is very powerful in distinguishing between low-dimensional chaos and randomness, but does not have the limitations observed with the correlation dimension. A n example of the use of the short-term prediction error for such a purpose is provided by the nonlinear forecasting algorithm of Farmer and Sidorowich [29]. In this algorithm, the error of a local linear predictor is plotted against a variable smoothing parameter. If maximum smoothing allows to get the best forecasts, then the structure in the time series is described by a linear stochastic model. If the smallest forecasting errors are obtained for minimum smoothing, then the time series is generated from a low-dimensional chaotic system. If any amount of smoothing less than the maximum but more than the minimum is needed to get the best forecasts, then a nonlinear stochastic model describes the structure in the time series. We will outline in subsection 4.2.4 an algorithm which can be used to construct these different types of models. For now, we mention that linear stochastic models have a long history, going back at least to Yule [116] and more recently Box and Jenkins [8]. On the other hand, nonlinear deterministic models have been constructed by the dynamical systems community since about 1987 (see [23] for instance) and independently, the statistics community has constructed nonlinear stochastic models since about 1980 (for a review see [103]). In this work, we assume that a good model from first-principles can not be found, Chapter 4. Time Series 76 Analyses and so we must find it directly from the time series obtained from the observation of the dynamical system. This situation is what happens for a dynamical system studied with the model of X u et al. since the collective behavior of the large number of elements in the system does not arise from a small set of equations. (Generally, self-organized critical systems are believed to have many D O F [14]). The time series generated in this thesis are all scalar (see section 4.1), i.e. they each involve only one variable, which is less than the number of variables necessary to fully describe the dynamical system. These time series are not contaminated by observational noise. However, dynamical noise affects the time series because the dynamics is in general not deterministic. To this effect, it is relevant to note that because time series for systems with both random o h and constant o h are t t considered, then it is possible to study the effect of the dynamical noise on the system dynamics. In any case, if neither observational noise nor dynamical noise are too large and the randomness in the time series is caused by low-dimensional chaotic behavior, then building a dynamical model directly from the time series involves two steps: state space reconstruction and nonlinear function approximation. These are the two important ingredients entering into a nonlinear forecasting algorithm. Before we present the version of the nonlinear forecasting algorithm that we use, we describe in some details these two ingredients. 4.2.2 State space reconstruction In most physical situations, only a single scalar time series {xi} (i = 1 , N ) obtained by making measurements on a dynamical system is available. The time series has in general fewer variables than the number of variables necessary to fully describe the system. A n illustration of this situation is provided by a fluid flow experiment, which in principle can be modeled by the Navier-Stokes equations. In the case in which only a single measure of a given velocity component at a given point in space is performed, then the resulting Chapter 4. Time Series 77 Analyses time series is inadequate to provide the initial conditions for the Navier-Stokes equations. To build a model directly from the time series, the state space of the dynamical system must be reconstructed from this unique time series. A popular method for reconstructing the state space of a dynamical system was introduced by Packard et al. [80] and put on a firm mathematical foundation by Takens [101]. In this method, the past behavior of a single scalar time series is used to get information about the present state of the system. This information can be represented as a delay vector of embedding dimension m and delay time r: X-i — (Xj, X j _ , r Xj_( _i) )^, m T where the components of a delay vector are called the delay coordinates. (4-2) The t symbol insures that the delay vectors are column vectors, which is the standard way of representing them in the literature. The optimal choice of parameters m and r has been the subject of many theoretical studies (see [34] for instance). In practice, r is arbitrary but it should not be too small or too large. Also, m is often chosen by trial and error, starting with a small value and increasing it, searching for optimal results. Other methods of state space reconstruction of common use include the principal component technique and reconstruction using derivative coordinates (see [30] for a brief outline of these methods). Also, filtering is sometimes used in combination with any reconstruction method. It is of interest to note that principal component, derivative and delay coordinates, and filtered versions thereof, are all related to each other by linear transformations whereas the transformation from delay coordinates to the original coordinates (i.e. the data of the time series) is nonlinear. As demonstrated by Fraser [33], nonlinear coordinate transformations can be greatly superior. It is not clear which method of state space reconstruction works best and what causes one method to be better than another though. In the forecasting algorithm presented in subsection 4.2.4, delay Chapter 4. Time Series 78 Analyses coordinates are used to reconstruct the state space of the dynamical system. 4.2.3 Nonlinear function approximation Once the state space of the dynamical system has been reconstructed, the next task is to predict the future behavior of the scalar time series obtained by making measurements on this system. This is called forecasting and it involves an extrapolation in time, in the sense that one uses data from one domain (the past) to extrapolate to the behavior of data in a disjoint domain (the future). State space reconstruction makes it possible to convert the extrapolation of a time series into a problem of interpolation in the state space of the system. This problem can be expressed mathematically in a simple way. Indeed, using the fact that the delay vectors [see (4.2)] contain information about the states of the system at different times, we have X = f Xi T i+T = f {x , T l £ _( _i) ) , Xi- , (4.3) t i T m T where, for simplicity, we assume that only the first coordinate of X{ T is forecasted, + T being called the forecasting function f T time. The function f T : 7c m —> TZ holds for all i. The can, in general, only be determined approximately. The most commonly used approach to time series forecasting assumes that f T is a linear function [8]. However, almost any real system contains nonlinearities, which must be taken care of by a nonlinear f. T A particular and important kind of nonlinearity is one due to chaos, for which the scalar time series is generated by motion on a L>-dimensional attractor. Sauer et al. [87], it is sufficient to take m > 2D for a smooth f T According to to hold for all i and if D < m < 2D, then (4.3) holds for almost all i. In contrast, if the time series is generated by a stochastic process then, for all m, a noise term is expected to appear on the right-hand side of (4.3) and the delay vectors are expected to fill out a set of dimension m [87]. This situation also prevails if the time series is generated by motion Chapter 4. Time Series 79 Analyses on an attractor with dimension D large. Indeed, in such a case, there will be insufficient data to approximate f T with a deterministic model with m > D and in this sense, a high- dimensional deterministic system is equivalent to a low-dimensional nonlinear stochastic system [14]. On the other hand, when D is small, then an accurate approximation of f T can always be obtained with a modest number of data (N ~ 10°). The approximation of function f T requires that a representation as well as a method of selecting the parameters of the representation be chosen. When f T is nonlinear, there is an infinite variety of possible representations, and in the absence of any prior knowledge, no reason to prefer one representation over another. In many cases, a given representation is chosen because it allows a fast algorithm for fitting the parameters. It is beyond the scope of this thesis to review the most popular nonlinear function representations (see however [30]), but it is relevant to mention that representations can be both global and local. We describe next one one way of obtaining a local approximation of f since T this method is used in the nonlinear forecasting algorithm considered below. A local approximation usually enables one to get a good fit since the variations are less extreme than when a global representation is used. The basic idea of a local approximation of f T is to break up the domain of f T into local neighborhoods and fit different parameters in each one. There are several ways to assign the neighborhoods (see [30] for a review). One way of doing this was introduced by Farmer and Sidorowich [29] and is now commonly used: they impose a metric on the state space, denoted by || ||, and find the k nearest neighbors of Xi [see (4.3)], i.e. the k delay vectors Xj with j < i that minimize || Xi — Xj ||. Afterwards, in each neighborhood, a local predictor for x T i+ is constructed. The simplest approach to do so is an approximation by the nearest neighbor (k — 1), i.e. X?™T — j+T- X For example, to predict tomorrow's weather, one would search the historical record and find the weather pattern the most similar to that of the last few days, and predict that Chapter 4. Time Series Analyses 80 tomorrow's weather pattern will be the same as the neighboring pattern one day later. A superior approach is the local linear approximation in fitting a linear polynomial to the pairs ( X J , X J + T ) . (k > 1). T h i s approach consists T h e fit can be made in any of several ways. Often, least squares by singular-value decomposition (see [85] for instance) is realized. In principle, approximations using higher-order polynomials can be performed but in practice, this may or may not give an improvement depending on the time series analyzed. T h e method of Farmer and Sidorowich has the advantage of being quick and accurate, but the disadvantage that the approximations are discontinuous. Continuity may be guaranteed by weighted local fitting (see an illustration of this in a particular case in [98]) but for simplicity, simple local fitting is often considered. In the nonlinear forecasting algorithm presented next, simple local linear fitting is implemented. 4.2.4 Nonlinear forecasting algorithm considered T h e nonlinear forecasting algorithm presented below investigates whether the data in a scalar time series generated from a dynamical system exhibit a low-dimensional chaotic behavior, as opposed to a high-dimensional or stochastic behavior. T h e algorithm constructs models with a variable smoothing parameter (k) which at one extreme defines a nonlinear deterministic model, and at the other extreme, a linear stochastic model. T h e accuracy of the short-term forecasting error as a function of the smoothing parameter reveals much about the underlying dynamics generating the aperiodic time series. In particular, if any amount of smoothing less than the maximum allows one to get the best forecasts, then the structure in the time series is described by a nonlinear stochastic model. T h e smoothing parameter (k) represents the number of nearest-neighbor delay vectors included in the local linear fitting of the function f T (see the discussion in the previous subsection). In addition to its use for detecting low-dimensional chaos, the nonlinear forecasting Chapter 4. Time Series Analyses 81 algorithm considered below can be used to investigate whether time series from highdimensional deterministic systems can be approximately modeled with low-dimensional nonlinear stochastic models. This problem has recently attracted the attention of some physicists who study collective phenomena such as spatiotemporal chaos [22], robust intermittency [53] and Self-Organized Criticality (SOC) [3]. A version of a nonlinear forecasting algorithm is now outlined. It was mainly inspired from the algorithm introduced by Casdagli [14] and described by Ward [109]. The implementation of this algorithm involves the following steps: (a) Divide the analyzed time series into a fitting set X N / + I , xi,...,XN F and a testing set ...x . As suggested in [102], the fitting set should be appreciably larger N than the testing set. (b) Choose an embedding dimension m and set the delay time r as well as the forecasting time T to 1. (c) Choose a delay vector Xi with Nf + m < i < N — 1 for a 1-step ahead forecasting test. (d) Compute the distances dij of test vector Xi from the delay vectors Xj (m < j < Nf — 1) in the fitting set using the ordinary Euclidean norm. (e) Order the distances d^, find the k nearest neighbors Xj(i), ^j{2), ^j(k) of Xi and fit a model of the following form: m w a + d x iy i 0 n j{ n+ (I = 1, ...k). (4.4) 71 = 1 The heapsort algorithm described in Press et al. [85] is used to order the distances dij. The parameters exo,a m are obtained by solving the normal equations for the linear system of equations (4.4) using singular-value decomposition (see [85]). Chapter 4. 82 Time Series Analyses (f) Use the fitted model (4.4) to estimate a 1-step ahead forecast Xi+i(k) for test vector Xi and compute its error ei(k) - \ x i(k) i+ - x i+l |. (g) Repeat steps (d)-(f) for all allowed i in the testing set and compute the normalized Root-Mean-Square (RMS) forecasting error E {k) m = 1 Ntxo- where a is the standard deviation of the time series and N = N — Nf — m is the t number of allowed vectors in the testing set. The time series is forecastable only if E (k) m < 1. (h) Repeat step (g) for several representative values of k between 2(m + 1) and Nf — m. (i) Vary the embedding dimension m. This algorithm can readily be applied to the analysis of time series by studying the curves E (k) as a function of k for a few values of m. We will do that first for time series m having a known structure in order to illustrate the different types of curves that result from the analysis. Afterwards, we will apply this algorithm to the analysis of the time series generated by the model of X u et al. (see section 4.1). In the literature, it is not always clear whether a nonlinear forecasting algorithm should be applied to the raw time series or the first-differenced time series. On the one hand, Sugihara and May [99] have analyzed first-differenced time series, where the first-difference operation was performed, they say, to clarify the suspected nonlinearities by reducing the effects of any short-term linear autocorrelations. On the other hand, Casdagli [14] has analyzed raw time series that he assumed stationary, this assumption 83 Chapter 4. Time Series Analyses being important since only stationary time series should be analyzed with the algorithm. According to Kendall and Ord [54], a time series is stationary if its mean and variance are constant over time and if the structure of the time series depends only upon the relative position in time of the observations. These conditions must be satisfied at least to a reasonable degree of approximation for a time series to qualify as stationary. In the case of the time series generated by the model of X u et al., we have checked that they are reasonably stationary. This means that it is not necessary to first-difference the time series since they are stationary in the mean (i.e. they do not exhibit a linear drift in the mean). However, the first-difference operation performed in [99] was shown to improve the results. For this reason, we will take no chance and consider both raw time series and their first-differenced counterparts. This obviously applies only to the time series generated by the model of X u et al. since the ones considered in the next subsection have a known structure and we therefore know what results to expect from the algorithm. 4.2.5 Tests of the algorithm with time series having a known structure In this subsection, we test the forecasting algorithm outlined in subsection 4.2.4 by analyzing four time series having N = 1000 data and a known structure. Before we do this, we describe how we generate these time series. The first two time series are generated by low-dimensional chaotic maps, namely the logistic map and the Henon map. The logistic map is univariate and is defined as [65] where p is the parameter of the map. For p > 3.58, the map exhibits chaos for almost all values of p. The time series was obtained for the value p = 4.0, a value in the chaotic regime. The Henon map is bivariate and is described by the pair of equations [44] Xi = 1 - Ax\_ + x yi-i 84 Chapter 4. Time Series Analyses which can be combined to get the single equation . Xi = 1 - Ax}_ x 4- Bxi-2- For the special values A = 1.4 and B = 0.3, the map is known to have a chaotic behavior. The last two time series have a substantial stochastic component. One of them is generated by a nonlinear stochastic model, the Threshold AutoRegressive (TAR) model. The T A R model is expressed mathematically as [104] Xi = ax^i + Ui (XJ_I < 1) Xi = fiXi-i + Ui > 1), where Ui is distributed according to a Gaussian with zero-mean and standard deviation equal to 1 (see [85] for the description of an algorithm to produce Gaussian deviates). To obtain the time series, we used the values a = —0.4 and 0 — 0.5. The last time series is generated by a linear stochastic model, a AR(1) (first-order autoregressive) model (see [54] for instance), i.e. Xi = (f)Xi-i + Ui, where <f> is a constant. To obtain the time series, we used the value <> / = 0.5. For the logistic map, the T A R model and the AR(1) model, it is necessary to choose one starting point for the iterations, whereas for the Henon map, it is necessary to choose two. The starting point(s) has (have) some effect on the resulting time series. However, if, say, the first 500 data generated in an iteration process are discarded, then the resulting time series will not be sensitive to the starting point(s). The four time series can then be readily analyzed using the nonlinear forecasting algorithm presented in the previous subsection. Chapter 4. Time Series 85 Analyses The two time series generated by low-dimensional chaotic maps were first analyzed with the forecasting algorithm and the results are displayed on figure 4.6 for an embedding dimension (m) equal to 2. Comparable results are obtained for m > 2 and are not shown. As expected, the best forecasts (i.e. the smallest E (k)) m are obtained, for both time series, when the smoothing parameter (k) is small. It is relevant to note also from the figure that excellent forecasts (i.e. E (k) m -C 1) are obtained at small k and that they are much better than the ones obtained at large k. This reflects the fact that no stochastic component is present in the time series. It is also possible from figure 4.6 to compare the forecasting errors for the logistic map and the Henon map. Indeed, at small k, E (k) m is always smaller for the logistic map than for the Henon map. This feature reflects the fact that the attractor for the Henon map has a more complex structure than the one for the logistic map. If we now turn to the time series having a stochastic component, the results of the analysis are plotted on figure 4.7 (time series for the T A R model) and figure 4.8 (time series for the AR(1) model) for three values of m in both cases. On figure 4.7, we see that the best forecasts are obtained at intermediate values of k (k ~ 400 — 500) and that the improvement in accuracy over a linear stochastic hypothesis (corresponding to the largest values of k) is only of about 2 — 3%. This is typical of a nonlinear time series having a substantial stochastic component. On figure 4.8, we notice that the forecasts get better as k increases, but no upturn in E (k) m at large k is observed. These features are typical of a linear time series having a substantial stochastic component. Finally, it is interesting to compare figure 4.6 to figures 4.7 and 4.8. We note that E (k) m is always much smaller for the time series generated by the low-dimensional chaotic maps than for the stochastic time series, as it should. Now that we have obtained typical curves of E (k) m versus k for time series having a known structure, we turn to the analysis of the time series introduced in section 4.1 and Chapter 4. Time Series Analyses 86 1 0 100 200 k 300 400 500 600 Figure 4.6: Normalized RMS forecasting error E (k) as a function of k for time series having N = 1000 data generated by the logistic map (with p = 4.0) and the Ffenon map (with A = 1.4 and B = 0.3). The solid line is for the logistic map time series and the dashed line for the Henon map time series. The parameters of the nonlinear forecasting algorithm are Nj = 800 and m = 2. m Chapter 4. Time Series Analyses 87 0.96 0.95- 0.89-—i 0 1 1 100 1 200 1 300 r 400 k 1 500 1 600 r700 Figure 4.7: Normalized RMS forecasting error E (k) as a function of k for a time series having N = 1000 data generated by the T A R model (with a = -0.4 and 0 = 0.5). The parameters of the nonlinear forecasting algorithm are Nf = 800 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line). m Chapter 4. Time Series Analyses 88 Figure 4.8: Normalized RMS forecasting error E (k) as a function of k for a time series having N = 1000 data generated by the AR(1) model (with <f> = 0.5). The parameters of the algorithm are Nf = 800 and m = 2 (solid line), m — 3 (dashed line), m — 5 (dotted line). m Chapter 4. Time Series 89 Analyses generated by the earthquake model of X u et al. 4.2.6 Results obtained w i t h the time series of interest Unless otherwise specified, all the time series considered in the analysis have N — 2000 data. As discussed in subsection 4.2.4, we consider both raw time series and first- differenced versions of them. The two types of time series outlined in section 4.1 are analyzed. Also, for a given type of time series, time series for systems with both random Oth and constant o h are considered. t We first analyze the time series of type 1. As discussed in section 4.1, the time series of type 1 exist in three varieties, namely (a), (b) and (c). We have analyzed the three varieties of time series of type 1 and found, for a 100 x 100 system, that only the time series of variety (c) are forecastable (i.e. have E (k) m < 1). The results for a time series of variety (c) for a system with random o h are plotted on figure 4.9. Similar curves are t obtained with a time series for a 100 x 100 system with constant o h (results not shown). t On figure 4.9, we observe a monotonic decrease of E [k) m with k. Also, we notice that the smallest forecasting errors are at most 5% better than the standard deviation of the corresponding time series. This suggests that the time series of variety (c) are weakly forecastable, probably because they have a large stochastic component. In addition, no nonlinear structure exists in the time series of variety (c). This was also found to be true for the time series of varieties (a) and (b) and whether or not the time series are first-differenced. Now, we turn to the time series of type 1 for a one-dimensional system. Contrary to the case of the time series for a 100 x 100 system, all three varieties of time series are forecastable, but they still have a large stochastic component and no nonlinear structure (results not shown). A time series which is of particular interest is the one consisting in the time intervals between large earthquakes. As we saw in chapter 3, the one-dimensional systems with constant o h are particularly interesting since they display t Chapter 4. Time Series 90 Analyses stress ordering. For this reason, we included only the results for a time series of variety (a) (S > 1000) for a 1000 x 1 system with constant a h (see figure 4.10). The curve for t m = 2 in particular seems to show that the best forecasts are obtained at intermediate values of k, but it should not be concluded that the time series has a nonlinear stochastic structure since it has not enough data. Instead, we prefer to safely state that the time series has probably no nonlinear structure but a large stochastic component. The best forecasting errors are 25% better than the standard deviation of the data. This may not be sufficient to allow the prediction of a future large earthquake in the model of X u et al. To further confirm this possibility, we have plotted on figure 4.11 the histogram N(5 ) of the time intervals 5 between large earthquakes. t t The plot on figure 4.11 is a log-linear one in order to better illustrate that the distribution is roughly exponential, i.e. N(5 ) t ~ e~ / °. 5t 5 Obviously, because of the small number of data used to get the histogram, then the straight line on figure 4.11 does not cover the whole range of 5 . We t can then conclude from figure 4.11 that the large events generated in the model seem to be Poisson events. (As a result, we have confirmed explicitly that it was correct above to conclude from the curves on figure 4.10 that no nonlinear structure was present in the data). We got similar results with the time intervals between large events generated in a 1000 x 1 system with random a h as well as in 100 x 100 systems. This result is t surprising and seems to contradict a conclusion of Chen et al. [16] obtained in a study of the model of X u et al. that the large events are temporally correlated. However, we believe that our conclusion is not in contradiction with theirs simply because we studied a series version of the model and they studied a parallel version of it (see the discussion in section 2.4 and [35]). Consequently, it seems that only the more complicated to study parallel models can give rise to temporally correlated large events. We next turn to the analysis of the time series of type 2, which are unevenly sampled. The E (k) versus k curves obtained for a 100 x 100 system with random a h are plotted m t Chapter 4. Time Series Analyses 91 0.93- 0.92 * i 0 - 1 100 1 200 1 300 1 400 1 500 k 1 600 H 700 Figure 4.9: Normalized RMS forecasting error E (k) as a function of A; for a time series of type 1 [variety (c)] for a 100 x 100 system with random o h- The sampling time interval A t is equal to 10 times the average time interval between the earthquake sequences. The parameters of the nonlinear forecasting algorithm are Nf = 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line). m t Chapter 4. Time Series Analyses 92 Figure 4.10: Normalized RMS forecasting error E (k) as a function of k for a time series of type 1 [variety (a)] for a 1000 x 1 system with constant o h- The data in the time series (N = 331) are the time intervals between large earthquakes (S > 1000). The parameters of the nonlinear forecasting algorithm are Nf = 264 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line). m t Chapter 4. Time Series Analyses 93 100 Figure 4.11: Histogram plot of the time intervals between large earthquake sequences (S > 1000) generated in a 1000 x 1 system with constant u h- The data set has 331 data and 20 bins were used to get the histogram. The vertical axis is logarithmic in order to better illustrate that the distribution is exponential. t Chapter 4. Time Series 94 Analyses on figure 4.12 (raw time series) and figure 4.13 (first-differenced time series). The curves obtained on both figures reveal a weak nonlinear stochastic structure (see figure 4.7 for comparison). However, E (k) m is appreciably smaller and the curves are less noisy for the raw time series, showing that the first-difference operation has blown up the stochastic component. It is therefore disadvantageous to perform first-differencing on the time series of type 2. On figure 4.12, the improvement in predictive accuracy of the nonlinear stochastic models (intermediate k) over the linear stochastic models (large k) is at most of 4% and it is at most of 3% on figure 4.13. As a matter of comparison, we show on figure 4.14 the E (k) m versus k curves for a time series for a 100 x 100 system with constant o h (without first-differencing). Again, we see that the time series has a t weak nonlinear stochastic structure, where the improvement in predictive accuracy of the nonlinear stochastic models (intermediate k) over the linear stochastic models (large k) is about 8 — 9%. We believe that this 8 — 9% improvement for the time series for a system with constant o h as compared to at most 4% for the time series for a system t with random o h reflects the fact that the latter system has a more noisy dynamics than t the former. We have also considered a time series for a system with random o th when the parameter a controlling the level of stress conservation (see chapter 3) in the system is set to 0.9 (rather than 1.0 for the time series of figure 4.12). The results obtained by analyzing this time series are plotted on figure 4.15. A nonlinear stochastic structure is still observed though weaker than in the case a = 1.0 (see figure 4.12). As compared to the case a = 1.0, E (k) m is now about 30% smaller on figure 4.15, showing that the introduction of stress nonconservation in the system allows for better forecasts. We believe that this is related to the fact that, as seen in chapter 3, stress nonconservation destroys the scale invariance and introduces a characteristic earthquake size in the system. For completeness, we finally mention that the analysis of time series for one-dimensional systems leads to similar features as the ones observed with two-dimensional systems. Chapter 4. Time Series Analyses 95 Figure 4.12: Normalized RMS forecasting error E (k) as a function of A; for a time series of type 2 for a 100 x 100 system with random u h- The parameters of the nonlinear forecasting algorithm are Nf = 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line). m t Chapter 4. Time Series Analyses 96 0.87- 0.86- —i 0 1 1 100 ; —i 200 1 300 1 400 1 500 k 1 600 H 700 Figure 4.13: Normalized R M S forecasting error E (k) as a function of k for a first-differenced time series of type 2 for a 100 x 100 system with random o h- The parameters of the nonlinear forecasting algorithm are Nf = 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line). m t Chapter 4. Time Series Analyses 97 «4£ Figure 4.14: Normalized RMS forecasting error E (k) as a function of k for a time series of type 2 for a 100 x 100 system with constant o h- The parameters of the nonlinear forecasting algorithm are Nf = 1600 and m = 2 (solid line), m — 3 (dashed line), m = 5 (dotted line). m t Chapter 4. Time Series Analyses 98 0.25- 0.245- 0.24- 0.235' 0.23- 0.225' 0.22' 100 200 300 400 500 600 700 fc Figure 4.15: Normalized RMS forecasting error E (k) as a function of fc for a time series of type 2 for a 100 x 100 system with random o h when the parameter a controlling the level of stress conservation in the system (see chapter 3) is set to 0.9. The parameters of the nonlinear forecasting algorithm are Nf = 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line). m t Chapter 4. Time Series 99 Analyses The appearance of nonlinear structure on figures 4.12, 4.14 and 4.15 is connected to the fact that, when the stress releases due to earthquakes happen, then the system stops following linearly the external drive (threshold dynamics!). In addition, the presence of a stochastic component in the time series simply reflects the fact that the response of the various elements in the system varies from one element to the other: not all elements t are triggered simultaneously as one element ruptures. The main consequence of these findings is that the phenomenon of SOC that the model of X u et al. was shown to I display in chapter 3 is inherently high-dimensional. This conclusion 'might be specific i to the model studied in this thesis though. It could have been guessed from the start since the systems we study are made of a large number of elements and rarely is lowdimensional chaos found in systems with many D O F , away from the transition to chaos. This was illustrated by Casdagli [14] through a nonlinear forecasting analysis of a time series consisting in a component of the velocity at a point in a fluid in a highly turbulent regime (i.e. the Reynolds number is high). We have checked for the robustness of the nonlinear stochastic structure found above by analyzing evenly sampled time series of type 2. The results we show are for an 1 evenly sampled time series of type 2 for a 100 x 100 system with random a h, but similar t results are obtained for a system with constant a h- Figure 4.16 displays the t E (k) versus k curves for such a time series. The comparison of the magnitude of E (k) m m on figures 4.12 and 4.16 leads us to conclude that the sampling procedure has weakened i even more the nonlinear stochastic structure: the improvement in predictive accuracy of • the nonlinear stochastic models over the linear stochastic models is about 2% on figure 4.16 as compared to 4% on figure 4.12. Now, to show that the observation of a nonlinear stochastic structure in the time series of type 2 is not an artifact of the algorithm, we have randomized the phases of the frequency components of the evenly sampled time series of type 2 whose results are plotted on figure 4.16, but without altering its power Chapter 4. Time Series 100 Analyses spectrum. The resulting time series is called a surrogate data set of the original data set. We generated the latter using the method of Theiler et al. [102]. In this method, the frequency components are multiplied by a factor e *, where cp is uniformly distributed 2 between 0 and 27r. Afterwards, the phases are symmetrized to take into account the fact that the original time series is real. The E {k) versus k curves for one such surrogate data m set are plotted on figure 4.17. We note from this figure that the weak nonlinear stochastic structure apparent on figure 4.16 has disappeared and that the resulting curves look like the ones for a linear stochastic time series (see figure 4.8). As a result, we have proved that the nonlinear stochastic structure observed on figure 4.16, even if it is weak, is real and some weakly coherent phase relationships exist between the frequency components of the system. In summary, we have found that the time series of type 1 are weakly forecastable or not at all in most of the cases considered. Only the time series of the time;intervals between large earthquakes generated in a one-dimensional system was found to have some degree of forecastability. For all the time series of type 1, no nonlinear structure was found. On the other hand, the time series of type 2 were found to be weakly nonlinear but with a substantial stochastic component. These results suggest that S O C is inherently a highdimensional phenomenon. Note however that it is possible that, with a parallel version of the model (see chapter 2), this conclusion does not hold true, but this remains to be investigated. Also, we have checked that there are weakly coherent phase relationships between the frequency components. We will see in the next section that the absence of strong phase relationships between the Fourier components in the evenly sampled time series of type 2 enables us to implement a rescaled range analysis. The information that 1 we will get from this second type of analysis is complementary to the one extracted in this section. In the next section, we will restrict our attention to the time series of type 2 since, as we saw in the present section, the time series of type 1, whatever the variety, Chapter 4. Time Series Analyses 101 0.52- 0.515- 0.51 0.505- 0.5- 0.495 0.49 0.485- 0.48 100 200 300 k 400 "r300~ 500 700 Figure 4.16: Normalized R M S forecasting error E (k) as a function of k for an evenly sampled time series of type 2 for a 100 x 100 system with random o h- The sampling time interval A i is about equal to 1.5 times the average time interval between the earthquake sequences. The parameters of the nonlinear forecasting algorithm arej Nf = 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line). m t 102 Chapter 4. Time Series Analyses 0.58 j - j r - r - 0.57- j 0.56- j 0.55- i jr. 0.51- 0.5 i 0 ; 1 100 1 200 1 300 1 1 400 500 k '—|i 600 H 700 Figure 4.17: Normalized RMS forecasting error E (k) as a function of A; for a time series which has the same power spectrum as the time series of figure 4.16, jbut randomized phases. The parameters of the nonlinear forecasting algorithm are Nf —, 1600 and m = 2 (solid line), m = 3 (dashed line), m = 5 (dotted line). m Chapter 4. Time Series 103 Analyses very much look like white noise and so, are not of interest anymore. 4.3 Rescaled Range Analysis In this section, we perform a second type of time series analysis, namely a rescaled range analysis. We begin with a general introductory discussion in subsection 4.3.1, followed by a brief discussion of the fractal analyses in subsection 4.3.2. We present the algorithm used to do the rescaled range analysis in subsection 4.3.3. This algorithm is tested on synthetic time series generated by the method of successive random addition in subsection 4.3.4 and applied to the evenly sampled time series of type 2 in subsection 4.3.5. In the latter subsection, we also discuss the results obtained. 4.3.1 ; Introductory discussion The power spectrum analysis is often used to examine irregular time series. It is rather common in nature to observe power spectra which have no eminent peaks and exhibit instead appreciably more power at low frequencies than at high frequencies. In such cases, the power spectrum density P(f) is in most cases approximated as a power-law function of the frequency (/): P(f) oc where (3 is called the spectral exponent and is positive. In this approximation, (3 is one of the quantities describing the irregularity of the time series. Time series which exhibit power-law behavior in their power spectrum are often termed fractal time series [115], but as shown by Higuchi [46], this is only a sufficient condition. Indeed, it is necessary for a time series to qualify as fractal that i the phase distribution of the frequency components be random. To illustrate this more clearly, we follow the ideas of Higuchi [46]. When the total number of data (N) in the i evenly sampled time series {x(j)} (j = 1, 2 , N ) is even, x(j) can be expressed in terms Chapter 4. Time Series 104 Analyses of the Fourier series representation as follows: N/2-1 N/2 x{j) = J2 kCOs{2njk/N)+ A fc=0 N/2 Yl B sm{2njk/N) k=l k k i (k ^ 0, N/2). 1 k k Y C cos(2ijk/N-6 ), k k k=0 where C = y/A\ + B\ and 9 = tan" (B /A ) C = A . = k k k For k = 0 and k = N/2, Both 9Q and 9 / are set to zero. The variables C and; 9 represent the N 2 k k amplitude and phase for this given value of k respectively. The power [spectrum density for a given k is P(k) = Cl/N, (4.5) where the phase 9 is absent from the relation. If the time series has a power-law behavior k (with exponent (3) in the power spectrum, i.e. if P(k) oc k~^ , then (4.5) can be used to obtain C for all k except k — 0. Higuchi has fixed Co to zero. Now, to construct k simulated data having a power-law spectrum, the phase 9 must be determined for a given k k. For simplicity, Higuchi assumed that 9 is given by a random variable with a uniform k distribution in the range [0,# maa; ], where 9 correlation between the phases is set to be m a x {9 9 i) k k can vary within the range [0,360°]. The = 8 t, kk where Sij denotes the Kronecker delta. Higuchi has considered the case (3 = 2 and varied the value of 9 m a x in the range i [0,360°]. His results have shown that the resulting time series more irregular (fractal) as 9 m a x {x(j)} becomes more and increases, but definitely that for small 9 m a x , the time series is not fractal. This shows that a quantitative analysis using only the power-law i exponent (3 is inappropriate for describing the irregularity of a fractal time series: the phase distribution of the frequency components must also be random.' When (3 is in the interval 1 < (3 < 3 and 9 m a x = 360°, the time series i {x(j)} obtained by the above procedure is a fractional Brownian motion (fBm) time series. Fractional Brownian motion has been introduced by Mandelbrot and I Van Ness [62] as a generalization of the ordinary Brownian motion (corresponding to: (3 = 2). When (3 = 2, x represents the positions of a particle undergoing a one-dimensional random Chapter 4. Time Series 105 Analyses walk. Fractional Brownian motion (fBm) time series satisfy the relationship formalized by Mandelbrot and Van Ness x(hj + jo) - x(j ) = h [x(j H 0 + j) 0 - a;(jo)] for any integer h > 0 and jo, where = implies that both sides of the equation have the same distribution function. The exponent H is called the Hurst exponent (0 < H < 1) and it is related to the spectral exponent (3 through the relation j3 = 2H + 1. It is relevant to mention that self-affinity and not self-similarity applies to the fBm time series since in general the stochastic variable (x) does not have the dimension of tim'e. However, since in the literature fractal time series is often used to define both fractal ahd self-affine time series, we will stick to this expression in this thesis. Now, a property which is particular to fBm is that it has long-run correlations, in the sense that past increments are correlated with future increments. This property can be illustrated more easily by shifting the data index (j), so that j now runs from —N/2 +1 to N/2 (the data with j = 0 had previously j = N/2). Also, for convenience, the new x(0) is set to zero and all the other x(j) are shifted accordingly. Then, the correlation function of future increments x(j) — x(0) with past increments x(0) — x(—j) normalized by the variance of x(j) can be expressed as [31] c ( j ) = < - i . ( : j . y ) > = 2 — - 1 . (x(j) ) 2 As can be easily checked, C(j) vanishes for all j in the case H — 0.5. On the other hand, C(j) is non-zero independent of j when H / 0.5. This remarkable feature of fBm leads to persistence for H > 0.5 and to antipersistence for H < 0.5. Persistence means that if we have in the past a positive (or negative) increment, then we also have on average a positive (or negative) increment in the future. On the other hand, antipersistence signifies that an increasing (or decreasing) trend in the past implies a decreasing (or increasing) trend in the future. Chapter 4. Time Series 106 Analyses The generation of fBm curves requires special methods. The two most popular methods are successive random addition [107] and spectral synthesis [88]. ' It is beyond the scope of this thesis to outline one or the other of these methods. We have generated fBm curves on figure 4.18 for three values of i f by the successive random addition method (see [31] for a quick algorithm to achieve this task). We see from the figure that as H gets larger, then the fBm curves get smoother and vice versa. Indeed, the curve for H = 0.1 exhibits local noise of the same order of magnitude as the total excursions of the curve whereas the one for H — 0.9 displays rather clear trends with relatively little noise. The i curve for H = 0.5 represents a crossover between these two behaviors. As a result, we see that exponent H characterizes the degree of roughness of a fBm curve. At this point, we would like to show the relevance of the concepts introduced above for fractal times series for the evenly sampled time series of type 2. The stochastic variable of the time series of type 2, i.e. the average stress (a), is like the water level in hydrology. The accumulated differences from the mean water level is called the volume accumulation in that field. By analogy with hydrology, we can introduce the volume accumulation of the average stress in the model of X u et al. Such a quantity is plotted on figure 4.19. The curve on this figure was obtained by a running sum of an evenly sampled time series of type 2 from which we originally removed the trend and the mean. A similar curve would be obtained with an unevenly sampled time series of type 2. (Note that, for convenience, we restrict the discussion to evenly sampled time series of type 2, but since from section 4.2 we know that the original time series of type 2 and their evenly sampled counterparts have a similar structure, then this restriction does not have serious implications). The curve on figure 4.19 is reminiscent of a fBm curve (for comparison, see figure 4.18), but this does not prove that it is actually a fBm curve! We have also obtained the power spectrum for the time series plotted on figure 4.19. We now outline the procedure that we have used to achieve this task. Chapter 4. Q Time Series i 0 1000 107 Analyses ___ ^ i 2000 3000 4000 ^ • 5000 3 I \ 6000 7000 ; 8000 Figure 4.18: Fractional Brownian motion (fBm) signals having 8193 data generated by the successive random addition method. From top to bottom, the curves have H = 0.1, H = 0.5 and H = 0.9 and have been shifted for clarity. Chapter 4. Time Series 108 Analyses The discrete Fourier transform of an evenly sampled time series (with sampling time interval A i ) {x(j)} (j = 0 , 1 , N — 1) is defined as X ( k ) = j : x ( 3=0 3 ) e 2 ^ N (4.6) , where k — 0,1, 2 , N — 1. The discrete Fourier transform can be easily obtained using the Fast-Fourier Transform (FFT) technique [85]. Afterwards, the periodogram estimate of the power spectrum can be obtained at N/2 + 1 frequencies using (see [85] for instance) (fk) p where the frequencies = ^ (I X{k) | + | X(N 2 N 2 - k) |)=A | 2 (k) X (k = 0,1, 2 , A ^ / 2 ) are defined as k fk = NAt and the time series is assumed real (i.e. not complex). To prevent the leakage of certain frequencies to adjacent frequencies, it is necessary to use a technique called data windowing. The window function that we use is of the form u;0-) = - { l - c o s [ 0 - + l/2)107r/7V]} w(N - j - 1) = w(j) where the values of j considered are the ones for which the argument of the cosine function is smaller than ir (otherwise w(j) is set to 1). In other words, the window function influences only the first 10% of the data at both ends of thej time series. The way the window function enters into the calculation of the power spectrum is by replacing in (4.6) x(j) by w(j)x(j). We have obtained the power spectrum for the data set plotted Chapter 4. Time Series Analyses 109 on figure 4.19 by the above calculation (see figure 4.20). We have not performed such a calculation for the running sum of an original time series of type 2 because it is not evenly sampled. The Lomb normalized periodogram (see [85] for instance) could be used to estimate the power spectrum in the latter case, but this computation has not been accomplished here. A power-law decrease of the power spectrum is observed on figure 4.20. However, there seems to have different j3 for different frequency intervals: we fitted, for f At k < 0.04 and f At k £ [0.04, 0.3], (3 = 2.94 ± 0.05 and (3 = 3.46 ± 0 . 0 3 respectively. As a result, the running sum of the time series of type 2 has a more complicated power spectrum than a fBm time series. The change in slope on figure 4.20 is probably related to finite-size effects in time. We will return to these considerations below. Despite the fact that figures 4.19 and 4.20 seem to suggest that the evenly sampled time series of type 2 are fractal time series, we know from section 4.2 that they are not pure fractal time series since they have a nonlinear component. However, they surely have an important fractal component because the power spectrum on figure 4.20 can be fitted to two power laws and from section 4.2, we know that they have only weakly coherent phase relationships between the frequency components. The!presence of more than one component in a time series is commonly observed in nature. For example, while the human heartbeat intervals have been hypothesized to be fractal in nature (this was based solely on the observation of a power-law dependence in the power spectrum!), there are also well-defined harmonic oscillations due to the feed-back regulation of blood pressure and the external forcing by respiration [111]. Now, despite the fact that the phase relationships between the frequency components are not fully random in the evenly sampled time series of type 2, it is possible to use fractal time series, analyses, even if in principle the latter are applicable only to pure fractal time series. This is what we want to do in this section by focusing on one fractal time series analysis, the rescaled range analysis. The goal of this analysis is to calculate the Hurst exponent H, which Chapter Time 4. Series 110 Analyses 2 •3" ! 0 -i 5 1 1 10 1 15 : U Figure 4.19: Running sum of an evenly sampled time series of type 2 for a 100 x 100 system with random o h having TV = 16,384 data. The sampling time interval At is equal to 1.5 times the average time interval between the earthquake sequences. The trend and the mean have been removed from the time series before computing the running sum. t Chapter 4. Time Series 111 Analyses 0.001 0.0001- 16-05- 16-06 1e-07 le-08 16-09 1e-10- 1e-11 Figure 4.20: Periodogram estimate of the power spectrum of the running sum of an evenly sampled time series of type 2 for a 100 x 100 system with random o h- The sampling time interval of the time series At is equal to 1.5 times the average time interval between the earthquake sequences, where the time series has 16,384 data. We have removed the trend and the mean before calculating the power spectrum. To reduce the fluctuations, we have computed an average power spectrum over 8 non-overlapping subsets (with 2048 data) of the time series resulting from the running sum operation. The dashed lines were obtained by doing linear regressions and have slope (3 — 2.94 ± 0.05 (left segment) and P = 3.46 ± 0.03 (right segment). t Chapter 4. Time Series Analyses 112 characterizes the time correlations in a time series. Note that the use of the rescaled range analysis in cases where the time series has a weak fractal component is unreliable. (This shows the danger of applying blindly the rescaled range analysis to any time series because the latter will always return a value of H whatever the importance of the fractal component). In cases like that, the only technique that we think could be used is the coarse-grained spectral analysis of Yamamoto and Hughson [115]. This technique extracts the fractal component from a time series and estimates what percentage of the power spectrum is due to the fractal component (this percentage is 100 for a fBm time series and zero for a pure harmonic signal). 4.3.2 Fractal time series analyses Fractal time series analyses assume that the moments (such as the mean and variance) of the processes studied depend on the number of pieces included at a| given resolution. Because the number of pieces depends on the resolution, these moments do not have fixed, limiting values. For example, when studying the fluctuations of a fractal (or self-affine) time series, the difference between the maximum and minimum values of the time series (i.e. the range) increases with time. Therefore, as more data are analyzed, the variance increases continually. A few fractal techniques are commonly used. Among them are the dispersional analysis and the rescaled range analysis. The dispersional analysis involves the measurement of the variance of a signal at a succession of different levels of resolution [5]. This technique has been widely used in experimental physics (e.g. in molecular beam epitaxy) but also in studies related to medecine and physiology. It seems that the general physics community as well as the medical physics community have used different names to describe the same analysis. The dispersional analysis is considered as a relatively recent technique. On the other hand, the rescaled range analysis was introduced by Hurst [48] in the 1950's. Its use is still Chapter 4. Time Series 113 Analyses mostly concentrated in the earth sciences [75], but it spreads slowly to other disciplines [18, 31]. The rescaled range analysis has benefited over the years from the works of Mandelbrot and Wallis [63, 64] and Bassingthwaighte and Raymond [4] among others. Basically, this analysis studies how the range of cumulative deviations of the data from the mean (R) divided by the standard deviation of the data (S) (or the ration R/S) varies with the size of the time window (or the number of data in the window M for an evenly sampled time series). For self-affine (or fractal) time series, R/S is a power-law function of M [R/S oc M ], where H is the Hurst exponent. The rescaled range analysis H is the only fractal analysis that we consider in this thesis. 1 Before we describe the algorithm to implement the rescaled range analysis, we summarize the main ideas brought about by the pioneering work of Hurst [48]. For the Nile river, Hurst observed that records of levels at a gauge did not vary randomly, but showed series of low flow years and high flow years (i.e. memory effects!). Because the flow in a large river system such as the Nile depends on the water content in a large drainage area, the amount of water stored in the drainage area will increase in prolonged periods of higher than average precipitation and the excess amount of water stored will contribute to the discharge in drier years, and vice versa. Similar ideas were expressed later by Mandelbrot and Wallis [63]. The fBm concept is capable to take these memory effects into account. However, some authors [75] suggested that it not necessary to resort to a complicated concept such as fBm: autocorrelation (or a lack of serial; independence) is generally accepted as a feasible explanation for the existence of memory effects. 4.3.3 Rescaled range analysis a l g o r i t h m considered The rescaled range algorithms that are usually found in the literature proceeds by dividing the analyzed time series into subseries of length M: for each subseries, R(M) and S(M) l are computed. The manner in which the time series is divided into subseries, and the Chapter 4. Time Series 114 Analyses range of values of M over which the calculation is performed, control the determination of the Hurst exponent H. The procedures advocated in [108] are flawed since the subseries overlap for any M, and as a result, the values of R(M)/S(M) are not independent. The only way to achieve the independence of the calculated ratio values is to choose, for a given M, nonoverlapping subseries. This procedure was advocated in [64, 4] and used in [18]. It will be used in this work. For a given M, it is recommended in [4] to obtain an average value of Rf S over the various subseries. These are the values that we will use in the fit to obtain H. Now, the fit of R/S oc M can be done either using an optimization H procedure for parameterization of the nonlinear regression or a linear regression between \og(R/S) and log(M). The former is theoretically more correct but it was shown in [4] that it provides values of H indistinguishable from the ones obtained by the latter for data sets not too small. For this reason, we will use the linear regression below to determine H. The algorithm that we implement is similar to the one in [18]. It does not include trend correction [4] for simplicity, but this choice should not be seen as limiting since trend correction does not always have a positive effect on the results as we have checked. Hurst did not consider trend correction in his original method [48]. The times series {x(j)} (j = 1, 2 , J V ) analyzed will be such that they have a number of data (N) which is an integral power of two. They are divided into non-overlapping subseries of M data, where M is an integral power of two ranging between 2 and N/2. The number of subseries of length M is N(M) = N / M and is an integer whatever M. We compute the quantity R/S (R/S)M for each where the latter is an average value of M, over the non-overlapping subseries. We now outline the algorithm that we will use. We start by calculating the mean, (x) , of the n-th subseries of length M, which is n M nM j=z(n-l)M+l We also calculate the standard deviation SM UI of the n-th subseries of length M, which Chapter 4. Time Series 115 Analyses where M and not M — 1 is used since the population standard deviation rather than the sample standard deviation is considered. This choice has only little effect on the results as we have checked. Our use of the population standard deviation is guided by the fact that in the literature, it is this quantity that is usually considered. In addition, we compute for all the allowed values of j the quantity j k=(n-l)M+l where (n — 1)M + 1 < j < nM. subtracting the least value of The range R ,M in the n-th subseries is obtained by n from the greatest value of Y M(j) n> Y M(j)nt The rescaled i range of the n-th subseries of length M is determined by the relation {R/S) M = n> A n average rescaled range, denoted by / { R , S ) M = (RTI,M) \NW)) I'(STI.M)- is then defined by (R/S)M, 1 1 \ N(M) £ { R / S ) N ' M i and computed for each allowed value of M. Afterwards, a straight line is fitted to the plot of log[(jR/5)M] versus log(M), the slope of which is precisely the Hurst exponent H. We are now ready to apply this algorithm to the analysis of times series. We will do this first for time series having a known value of H. This preliminary istudy will enable us to test whether the above algorithm gives reliable results. In a second time, we will apply the algorithm to the analysis of the evenly sampled time series of type 2 introduced in section 4.1. Chapter 4.3.4 4. Time Series 116 Analyses T e s t s o f t h e a l g o r i t h m w i t h t i m e series h a v i n g a k n o w n v a l u e o f H The time series that we consider in this subsection were generated by the successive random addition method [107] using the algorithm presented by Feder [31]. As discussed by Churilla et al. [18], it is clear that the time series to be analyzed with the rescaled range algorithm are the time series of the increments, not the time series of the accumulated increments, where the former is obtained from the latter by a first-difference operation. This is because the rescaled range analysis computes the range of cumulative of the data from the mean (R). deviations • We have analyzed with the algorithm of subsection 4.3.3 synthetic time series with a known value of H having N = 2 13 data. The results are plotted on figure 4.21 for H = 0.1, H = 0.5 and H = 0.9. The curves are almost perfect straight lines, except the one for H = 0.1. Also, note that even at large M, the slope of the latter is still larger than 0.1, showing that the slope value 0.1 is reached only asymptotically. The method of generation of the time series is not responsible for the absence of a linear segment ! for the case H = 0.1 on figure 4.21 since we have obtained rather straight lines for the other two cases. It is now well documented (see [4] for instance) that the rescaled range analysis does not work well on time series having negative correlations (i.e. H < 0.5) and our results are simply in agreement with this fact. It is fortunate that the rescaled range analysis works fine when applied to signals having H > 0.5 because, in the real world, the latter are common while the signals with H < 0.5 are relatively rare. Now, we have fitted straight lines to the curves having H > 0.5 on figure 4.21 and found the values H RR « 0.55 and H RR « 0.83 for the curves for H = 0.5 and H = 0.9 respectively. The discrepancies between H RR and H are of about +10% for H = 0.5 and —8% for H = 0.9. As a result, the rescaled range analysis is a biased method which overestimates H when H ~ 0.5, underestimates H when H ~ 0.9 and provides wrong estimates of H Chapter 4. Time Series 117 Analyses when H is small. These findings are in agreement with the ones in [4] iii particular. Note that these limitations of the rescaled range analysis are not so severe if one is interested only in the order of magnitude of H for an unknown signal though. The reliability of the i estimate of H provided by the analysis can be checked for instance by shuffling the order of the data and reapplying the algorithm to the resulting time series. If largely different estimates of H are obtained for the original time series and the shuffled time series, then one can safely conclude that the value of H obtained for the former is representative of a real physical effect. We can now safely apply the algorithm of subsection 4.3.3 to the analysis of the evenly sampled time series of type 2 since it works rather well on synthetic time series having a known value of H, at least if H > 0.5. Our hope is that the time correlations in the evenly sampled time series of type 2 are positive. 4.3.5 R e s u l t s o b t a i n e d w i t h t h e t i m e series o f i n t e r e s t Evenly sampled time series of type 2 for systems with both random o h and constant o h t t are analyzed in this subsection. As mentioned in subsection 4.3.1, a running sum of an evenly sampled time series of type 2 is reminiscent of a fBm time series, even if in fact we 1 have found that it is not a pure fBm time series. The evenly sampled time series of type 2 are therefore suitable to be analyzed with the algorithm of subsection 4.3.3. However, since the phase relationships between the frequency components in these time series are not fully random (due to the nonlinear component), we will investigate if this has an effect on the results of the rescaled range analysis. The results of the analysis for an evenly sampled time series of type 2 for a system with random o th having N = 2 14 data are plotted on figure 4.22. The sampling time interval At is equal to 1.5 times the average time interval between the earthquake sequences. We got similar results with larger values of At but in these cases, the (R/S)M versus Chapter 4. Time Series 118 Analyses 1000- 0 100- 10- ^1 10 I I I I I I I 100 I I I I I I I 1000 I I I I I I 10000 M Figure 4.21: Log-log plot of (R/S)M versus M for three synthetic time series having N = 2 data generated by the successive random addition method. From top to bottom, the curves correspond to the cases H = 0.9, H = 0.5 and H = 0.1. The dashed line at the right of the figure has a slope of 0.1 and is there to illustrate that the curve for H — 0.1 will only slowly reach the slope value H = 0.1 as M increases. Note that each curve was obtained by averaging (R/S)M (for a given M) over 5 realizations corresponding to the same value of H. For each curve, we have excluded the point corresponding to M = 2 since it is trivial and always has (R/S)M=2 — 113 R R Chapter Time 4. Series 119 Analyses M curves cut-off at increasingly smaller M (results not shown) since there are less data. For M = 2, the ratio is always equal to 1, so we did not include this trivial point on figure 4.22. What we observe on the figure is a segment of a given slope for the smallest values of M and another one with a lower slope for the largest values of M. The slope for the segment M € [2 ,2 ] is H ^ 0.90 (see figure 4.22). To check if the decrease of 2 7 the slope at large M is not an artifact of the time series under consideration, we have also analyzed evenly sampled time series of type 2 for 50 x 50 and 25 x 25 systems with random o~ ht The (R/S)M versus M curves for these two time series as well as for the 100 x 100 time series (see figure 4.22) are plotted on figure 4.23, where the sampling time interval At is the same for all three time series. From this figure, we observe that, as the system size gets smaller, then the change in slope tends to occur at smaller M. As a result, the change in slope on figure 4.22 is a real effect, and represents some kind of finite-size effect in time. Note that since all three curves on figure 4.23 span the same M interval, then the change in slope can not be due to the finite number of data in the time series: the system size solely determines at what value of M it will take place. Now, the presence of positive correlations in the 100 x 100 system with random o th is real as we have checked by randomizing the order of the data. The curve obtained by averaging over 5 random shuffles is plotted on figure 4.22 (bottom curve). This curve is rather straight, with a slope equal to H ~ 0.55. Therefore, randomizing the order of the data has caused the destruction of the positive correlations that exist in the system. We have also found similar results for an evenly sampled time series of type 2 for a 100 x 100 system with constant a (see figure 4.24). Again, we notice that the linear segment M G [2 ,2 ] 2 th of the curve has a slope H 7 0.90. There is still a change in slope due to finite-size effects. In addition, we noted that randomizing the order of the data has destroyed the time correlations (see the curve at the bottom of figure 4.24). Thus, we can say that the evenly sampled time series of type 2 for systems with both random errand constant a h t Chapter 4. Time Series 120 Analyses have positive correlations that are affected by finite-size effects. We have also analyzed evenly sampled time series of type 2 for one-dimensional systems and found results (results not shown) which are qualitatively similar to the ones obtained with two-dimensional systems. In addition, time series for which parameter a controlling the degree of stress conservation in the system (see chapter 3) is set to values smaller than 1 provide results again qualitatively similar to the ones obtained with time series for which a = 1. •' In section 4.2, we have conjectured that the evenly sampled time series of type 2 have weakly coherent phase relationships between their frequency components. The latter are responsible for the observation of a weak nonlinear stochastic structure. Now, we would like to investigate if these weakly coherent phase relationships have any effect on the (R/S)M versus M curves. This can be done by obtaining the (R/S)M\ versus M curves for a given time series and for a surrogate time series having the same power spectrum as the original time series but randomized phases (see the algorithm of Theiler et al. [102] for how to generate such surrogate data sets). We have performed such a calculation for a time series for a 100 x 100 system with random o~ h as well as for 5 surrogate time t series of the above type. The results are plotted on figure 4.25 and clearly show that the weakly coherent phase relationships between the frequency components have a minor effect on the (R/S)M versus M curve. Therefore, this confirms that it was correct in the first place to analyze the evenly sampled time series of type 2 with the rescaled range analysis because the nonlinear component is sufficiently weak, so the fractal component of the time series is clearly dominant. Finally, to conclude this section, we give our arguments to justify the results obtained with the evenly sampled time series of type 2. We believe that the existence of strong positive correlations is due to the nature of the driving in the model of X u et al. because the amount of stress added to the system between the earthquake sequences depends Chapter 4. Time Series 121 Analyses 1000 ] - i 1 i i i i i 111 i i i i 1111 i 10 i 100 i i i i 1111 1000 I I i i i 111 10000 M Figure 4.22: Log-log plot of (R/S)M versus M for an evenly sampled time series of type 2 for a 100 x 100 system with random a h having A t equal to 1.5 times the average time interval between the earthquake sequences (N = 2 ). The curve corresponding to this time series is plotted using a solid line and lozenges. The dotted line (H = 0.9) is a fit of the linear segment M G [2 ,2 ] of the latter curve. The curve at;the bottom was obtained with 5 random shuffles of the time series: the curve of (R/S)M versus M was calculated for each realization and an average value of (R/S)M over the 5 realizations was computed for each M. t 14 2 7 Chapter 4. Time Series 122 Analyses 1000- 100- CO 10- I 10 I I I I I I -1—I I I I I I 1000 100 I I I I I I 10000 M Figure 4.23: Log-log plot of (R/S)M versus M for evenly sampled time series of type 2 for L x L (L = 25,50,100) systems with random o h- The sampling time interval A i is equal to 1.5 times the average time interval between the earthquake sequences and the number of data is N = 2 for all three time series. From top to bottom, the curves refer to the 100 x 100, 50 x 50 and 25 x 25 systems respectively. t 14 Chapter 4. Time Series 123 Analyses 1000 10000 Figure 4.24: Log-log plot of (R/S)M versus M for an evenly sampled time series of type 2 for a 100 x 100 system with constant o h having At equal to 1.5 times the average time interval between the earthquake sequences (N = 2 ). The curve corresponding to this time series is plotted using a solid line and lozenges. The dotted line is a fit (H = 0.9) of the linear segment M G [2 ,2 ] of the latter curve. The curve at the bottom was obtained with 5 random shuffles of the time series: the curve of (R/S)M versus M was calculated for each realization and an average value of (R/S)M over the 5 realizations was computed for each M. t 14 2 7 Chapter 4. Time Series 124 Analyses 1000- 100- CO 10 I I I I I — I I I I I I I I 10 100 I 1000 I I I II 10000 M Figure 4.25: Log-log plot of (R/S)M versus M for an evenly sampled time series of type 2 for a 100 x 100 system with random cr . The sampling time interval At is equal to 1.5 times the average time interval between the earthquake sequences and the number of data is N = 2 . The curve plotted using a solid line and lozenges was Obtained with the original time series whereas the one plotted using a dashed line and pluses was obtained with 5 surrogate time series having the same power spectrum as the original time series but randomized phases. For the latter curve, {R/S)M was calculated for each realization, for a given M, and an average value of (R/S)M over the 5 realizations was computed for this value of M. (/l 14 Chapter 4. Time Series 125 Analyses on the past history (see chapter 2). In addition, the correlations are positive (H 0.5) > and not negative (H < 0.5) because the system tends to accumulate stress until a large earthquake sequence occurs and releases a sizable fraction of that stress. This can be seen on figure 4.3 for instance where sometimes the average stress added to the system is much less than the average stress released from it, but this large stress release is due to the past accumulation of stress in the system (stress reservoir is full!). Now, we have also noticed that these positive correlations are affected by finite-size effects. We attempt here an explanation for the latter. We have seen in chapter 3 that the cumulative sizefrequency distribution is affected by the finite size of the system in the sense that the power law for small earthquake sizes does not extend to large earthquake sizes (see figure 3.1 for instance). On the other hand, we have checked in section 4.2 that the histogram (N(5 )) t of the time intervals (5 ) between large earthquakes is roughly exponential, i.e. t proportional to e~ ^°, where 5 is a characteristic time interval and can be determined 5t 0 by a linear regression of log[N(5 )} t versus 6 . We have performed such a linear regression t for the large earthquakes (S > 200) generated in a 100 x 100 system with random o~ h and t found the value So ~ 0.154. We conjecture that this characteristic time interval between large earthquakes is responsible for the decrease of the slope at large M on figure 4.22. Indeed, the cut-off happens at M of size MAt. 128 on this figure, which corresponds to a window Therefore, using the value of At for the corresponding time series, we get MAt — 0.167 ~ 5 . As a result, we have demonstrated that the finite-size effects 0 affecting the (R/S)M large earthquakes. versus M curves are due to a characteristic time interval between These large earthquakes have the effect of resetting the system by destroying the memory effects in it. Since, as we have checked, 8 gets smaller as the 0 system size decreases, then this explains why on figure 4.23 the change in slope occurs at increasingly smaller M as the system size gets smaller. Note in addition that if we would go to arbitrarily large M, then we would see that the slope past the cut-off point Chapter 4. Time Series Analyses 126 approaches 0.5, the value in the absence of any memory effects. 4.4 S u m m a r y of the Important Results and Conclusions In this chapter, we have performed two types of analysis of time series generated by the model of X u et al. These time series are of types 1 and 2. The time series of type 1 involve the size of the earthquake sequences (S) as a function of the occurrence time. They exist in three varieties: varieties (a), (b) and (c) consist in the time intervals between consecutive earthquake sequences with S = So (or S > So), the time intervals between consecutive earthquake sequences whatever S and the sizes of earthquake sequences about evenly sampled in time respectively. The time series of type 2 involve the average stress in the system as a function of the long time. A nonlinear forecasting analysis was first implemented. It investigates whether the data in a scalar time series generated from a dynamical system exhibit low-dimensional chaotic behavior as opposed to high-dimensional (or stochastic) behavior. The analysis was performed by using an algorithm which constructs models for the structure in the time series with a variable smoothing parameter (k). The accuracy of the short-term forecasting error as a function of k reveals much about the underlying dynamics generating the time series. We have applied a nonlinear forecasting algorithm to the analysis of the time series of types 1 and 2. The time series of type 1, whatever the variety, were found to be weakly or not at all forecastable. The only time series of type 1 that have some degree of forecastability are the ones involving the time intervals between large earthquakes generated in one-dimensional systems. For these time series as well as those involving the time intervals between large earthquakes generated in two-dimensional systems, we have obtained histograms of the data and found that they have an exponential distribution. This is in contradiction with the finding that large earthquakes are temporally correlated Chapter 4. Time Series 127 Analyses in the model of X u et al. [16]. We argued that this is because the model of X u et al. in the form used in this thesis is a series model. For all the time series of type 1, no nonlinear structure was found, but they all have a large stochastic component. For this reason, they are of little interest since they seem to be indistinguishable from white noise. The time series of type 2 were found to be forecastable, with a nonlinear but high-dimensional structure. For the evenly sampled time series of type 2, we have checked that the weakly coherent phase relationships between the Fourier components are responsible for the observation of this structure. We suggested that the appearance of nonlinear structure was due to the fact that the system stops following linearly the external drive (threshold dynamics!) when it releases stress through earthquakes. A direct consequence of these results is that self-organized criticality is inherently a high-dimensional phenomenon, as previously suspected by other authors. Note that this conclusion might be particular to the model studied in this thesis. Fractal time series have power-law behavior (/~ ) in their power spectrum. This is /3 however only a sufficient condition since, as shown by Higuchi [46], fractal time series must have random phase relationships between their frequency components. A common example of a fractal time series is provided by the record of the positions of a Brownian particle moving in a one-dimensional Euclidean space. Fractional iBrownian motion (fBm) is a generalization of the ordinary Brownian motion exhibiting long-run time correlations. The latter are positive if H > 0.5, negative if H < 0.5 and zero if H = 0.5 (ordinary Brownian motion), where H is called the Hurst exponent and it characterizes the degree of roughness of the record of the motion. From the results of the nonlinear forecasting analysis, we know that the evenly sampled time series of type 2 are not pure fractal time series since they have weakly coherent (and not truly random) phase relationships between the frequency components, but at least they have an important fractal component. The determination of the Hurst exponent H characterizing the time Chapter 4. Time Series 128 Analyses correlations in time series with an important fractal component can be accomplished by implementing a rescaled range analysis. This analysis consists in studying how the range of cumulative deviations of the data from the mean (R) divided by the standard deviation of the data (S) (or the ratio R/S) varies with the size of the time window (or the number of data M in the window for evenly sampled time series). The' plot of \og{R/S) versus log(M) is linear with a slope H for a pure fractal time series. The quantity we calculated for each M was (R/S)M, which is an average value of R/S over non-overlapping subseries of M consecutive data. The values of M were chosen to be equally spaced on a logarithmic scale. The analysis of the evenly sampled time series of type 2 has led to the I following results: the time series of type 2 have strong positive correlations (H ~ 0.9) which are affected by finite-size effects since H decreases for large time windows. The weakly coherent phase relationships between the frequency components were shown to have little impact on our results. The presence of strong time correlations in the time series of type 2 is due to the nature of the driving in the model because the amount of stress added to the system between the earthquake sequences depends on the past history. The correlations are positive (and not negative) since the system tends to accumulate stress until a large earthquake occurs and releases a sizable fraction of that stress. A characteristic time interval between large earthquakes is responsible for the appearance of finite-size effects in time. The large earthquakes cause the destruction of the memory effects in the system by resetting the latter. Chapter 5 Conclusion The research on complex phenomena has led recently to the introduction of the concept of Self-Organized Criticality (SOC) by Bak, Tang and Wiesenfeld [3]. ;This concept has been presented in order to understand the origin of the power laws observed in nature. Acccording to SOC, certain driven systems consisting of a large number of elements evolve towards a critical state with no characteristic time or length scales. The critical state is reached without fine tuning of a control parameter, unlike for phase transitions in equilibrium statistical physics. A simplified sandpile model illustrates the ideas of SOC. In this model, sand is added grain by grain at a random site on the pile. This random driving process lasts until one site becomes unstable, i.e. when the local slope becomes larger than a threshold value. Then, grains are redistributed to the nearest neighbors of the unstable site and might themselves become unstable as well, and j so on, leading to the creation of an avalanche. After a while, the pile eventually reaches a critical state, characterized by a critical slope, in which additional grains fall off the pile via avalanches of all sizes, distributed in size and lifetime according to a power law. According to Sornette [93], the class of phenomena described by SOC relies on the condition of a very slow driving of the systems and the existence of fast burst-like responses of them. Earthquakes are one of the first example of phenomena that comes to mind for which such a time scale separation is apparent. Our interest in the SOC concept has been oriented towards the understanding of earthquakes in this thesis. The earthquake occurrence is related to the appearance of fractures as well as to 129 Chapter 5. 130 Conclusion friction. In principle, both fracture and friction should be incorporated in an earthquake model. However, in practice, this is a formidable task and realistic situations have yet to be correctly modeled. The model that we have used throughout this thesis studies the creation of faults in a planar elastic region and was introduced by Xujet al. [113, 114]. This model has its validity at much larger length scales than the atomistic one, where the medium can be described by continuous fields. The model of X u et al. involves the following: • it focuses on a one-layer planar region which is discretized into L \ x L squares, 2 • it assigns the stress tensor directly to each square, ! i • it solves the elastic equations for all the squares before and after the rupture of a given square. In the model, a square ruptures when a > a h, where a is the stress on that square t (a = a xy if there is only one shear mode of rupture, a xy being one component of the stress tensor) and o~ h the stress threshold, and this is followed by a double couple stress t redistribution. We have applied the model of X u et al. to the study of earthquakes by implementing the following: • embed the region of interest in an infinite medium and neglect'the effect of any activity originating outside the region (i.e. we use open boundary conditions), • uniformly drive the region by adding stress amount ACT to all the squares (this is enough to cause the rupture of the square which has a the closest to a ), th the surrounding medium provides the necessary external drive, • increment the long time scale (defined by variable t) by the amount ACT, where Chapter 5. 131 Conclusion • once an earthquake sequence has begun, freeze t since earthquakes are instantaneous events on the long time scale, • reset o th of a ruptured square to a new value (for a system with random cr ) or set tfl it to a constant value (for a system with constant a h), t • when more than one square break simultaneously, use the independent redistribution procedure, • count the total number of broken squares (S) during the sequence to estimate the size of the earthquake. The parameters of the model are: • 9i, which controls the width of the uniform (random) distribution of stress thresholds (a th e [0*,i.o]), • X, which fixes the fraction of stress retained by a ruptured square, • a, which introduces stress nonconservation in the redistribution of stress (a = 1 corresponds to the conservative case). Once the model has reached the stationary state, we have checked that the GutenbergRichter power law C(S) oc S~ is satisfied, where C(S) is the cumulative fraction of T earthquakes with S or more broken squares. However, this power law is affected by finite-size effects: there are the effects of the finite size of the system and of the finite size of the squares. To take these effects into account, we have performed a.finite-size-scaling analysis. For simplicity, we have considered systems with L x 1 squares (i.e. lattices arranged in a strip). The finite-size-scaling function that we have considered includes a factor which corrects for the finite-size effects affecting the power law:; C(S, L) = S- g{(S r - l)L- /f(X, v 0 U a)}. Chapter 5. 132 Conclusion The results that we have obtained are, for a = 1: • for fixed parameter values, when L is allowed to vary, then the finite-size-scaling function works well over the whole range of S for systems with random a h and for t S > 50 for systems with constant a h- For systems with constant a h, there is a t t crossover region for S < 50, • for fixed L, when 8i or X is allowed to vary, then the above function works fine for all S for systems with random a h and for S > 50 for systems with constant o~ ht t These results imply that a system which has been discretized by a fine mesh and a small X (or a large 9{) can be approximated by a coarser mesh using a larger X (or a smaller #/). In all the cases, r ?s 0.30 gives the best results and is a robust value for L x 1 systems. However, for L x L systems, a preliminary check has convinced us that r ~ 0.4 provides optimal results. This shows that, at least quantitatively, the results obtained with L x 1 systems are particular to this lattice configuration. We have also considered a mean-field version of the model of X u et al. in which all the elements are equally coupled with each other and discovered that the finite-size-scaling function with r ^ 0.30 does not work well. This shows that r is also sensitive to the details of the model. Despite the sensitivity of r to the lattice configuration and the details of the model, we can state that the model of X u et al. displays SOC since the finite-size-scaling hypothesis works fine. This conclusion holds only in the conservative case (a = 1). In the non-conservative case (a < 1), the model is not self-organized critical since, as we have checked, a characteristic earthquake size So < L appears, So being a function of a and L. The quantity that we have considered in the finite-size-scaling analysis is the power law C(S) oc S~ , which is an average over time as well as over space of the activity in T Chapter 5. 133 Conclusion the system. In the stationary state, this is a static distribution. We have also studied in this thesis the dynamic quantities of the model, such as, for instance, the size of the earthquake sequences (S) and the average stress in the system (a), o referring to a spatial average over the squares of the system. The sets of values of S and a over the long time (t) make up time series which we have called time series of types 1 and 2 respectively. The time series of type 1 exist in three varieties: varieties (a), (b) and (c) consist in the time intervals between consecutive earthquake sequences with S = So (or S > So), the time intervals between consecutive earthquake sequences whatever S and the sizes of earthquake sequences about evenly sampled in time respectively. These time series have been analyzed using two different types of analysis. ' A nonlinear forecasting analysis was first implemented. It investigates whether the data in a scalar time series generated from a dynamical system exhibit low-dimensional chaotic behavior as opposed to high-dimensional (or stochastic) behavior. That analysis was performed by using an algorithm which constructs models for the structure in the time series with a variable smoothing parameter (k). The accuracy of the short-term forecasting error as a function of k reveals much about the underlying dynamics generating the time series. We have applied a nonlinear forecasting algorithm to the analysis of the time series of both types. The time series of type 1, whatever the variety, were i found to be weakly or not at all forecastable. The only time series of type 1 that have some degree of forecastability are the ones involving the time intervals between large earthquakes generated in L x 1 systems. For all the time series of type 1, no nonlinear structure was found, but they have a large stochastic component. For this reason, they are of little interest since they seem to be indistinguishable from white noise. On the other hand, the time series of type 2 were found to be forecastable, wifh a nonlinear but i high-dimensional structure. We suggested that the appearance of nonlinear structure in Chapter 5. 134 Conclusion these time series was due to the fact that the system stops following linearly the external drive (threshold dynamics!) when it releases stress through earthquakes. A direct consequence of these results is that S O C is inherently high-dimensional, at least in the model studied in this thesis. It might be possible to observe a low-dimensional behavior in other versions of the model, but this remains to be proven. Since it is usually the magnitude of the earthquakes which is known in practical situations, the results obtained with the time series of type 1 seem to close all doors towards the prediction of earthquakes. However, it is relevant to note that the definition of S used in this thesis is not tied to the state variable (i.e. the stress) in the model. If we would have considered as a definition of S the difference between the average stress in the system before and after the earthquakes (i.e. a^f — a / ), then the results obtained Q t with the time series of type 2 suggest that it might be possible to make some predictions on the latter quantity. On a related note, we mention that o\, t — a / e a t is of the same order of magnitude as the quantity called stress drop in the literature on seismology, the latter being related to the magnitude of the earthquakes (see [51] for instance). The second type of time series analysis that we have performed is a rescaled range analysis. In principle, it is applicable only to fractal time series, which have random phase relationships between their frequency components as well as a power-law behavior in their power spectrum [46]. However, in practice, it is usually sufficient for the time series to have an important fractal component for the rescaled range analysis to apply. This is the case for the evenly sampled times series of type 2 which, in addition toj a weak nonlinear component, have an important fractal component as we have checked. The rescaled range analysis determines an exponent, the Hurst exponent H (0 < H < 1), characterizing the time correlations in time series with an important fractal component (and in general i in fractal time series). This analysis consists in studying how the range of cumulative deviations of the data from the mean (R) divided by the standard deviation of the data Chapter 5. 135 Conclusion (S) (or the ratio R/S) varies with the size of the time window (or the number of data M in the window for evenly sampled time series). The plot of \og(R/S) versus log(M) is a straight line with a slope H for fractal time series. We have used an algorithm to achieve that task and found the following results for the evenly sampled time series of type 2: they have strong positive correlations (H « 0.9 > 0.5) which are affected by finite-size effects since H decreases for large time windows. The presence of strong time correlations in the time series of type 2 is due to the nature of the driving in the model because the amount of stress added to the system between the earthquake sequences depends on the i past history. The correlations are positive and not negative (when H < 0.5) since the system tends to accumulate stress until a large earthquake occurs and releases a sizable fraction of that stress. We have checked that a characteristic time interval between large earthquakes causes the appearance of the finite-size effects in time. The large earthquakes are responsible for the destruction of the memory effects in the system because they reset the latter. - The work performed in this thesis has enabled us to better understand the earthquake phenomenon but is in no way exhaustive. In the future, we suggest considering more realistic versions of the model of X u et al., for instance a parallel version of it. A parallel model treats the situations where more than one square rupture simultaneously in a selfconsistent manner, by assuming that the rupture at one site is affected by the ruptures at the other sites. Chen et al. [16] have initiated that task, but much work remains to be done though. For instance, it would be interesting to see if improyed forecastability is obtained with the time series of type 2 and if the memory effects are modified in this version of the model. A n extra suggestion for the future is to focus on a definition of the size of the earthquakes (S) which is more closely tied to the dynamics than the one used in this thesis: S = d j — d f . In particular, it would be interesting to investigate, by a be a t finite-size-scaling analysis, what is the effect on the power-law exponent r of this choice Chapter 5. Conclusion 136 for S. This should be done for both the series and the parallel versions of the model. Bibliography i [1] P. Bak and K . Sneppen, Phys. Rev. Lett. 71, 4083 (1993). ; i [2] P. Bak and C. Tang, J. Geophys. Res. B 94, 15,635 (1989). j [3] P. Bak, C. Tang and K . Wiesenfeld, Phys. Rev. Lett. 59, 381 (1987). [4] J. M . Bassingthwaighte and G. R. Raymond, Ann. J. Biomed. Eng. 22, 432 (1994). [5] J. M . Bassingthwaighte and G . R. Raymond, Ann. J. Biomed. Eng. 23, 491 (1995). [6] R. Bhagavatula, K . Chen and C. Jayaprakash, J. Phys. A:Math. Gen. 27, L155 (1994). ! [7] S. Bottani, Phys. Rev. Lett. 74, 4189 (1995). [8] G. E . P. Box and G. M . Jenkins, Time Series (Holden-Day,San Francisco, 1970). j Analysis Forecasting and Control 1 [9] W . A . Brock, W . D. Dechert and J. Scheinkman, A Test for Independence Based on the Correlation Dimension [Working Paper, Department of Economics, University of Wisconsin (Madison), University of Houston, University of Chicago, 1986]. [10] J. Buck, Quart. Rev. Biol. 63, 265 (1988). [11] R. Burridge and L . Knopoff, Bull. Seism. Soc. A m . 54, 1875 (1964). [12] R. Burridge and L . Knopoff, Bull. Seism. Soc. A m . 57, 341 (1967). [13] J. M . Carlson and J. S. Langer, Phys. Rev. Lett. 62, 2632 (1989): [14] M . Casdagli, J. Roy. Statist. Soc. B 54, 303 (1991). [15] K . Chen, P. Bak and S. P. Obukhov, Phys. Rev. A 43, 625 (1991). [16] K . Chen, R. Bhagavatula and C. Jayaprakash, J. Phys. A:Math. Gen. 30, 2297 (1997). j [17] K . Christensen and Z. Olami, Phys. Rev. A 46, 1829 (1992). ; [18] A . M . Churilla, W . A . Gottschalke, L . S. Liebovitch, L. Y . Selector, A . T. Todorov and S. Yeandle, Ann. J. Biomed. Eng. 24, 99 (1996). | 137 138 Bibliography [19] A . Corral, C. J. Perez, A . Diaz-Guilera and A . Arenas, Phys. Rev. Lett. 75, 3697 (1995); A . Corral, C. J. Perez, A . Diaz-Guilera and A . Arenas, Phys. Rev. Lett. 74, 118 (1995). •' [20] P. Cowie, C. Vanneste and D. Sornette, J. Geophys. Res. B 98, 21,809 (1993). [21] J. P. Crutchfield, J. D. Farmer, N . H . Packard and R. S. Shaw, Sci. A m . 255 (No 6), 46 (1986). [22] J. P. Crutchfield and K . Kaneko in Directions (World Scientific,Singapore,1988). in Chaos, volume: I, ed H . Bai-Lin [23] J. P. Crutchfield and B . S. McNamara, Complex Syst. 1, 417 (1987). [24] A . Di'az-Guilera, Europhys. Lett. 26, 177 (1994). [ [25] F. Donze, P. Mora and S.-A. Magnier, Geophys. J. Int. 116, 46 (1994). [26] B . Drossel and F. Schwabl, Phys. Rev. Lett. 69, 1629 (1992). [27] P. England and D. McKenzie, Geophys. J. Roy. Astr. Soc. 70, 295 (1982). [28] E. Fama, Management Science 11, 404 (1965). [29] J. D. Farmer and J. J. Sidorowich, Phys. Rev. Lett. 59, 845 (1987). i [30] J . D . Farmer and J. J. Sidorowich in Evolution, Learning and Cognition, Lee (World Scientific,Singapore,1988). [31] J. Feder, Fractals ed Y . C. (Plenum Press,New York, 1988). [32] S. Feng and P. N . Sen, Phys. Rev. Lett. 52, 216 (1984). [33] A . M . Fraser, Physica D 34, 391 (1989). [34] A . M . Fraser and H . L . Swinney, Phys. Rev. A 33, 1134 (1986). [35] A . Gabrielov, W . I. Newman and L . Knopoff, Phys. Rev. E 50, 188 (1994). [36] L . G i l and D. Sornette, Phys. Rev. Lett. 76, 3991 (1996). [37] B . V . Gnedenko and A . N . Kolmogorov, Limit Distributions Random Variables (Addison-Wesley,Reading, 1954). for Sum of [38] P. Grassberger, Phys. Rev. E 49, 2436 (1994). [39] P. Grassberger and I. Procaccia, Phys. Rev. Lett. 50, 346 (1983). : Independent 139 Bibliography [40] H . S. Greenside, A . Wolf, J. Swift and T. Pignataro, Phys. Rev. A 25, 3453 (1982). [41] G . Grinstein, D.-H. Lee and S. Sachdev, Phys. Rev. Lett. 64, 1927 (1990). [42] D. Groleau, B . Bergersen and H.-J. X u , J. Phys A : Math. Gen. 30, 3407 (1997). [43] B . Gutenberg and C. F. Richter, Bull. Seism. Soc. A m . 34, 185 (1944). [44] M . Henon, Commun. Math. Phys. 50, 69 (1976). [45] H . J. Herrmann and S. Roux (eds) Statistical Media Models for the Fracture (North-Holland,Amsterdam,1990). of Disordered ' [46] T. Higuchi, Physica D 46, 254 (1990). [47] N . F . Hunter Jr. in Nonlinear Modeling and Forecasting, Eubank (Addison-Wesley,Redwood City,1991). eds M . Casdagli and S. [48] H . E. Hurst, Trans. Amer. Soc. Civ. Engrs. 116, 770 (1951). [49] T. Hwa and M . Kardar, Phys. Rev. Lett. 62, 1813 (1989). [50] I. M . Janosi and J. Kertesz, Physica A 200, 179 (1993). [51] H . Kanamori and D. L . Anderson, Bull. Seism. Soc. A m . 65, 1073 (1975). [52] Y . Kantor and I. Webman, Phys. Rev. Lett. 52, 1891 (1984). [53] J. D. Keeler and J. D. Farmer, Physica D 23, 413 (1986). i [54] M . Kendall and J. K . Ord, Time Series (Edward Arnold,Kent,1990). [55] L . Knopoff in Disorder and Fracture, (Plenum Press, New York, 1990). eds J. C. Charmet, S. Roux and E. Guyon ' [56] L. Lam and V . Naroditsky (eds) in Modeling Complex Phenomena, Proceedings of the Third Woodward Conference, San Jose State University, April 12-13 1991 (Springer-Verlag,New York,1992). [57] P. Levy, Theorie de VAddition 1954). des Variables Aleatoires (Gauthier-Villars,Paris, 1937- [58] L . Lise and H . Jedtoft Jensen, Phys. Rev. Lett. 76, 2326 (1996). [59] A . E . H . Love, A Treatise on the Mathematical York,1944). Theory of Elasticity j (Dover,New 140 Bibliography [60] B . Malraison, P. Atten, P. Berge and M . Dubois, J. Phys. Lettre (France) 44, L897 (1983). B. B . Mandelbrot, The Fractal Geometry of Nature (Freeman,San Francisco, 1983). B. B . Mandelbrot and J . W . Van Ness, SIAM Rev. 10, 422 (1968). B. B . Mandelbrot and J . R. Wallis, Water Resours. Res. 4, 909 (1968). B. B . Mandelbrot and J . R. Wallis, Water Resours. Res. 5, 228-241; 242-259; 260267; 321-340; 967-988 (1969). R. M . May, Nature 261, 459 (1976). P. Meakin, G . L i , L. M . Sander, E . Louis and F . Guinea, J . Phys. A:Math. Gen. 22, 1393 (1989). A. A . Middleton and C. Tang, Phys. Rev. Lett. 74, 742 (1995). 1 P. Miltenberger, D. Sornette and C. Vanneste, Phys. Rev. Lett. 71, 3604 (1993). R. E . Mirollo and S. H . Strogatz, S I A M J. Appl. Math. 50, 1645 (1990). E. W . Montroll and M . Shlesinger in Nonequilibrium Phenomena II, From Stochastic to Hydrodynamics, Studies in statistical mechanics X I , eds J . L/Lebowitz and E . W. Montroll (North-Holland,Amsterdam,1984). P. Mora and D. Place, Int. J . Mod. Phys. C 4, 1059 (1993). K . Nagel and E. Raschke, Physica A 182, 519 (1992). H. Nakanishi, Phys. Rev. A 43, 6613 (1991). W . I. Newman and A . M . Gabrielov, Int. J . Fract. 50, 1 (1991). ; C. P. North and D. I. Halliwell, Math. Geol. 26, 531 (1994). S. Obukhov in Random Fluctuations and Pattern Growth, eds H . E . Stanley and N . Ostrowsky (Kluwer,Dordrecht,1988). Z. Olami, H . J. S. Feder and K . Christensen, Phys. Rev. Lett. 68, 177 (1991). F. Omiri, J . College Sci. Tokyo Imper. Univ. 7, 111 (1895). A . R. Osborne and A . Provenzale, Physica D 35, 357 (1989). N . H . Packard, J . P. Crutchfield, J . D . Farmer and R. S. Shaw, Phys. Rev. Lett. 45, 712 (1980). 141 Bibliography [81] V . Pareto, Cours d'economie politique. Reprinted as a volume of Oeuvres (Droz,Geneva,1896-1965). [82] C. S. Peskin in Mathematical Aspects of Heart Physiology Completes (Courant Institute of Mathematical Sciences,New York University,New York,1975). [83] L . Pietronero, A . Vespignani and S. Zapperi, Phys. Rev. Lett. 72, 1690 (1994). ; [84] M . Plischke and B . Bergersen, tific,Singapore,1994). Equilibrium Statistical Physics (World Scien- [85] W . H . Press, S. A . Teukolski, W . T. Vetterling and B . P. Flannery, Numerical in Fortran (Cambridge University Press,Cambridge,1992). Recipes [86] S. Roux and E. Guyon, J . Phys. Lettre (France) 46, L999 (1985). [87] T. Sauer, J . A . Yorke and M . Casdagli, Embedology Maryland,1990). [88] D . Saupe in The Science of Fractal Images, (Technical Report,University of eds H . O. Peitgen and D. Saupe (Springer-Verlag,New York,1988). [89] C. H . Scholz, The Mechanics Press,Cambridge, 1990). of Earthquakes and Faulting (Cambridge University [90] A . Sherman, J . Rinzel and J . Keizer, Biophys. J. 54, 411 (1988). [91] R. F . Smalley, D. L . Turcotte and S. Solla, J. Geophys. Res. B 9Q, 1894 (1985). [92] D. Sornette, J. Phys. I (France) 2, 2065 (1992). [93] D. Sornette, J. Phys. I (France) 4, 209 (1994). [94] D. Sornette, A . Johansen and I. Dornic, J . Phys. I (France) 5, 325 (1995). [95] D . Sornette, P. Miltenberger and C. Vanneste, Pageoph. 142, 491 (1994). [96] A . Sornette and D. Sornette, Europhys. Lett. 9, 197 (1989). [97] J . A . Steketee, Can. J . Phys. 36, 192 (1958); J. A . Steketee, Can. J . Phys. 36, 1168 (1958). [98] K . Stokbro and D . K . Umberger in Nonlinear Modeling and Forecasting, Casdagli and S. Eubank (Addison-Wesley,Redwood City, 1991). [99] G . Sugihara and R. M . May, Nature 344, 734 (1990). eds M . Bibliography 142 [100] H . Takayasu and H . Inaoka, Phys. Rev. Lett. 68, 966 (1992). [101] F. Takens in Dynamical Systems and Turbulence, eds D. Rand and L.-S. Young (Springer-Verlag,Berlin,1981). I [102] J. Theiler, B. Galdrikian, A . Longtin, S. Eubank and J. D.[ Farmer in Nonlinear Modeling and Forecasting, eds M . Casdagli and S. Eubank (AddisonWesley,Redwood City,1991). [103] H . Tong, Non-linear Time Series: A Dynamical System Approach (Clarendon Press,Oxford,1990). . [104] H . Tong and K . Lim, J. Roy. Statist. Soc. B 42, 245 (1980). [105] P. Tong and J. Rossetos, Finite Element Method: Basic Technique and Implemen- tation (MIT Press,Cambridge,1977). [106] J. P. Vilotte, M . Daignieres and R. Madariaga, J. Geophys. Res. B 87, 10,709 (1982). [107] R. F . Voss in The Science of Fractal Images, eds H . O. Peitgen and D. Saupe (Springer-Verlag,New York, 1988). ; [108] J. R. Wallis and N . C. Matalas, Water Resours. Res. 6, 1583 (1971). [109] L . M . Ward in press in Nonlinear Dynamics in Human Behavior, ed W . Sulis (World Scientific,Singapore). [110] W . Weibull, J. Appl. Mech. 8, 293 (1951). [Ill] B . J. West and A . L. Goldberger, A m . Scientist 75, 354 (1987). [112] A . Wolf, J. B . Swift, H . L . Swinney and J. A . Vastano, Physica D 16, 285 (1985). [113] H.-J. X u , B. Bergersen and K . Chen, J. Phys. A:Math. Gen. 25,iL1251 (1992). [114] H.-J. X u , B. Bergersen and K . Chen, J. Phys. I (France) 3, 2029 (1993). ! [115] Y . Yamamoto and R. L. Hughson, Physica D 68, 250 (1993). [116] G . U . Yule, Philos. Trans. Roy. Soc. London A 226, 267 (1927). ; Appendix A Lattice G r e e n Function Simplification and Properties Firstly, we want to illustrate how to simplify the calculation of the lattice Green function (G) introduced in the discretization scheme defined by X u et al. and presented in chapter 2. This will be done by reducing the number of integrals in (2.6) fromftwo to one. This i will allow a much more accurate evaluation of the value of G(f — f ), especially at large 0 r-f . ; 0 Secondly, we want to give some useful properties of the lattice Green function which are easy to prove. A.l Simplification of the Lattice G r e e n F u n c t i o n Referring to (2.6) and (2.8), the lattice Green function is expressed as G , _ f f o ) = [* dk^r 2TT J-IK dk sin k sin k ^. _ ) 2TT (1 — cos k cos k ) 2 y 2 x y (f fo ( A ^ 2 J-TT x y where f — f = (X, Y), X and Y being two integer numbers proportional to the lattice 0 constant a (defined as 1 in this work for simplicity) and which can be negative, zero or positive. If the lattice has L\ x L squares, then | X |= 0,1, 2 , L 2 | Y |= 0,1,2,...,Z, - 1. We will see below that G(X,Y) 2 x — | l and = G(-X,Y) = G{X,-Y) = G(—X, — Y) so it is sufficient to focus on X > 0 and Y > 0 throughout this appendix. Then, k • (r — r ) = k X + k Y and the exponential factor in ( A . l ) can be written as 0 x y iH?-r ) e 0 = ( C0S kxX + k Y) y 143 + isin(k X x + kyY). (A.2) Appendix A. Lattice Green Function Simplification and Properties 144 The right-hand side of (A.2) can be further developed using obvious trigonometric identities to give four terms. Out of these four terms, only one provides a non-zero contribution to G(f — f ) whatever the values of X and Y [cos(k X) cos(k Y)}. Q x The other three terms y are such that they render at least one of the integrals in ( A . l ) odd and' as a result, their contributions are zero since the integration is done over a complete cycle. Keeping only the relevant term of (A.2), ( A . l ) becomes 7-7T 2-7T J-n 2TT (1 — cos k cos k Y T (A.3) v Let us rewrite (A.3) as sin ky cos(kyY) dk (11 — c x cos ky) r7r G(X,Y) = - ^sin £; cos(£; X) x 'KJ-TT 2 / 2 x ZTT y JO 2 (A.4) where c = cos k is treated as a constant when the integration over k is done. Defining x y the integral over k as / , where y r n sin k cos(k Y) (1 — c x cos k ) ' 2 y Jo d k y y (A.5) 2 y (A.4) can now be written as 1 r n G(X, Y) = - rlk — sin k cos(k X) x (A.6) x /. 2 x The next step is to obtain an analytical expression for I for various: values of Y. The simplest case is Y = 0. Setting Y = 0 in (A.5), we have: 7T IY=O = / Jo r ' • ' sin L, 2 1 — c x cos k yJ Using a software such as Mathematica, an analytical expression for I =o can be easily Y derived: (l-VT^) IY=> Y=0 145 Appendix A. Lattice Green Function Simplification and Properties For values of Y different from 0, I can be simplified further by using software Mathematica. We thus obtain, starting from (A.5), 1 r* = 2c4 -2 cos{k Y) + cY cos [(1 - Y)k ] - cY cos [(1 + Y)k ] y y y — 1 + c x cos k d k y (A.8) v It is not possible to simplify I further without fixing Y > 0. For simplicity, we consider only the smallest possible values of Y > 0. Then, setting Y = 1 in (A.8), and using software Mathematica, we obtain, after some simplification of the resulting expression, r 7T iy=i = -2 + 2-c 2 (A-9) Repeating this procedure for Y = 2, Y = 3 and Y = 4, we obtain respectively IY=3 — 7T iV"=4 — 6 - 5c ' 6 IT 2- - + IY=2 c 2 c \/l- c 2 2 (A.10) 2 , . 16 - 18c + 3c 2 ( - 8 + 5c ) + 2 IT 4 22 (A.11) vT 40 - 56c + 17c 2 4(-10 + 9c - c ) + 2 4 4 (A.12) vT We could pursue this procedure for all the necessary values of Y but we stop here for brevity. Now, we can substitute back into (A.6) the analytical expression of I for the given value of Y and try to cast the expression in a convenient form. Let us illustrate this in the case Y = 0. Using c = cosk in (A.7) and substituting in (A.6), we have x 7T 1 — \J\ — cos k dk. —?-sm k cos(k X) 2 fjL~ x 2 / x x cos k y/\ — cos k 2 LIT -7r 2 x x Using obvious trigonometric identities, we can simplify the integrand; in the latter expression and cast the integral in the form r dk sin k cos(k X)(l— | sinA^ |) w 7T / -7T 2 x 2TT x x cos k I sin k 2 x x (A.13) 146 Appendix A. Lattice Green Function Simplification and Properties Now, k x the integrand of (A.13) is even in k and so, using the fact that sin£; > 0 for x x £ [0,7r], we can get rid of the absolute values and simplify (A.13) further: r dk — Jo f ^,-,r ^ G(X,0)= n sin k cos(k X) ^ - 1 , 1 + sin k x x , „ (A.14) x A x We can make the change of variables k —> k '+n/2 which leads to sin(k '+ir/2) = cos k and cos[(A; ' + ir/2)X] x x x x x cos(fX) - sin(k 'X) = cos(k 'X) x x sin(f X ) . Out;of the two terms in the latter identity, only the first one leads to an integrand even in M ' (the bounds of x are — TT/2 and 7r/2) and so, to a non-zero contribution to G(X, 0). the integral over k x As a result, (A.14) becomes rv- ^; °t-?. c 0 g(*,o) = ± ( g ) c \7r/ 7o 1 + cosk 1 (A.i5) x where we have used the fact that for an integer X, c o s ( | X ) = 1 for X = 0,4,8,... [(+) sign in (A.15)], cos(f X) = - 1 for X = 2, 6,10,... [(-) sign in (A.15)], and cos(fX) = 0 ! for X odd. We can repeat a similar procedure for G(X, 1), G(X,2), G(X,3) and G(X,4) using (A.9), (A.10), (A.11) and (A.12) respectively. For brevity, we just give the results: (-) G(X,i)=± ' G(X,2) r ' d Jo [2\ = ± ( - ) W2 & - ™ * * ' ) M * * ' 2 k x K G(X 3) = ± f-1 f eta /( /-./v ^ G(x,4) = 1 \TV) Jo r ( ±{-)j Q / 2 x x x x 3cos O( 1 CQS (1 + cos k ') x x ^ ) ( ^) , sin fc tan A; ' 2 x x j/ cos(k 'X) > V Y (1 + COS k \7T/ 70 ' • x (i - 2cosk ')cosk dk ' - /2 ) x (1 + cos kj) tan k x {l-4cosk ')(l-cosk ')cosk 'cos{k 'X) l dk x x ( 1 + x c o s i f c x 0 3 x x ; • When y is an odd integer, then the (+) sign applies to X = 1,5,9,... and the (—) sign to X = 3,7,11,... whereas when Y is an even integer, then the (j-) sign applies to X = 0, 4, 8,... and the (-) sign to X = 2, 6,10,.... : Appendix A. Lattice Green Function Simplification and Properties 147 To summarize, let us give the expressions for G(X, Y) when Y is fixed to a given value (odd or even integer) different from 0 (the case Y = 0 is special and then G(X, 0) is given by (A.15)). It can be easily checked that the above results for G(X, 1), G(X, 2), G(X, 3) and G(X, 4) are special cases of the following general expressions: • Y is an even integer different from 0 2^ /-"/a ,(1 - y c o s f c ' ) ( l - cosA: ')T- cos'A; 'cos(A; 'A:) 1 s I x (l + cosk ')$ :r ' +1 x (A.16) where the (+) sign is for X = 0,4, 8,... and the (—) sign for X = 2, 6,10,... • Y is an odd integer G(X,Y) = ± C-) Vvr/7o r^,/l-yco (1 + cos k ) x s f c ,-)(l-co^,-)^sinj(^) 2 tan£; ' x where the (+) sign is for X = 1, 5, 9,... and the (—) sign for X = 3, 7,11,... These expressions clearly show that if X + Y is an even integer, then G(X,Y) ^ 0 whereas if X + Y is an odd integer, then G(X, Y) = 0. We will see in appendix B how to overcome this unphysical property of G(X,Y). For now, let us describe some useful properties of G(X, Y). A.2 Properties of the Lattice G r e e n F u n c t i o n Using (A.3), it can be checked that the following property applies to G(X,Y): G(X,Y)=G(-X,Y). The proof of such a property can be done by replacing in the expression for G(—X, Y) k by — k and reorganizing the integrand, which will turn out to be identical to the one x x of G(X, Y). In a similar way, the following properties also apply: I G(X,Y) = G(X,-Y) Appendix A. Lattice Green Function Simplification G(X,Y) = and 148 Properties G(-X,-Y). Another useful property of G(X, Y) is G(X,Y) = G(Y,X) which can be proved by replacing in the expression for G(Y, X) [using (A.3)] k by k and x y k by k . These changes of variables can be done since k and k are dummy variables. y x x y That property greatly simplifies the task of calculating G(X, Y) for all X and Y. Indeed, we can then limit ourselves to the computation of G(X, Y) for X > Y. We have also checked that G(X,Y) falls off as 1/r 2 at large r = \/X 2 + Y, 2 where r is the distance from the rupture location to a given square of the lattice. The latter property can be useful when dealing with very large lattice sizes because the calculation of G(X,Y) using (A.16) and (A.17) is sometimes not accurate at large r. Finally, it is of interest to note that G(X, Y) can be multipled by a global constant without changing the redistributed stresses. This can be seen by referring to (2.8). Appendix B C o a r s e - G r a i n i n g of the Lattice G r e e n Function As we have seen in appendix A , G(X, Y) — 0 for X + Y equal to an odd integer and G(X,Y) ^ 0 for X + Y equal to an even integer, where X and Y are assumed to be positive or zero but not negative. (The latter restriction is not severe since from the properties of G(X, Y) outlined in appendix A, the signs of X and Y are irrelevant). This unphysical property of G(X, Y) can be overcome by a coarse-graining procedure. First, we illustrate this procedure on a one-dimensional strip having 2L\ squares. Such a lattice configuration is considered in chapter 3 for the study of the scaling properties of the model of X u et al. We imagine that one of the elementary squares (called the source square) breaks. During the stress redistribution which follows, we can mentally combine the squares pairwise, so that there are now L\ pairs of squares in the lattice. The lattice Green function associated with the pairs can be expressed as G\X,0) = G(2X,0), where 2X is the horizontal distance between the left (right) source square to a given left (right) square (in unit of a), and X = 0,1, 2 , L \ — 1 is the distance between the pairs in unit of 2a (a was fixed to 1 in chapter 2). For completeness, it can be mentioned that in the case of a two-dimensional (2Iq) x (2L ) lattice of squares having a = 1, we can, during the stress redistribution 2 which follows the rupture of an elementary square (again called the source square), generalize the above procedure by mentally combining four of them such as to form a bigger 149 Appendix B. Coarse-Graining of the Lattice Green Function . 150 square with lattice constant a = 2. The lattice now has L\ x L squares with a = 2. The 2 most natural coarse-grained lattice Green function can be expressed in this case as G (X,Y) = 2 G(2X,2Y) G(2X + ^[G(2X+ - 1,2Y + 1,2Y+ 1) + G(2X+ 1,2Y-1) 1) + G{2X - 1,2Y - 1)] where X = 0,1, 2 , 1 ^ — 1 and Y = 0,1, 2 , L 2 + \ — 1 refer to distances between the squares with a = 2. i G and G will be used in (2.8) to replace G when the lattice of squares is onej dimensional and two-dimensional respectively. 1 2 Now, using the properties of the elementary (not coarse-grained) lattice Green funci tion (G) outlined in appendix A , we can easily prove the following [properties of the coarse-grained lattice Green function G (it refers to both G and G ): 1 G(X, Y) = G(-X, Y) = G(X, G(X,Y) We have also checked that G(X,Y) = -Y) = 2 G(-X,-Y) G(Y,X). falls off as 1/r 2 at large r = s/X 2 + Y. 2 It is remarkable that the elementary lattice Green function (G) and the coarse-grained lattice Green function (G) have similar properties. This gives us confidence that the coarsegraining procedure that we use does not alter G that much.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Study of the scaling and temporal properties of a simplified...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Study of the scaling and temporal properties of a simplified earthquake model Groleau, Daniel 1997
pdf
Page Metadata
Item Metadata
Title | Study of the scaling and temporal properties of a simplified earthquake model |
Creator |
Groleau, Daniel |
Date Issued | 1997 |
Description | Certain driven systems consisting of a large number of elements evolve towards a critical state with no characteristic time or length scales. This class of phenomena is described as Self-Organized Criticality (SOC). SOC relies on the condition of a slow driving of the systems and the existence of fast burst-like responses of them and includes earthquakes. We employ a model proposed by Xu et al. for ruptures in an elastic medium, subject to shear stress, and apply it to the study of earthquakes. In the model, the size of an earthquake is defined as the number of ruptures occurring sequentially on the basic units of discretization (squares) of the medium. A histogram of the earthquake sizes shows that the model is not completely scale invariant due to finite-size effects. To take them into account, we implement a finite-size-scaling analysis. The results of this analysis show that the model is scale invariant only when there is stress conservation. So, the model displays SOC in the conservative case only. We also study the dynamic quantities of the model, in particular the average stress in the system. The sets of average stress values are analyzed using two types of time series analysis. The nonlinear forecasting analysis investigates whether time series exhibit low-dimensional chaotic behavior as opposed to high-dimensional (or stochastic) behavior. We find that the above time series have a nonlinear structure, but with a substantial stochastic component, so SOC is inherently high-dimensional. The appearance of nonlinear structure is due to the fact that the system stops following linearly the external drive when it releases stress through earthquakes. The rescaled range analysis characterizes the time correlations (or memory effects) in the time series. We find strong positive time correlations in the above time series. Their presence is due to the nature of the driving in the model. These memory effects are destroyed as soon as a large earthquake resets the system. |
Extent | 7093143 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-03-31 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0085652 |
URI | http://hdl.handle.net/2429/6675 |
Degree |
Doctor of Philosophy - PhD |
Program |
Physics |
Affiliation |
Science, Faculty of Physics and Astronomy, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1997-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1997-250628.pdf [ 6.76MB ]
- Metadata
- JSON: 831-1.0085652.json
- JSON-LD: 831-1.0085652-ld.json
- RDF/XML (Pretty): 831-1.0085652-rdf.xml
- RDF/JSON: 831-1.0085652-rdf.json
- Turtle: 831-1.0085652-turtle.txt
- N-Triples: 831-1.0085652-rdf-ntriples.txt
- Original Record: 831-1.0085652-source.json
- Full Text
- 831-1.0085652-fulltext.txt
- Citation
- 831-1.0085652.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0085652/manifest