- Email: [email protected]
Fractal, multifractal and sliding window correlation dimension analysis of sedimentary time series
Computers & Geosciences 25 (1999) 1009±1021
Fractal, multifractal and sliding window correlation p dimension analysis of sedimentary time series A. Prokoph* University of Ottawa, Department of Geology, P.O. Box 450, Stn. A, Ottawa, Ontario, Canada, K1N 6N5 Received 20 April 1998; received in revised form 22 August 1998; accepted 18 September 1998
Abstract Multifractal analysis, the box counting technique and sliding window correlation dimension analysis are methods for the detection and description of nonlinear variations in sedimentary time/depth series measurements such as well-log data. Combinations of these methods allow the detection of sedimentary discontinuities and fractal clustering of events such as bentonite layers or fossils horizons in sedimentary sections. The methods are used to analyze linear, periodic, fractal and random patterns with dierent sampling intervals and after setting distinctive measurement thresholds. Application of the methods is illustrated by means of computer simulation models and on the gamma-ray log from a Mid-Cretaceous sedimentary time series. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Time series analysis; Sedimentation models; Nonlinear; Fractal; Multifractal
1. Introduction Sedimentary time series lacking good lithological background are common in sedimentology and stratigraphy (e.g. logging data). Measurement of distinct sediment layers can be particularly time-consuming or impossible in weathered sections or in boreholes. Classical time series analysis methods, such as use of Markov chains, trend analysis and spectral analysis (Fourier or Walsh analysis) are widely used and explained in detail in Davis (1986) and Schwarzacher (1993). In recent years, wavelet analysis which permits the detection of signals in the frequency spectrum and
Code available at http://www.iamg.org/CGEditor/ index.htm * Tel.: +1-613-562-5800/6799; fax: +1-613-562-5192. E-mail address: [email protected] (A. Prokoph) p
in depth, (e.g. Prokoph and Barthelmes, 1996), and neural network analysis (Brescia et al., 1996) have been applied to geological time series. However, none of these methods are useful for the characterization of nonlinear systems. Sedimentary time series dier from physical measurements and experimental time series in several respects. 1. Only the preserved signals can be measured, not the original forcing due to physical processes. 2. The signal intervals can not be interpreted as timeintervals, because of ¯uctuating and often unknown sedimentation rates of the embedding rock (or ice). 3. Sedimentary processes are often nonstationary and discontinuous. This means that they can begin and end abruptly, e.g. by volcanic activity or subsidence in a sedimentary basin as illustrated in Fig. 1. Origin and preservation of the signals can also be in¯uenced by stochastic processes (e.g. Nicolis and Nicolis, 1991; Sanchez and Vos, 1994; Wiesenfeld
0098-3004/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 3 0 0 4 ( 9 9 ) 0 0 0 6 3 - 1
1010
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Fig. 1. Example of variations of simultaneous nonstationary processes with thresholds and dierent impacts on sedimentary successions. Signals of process 1 above threshold remain measurable but stops at certain time (e.g. vanishing volcanic activity); process 2 is measureable and increases its intensity linearily up to certain time and vanishes there abruptly; the intensity of process 3 declines continuously and it not measureable below threshold, but still exists.
and Moss, 1995) or postsedimentary erosion of the embedding rocks. 4. Data used for analysis are not equally spaced with respect to time in the original process, but set on the basis of economic restrictions and a priori considerations (e.g. about sedimentation rate). Therefore, time series data often lack sucient information for geological interpretation. 5. Data used may be completely or partially non-representative in regard to the underlying processes and patterns. Some of the nonlinear time series analysis methods are suitable for experimental time series with well-known parameters; for example, a nonlinear model can be advantageous for correlation dimension function analysis (Fowler and Roach, 1993). The objective of this paper is focused on the question whether or not nonlinear analysis techniques can detect, localize and parameterize deterministic patterns in complex sedimentary time series. Dierent methods of nonlinear analysis (multifractal spectrum, boxcounting and correlation dimension) are used to describe and detect characteristics of nonlinear, sinusoidal, random, and other nonstationary and discon-
Fig. 2. Experimental data sets from sedimentary time series models.
tinuous time series without relying on a priori assumptions about their dynamics. 2. Time series computer simulation models Generally, a process can be modeled by the dierential equation: dx=dt ÿdU=dx A sin ot F t where U is the potential bar, Asin(ot ) is the signal, and F is the noise (Wiesenfeld and Moss, 1995). Partially, the methods used allow for bar potentials as thresholds for binary (event) data in box-counting and multifractal analysis. The models including signal Asin(ot ) are without noise F, and those with noise
Table 1 Statistics of data sets generated by time series models and by gamma-ray log (BCCP group, 1994) time series Discontinuity Sinoidal function Random data Logistic function Gradual change Gamma-ray
UNCONF SINUS RANDOM FRAC39 MODSR GAMMA
number of data
mean
standard deviation
range
upper threshold value
1024 1000 1000 1024 1996 2275
0.7 0 0.5 0.5 0.01 74.79
20.902 20.708 20.287 20.296 21.06 29.36
ÿ1. . .+11 ÿ1. . .+1 0. . .+1 0. . .+1 ÿ3.72. . .+3.54 44.71±107.18
1.6 0.708 0.787 0.796 1.07 84.15
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
1011
Fig. 3. Distribution plots of modeled time series data.
without signal, to provide an easier understanding of the results. Five models were used to generate metrical as well as binary time series by selection of values greater than the mean plus one standard deviation in each case. Basic statistics are shown in Table 1. None of these models contain a superimposed linear trend. This is because trends generally can be removed by trend analysis (Davis, 1986). All model data are illustrated in Fig. 2. 2.1. SINUS model This model supports continuous periodic forcing of sedimentation at constant accumulation rates due to periodic ¯uctuations of the Earth's orbital parameters involving annual cycles and Milankovitch cycles. The model data simply are de®ned as X(i )=sin(2pt/250), i = 1 . . . 1000, with four sinusoidal curves. 2.2. RANDOM model Random processes (e.g. white noise) are common in sedimentary time series because of stochastic measurement errors and stochastic processes of sediment accumulation and composition. Results for this model are produced using the random generator in Microsoft EXCEL 5.0, with i = 1 . . .1000. Consequently, the values X(i ) are uniformly distributed throughout the entire data range (Fig. 3(A)).
2.3. FRAC39 model This model serves to illustrate the preservation of feedback and `self-organized criticality' processes in sedimentary sequences, which are related to the nonlinear models of May (1976) or Bak et al. (1988). The equation used is de®ned by the logistic map x(i + 1)=mx(i )(1ÿx(i )) with i = 1 . . . 1024. For values of m > 3.57, the process model becomes (deterministically) chaotic. The value of m=3.9 was used for application and the initial value was set equal to x(0)=0.6. The logistic model has been studied extensively in experiments with other nonlinear analysis methods (Fowler and Roach, 1993; Baker and Gollub, 1990). Consequently, most results obtained here are not new, but they are helpful in explaining the relative advantages and disadvantages of dierent analysis techniques. The single feedback process of the logistic map provides an almost uniform distribution throughout the entire data range with only a slight over-presentation of the very high (>0.9) and very low (<0.2) values (Fig.3(B)). 2.4. UNCONF model This model data may be used to describe abrupt temporal events in the sedimentary succession without in¯uence on later sedimentation (e.g. bentonite layers), abrupt increases/decreases in the sedimentation rate, which may be caused by abrupt changes in subsidence or sediment in¯ux and by abrupt facies changes (e.g.
1012
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Fig. 4. Correlation dimension analysis: power-law relationship between correlation C(r ) and eective data range r of (A) FRAC39 model and (B) MODSR model; saturation and depopulation limits used in estimation are marked by dotted lines; (C) correlation dimensions d due to embedding dimension (E ) for FRAC39 (q), MODSR (r) models and white noise (r).
from limestone/marl-cycles to sandstone/siltstonecycles), which may be caused by changes in the composition of terrigenous in¯ux or carbonate productivity. The frequency distribution of the data is characterized predominantly by four data classes at X(i )=(ÿ1,ÿ0.8), (0,0.2), (0.8,1) and (1.8,2) as shown in Fig. 3(C). The model used is described by
original data information. This eect has been described by Kowalewski (1996) as `time-averaging'. In the model, three high-frequency, equal-amplitude, sinusoidal signals (period=20, 40 and 100) are transformed with respect to a low-frequency sinusoidal transformation of accumulation (sedimentation) rate (period=1000 data). Therefore
x i sin 2pi=100
x t sin 2pt=20 sin 2pt=40 sin 2pt=100
x i sin 2pi=50
i 700 . . . 1024 i 300 . . . 699
x i sin 2pi=50 10 i 500 x i sin 2pi=50 ÿ 0:5 i 1, . . . 300: 2.5. MODSR model The complex model describes the preservation of signals transformed from time-dependency (dt ) to depthdependency (ds ). In its simplest form, y(s )=ax(t )dt. The resulting time series y(s ) was transformed to become equidistant with Ds = 1 by means of a simple weighted average calculation. Therefore, some of the new sample intervals contain an increased number of original data, whereas other intervals completely lack
and
y s Sx t cos 2pt=500 dt with t = 1 . . . 2000 and s = 1 . . . 1996. For the analysis, the variable s is replaced by discrete values i. After transformation, the data are quasi-normally (Gaussian) distributed (Fig. 3(D)). Because of the cosine function used, this model may describe gradual changes in sedimentary sections, such as thinning-upward or thickening-upward successions. The models UNCONF and MODSR have been previously used by Prokoph (1997) to demonstrate the preservation of periodic cycles and discontinuities by wavelet analysis. 3. Sliding window correlation dimension To evaluate the sensitivity of variables at dierent stages of the depositional process, the correlation
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
1013
dimension function plot (Grassberger and Procaccia, 1983; Nicholis and Nicholis, 1984) was used to estimate the dimensionality of the time series (or degrees of freedom of the system). The analytical procedure of correlation dimension analysis, including the transformation of a time series x(t ) into an array of E-dimensional vectors Xi, was already described by Nicolis and Nicolis (1984, p. 531), and Mudelsee and Stattegger (1994). The correlation function of the attractor C(r ) is estimated by
C r
N 1 X y r ÿ jXi ÿ Xj j 2 N i, j1
where y(x )=0 if x < 0 and y(x )=1 if x > 0. Relatively small C(r ) varies approximately as C(r )=rd implying that the dimensionality d of the attractor is given by the slope of log C(r ) versus log r within the restricted range of values of r where log C(r ) 0 dvlog rv. If d reaches its saturation limit beyond some relatively small embedding dimension E, then the value of E provides the minimum number of variables necessary to model the behavior represented by the attractor. For this study, the interval of distance r was selected empirically and set equal to a value between rmin greater than the probable measurement error and rmax equal to half of the data range R/2. Twenty distance values of r were chosen according to a logarithmical interval and the C(r )-interval from 0.001 . . .0.1 was automatically selected for calculation of d. This is because the C(r )-interval is relatively far from the saturation limit as well as from the `depopulation zone' without measurable power-law relationship (Fowler and Roach, 1993). Fig. 4 gives the dependence of C(r ) versus r in a log±log plot, and straight lines in the plot of model FRAC39 indicate a relatively strong power law relationship for scale invariance. Curved and discontinuous placements of C(r ) for model MODSR characterize absent or limited-range power law relationships between C(r ) and distance r. Correlation dimensions become saturated at d = 1.06 in FRAC39 and at d = 2.7 in MODSR, respectively. Notice that, initially, a single process controls the model FRAC39 and three processes control the model MODSR. These correlation dimension estimations are well established and allow the detection of the underlying nonlinear stationary processes. For the detection of nonstationary processes, a sliding window algorithm is applied where a distinctive analysis window size o < N is shifted in discrete analysis steps s over the whole data set. The ®nal sliding window correlation is de®ned as
Fig. 5. Sliding window correlation dimension plots for six sedimentary time series models.
C r, s
so=2 X 1 y r ÿ jXio ÿ Xjo j 2 o i, j1sÿo=2
The length of the time series analyzed is much less than the entire time series and the accuracy of the estimation is limited (cf. Essex, 1991). The analysis windows contain 200 data points for every calculation. Fig. 5 shows the temporary variations of the correlation dimension analysis of ®ve dierent time-series models. Notice that the discontinuities in the model UNCONF at i = 300 (`facies change') and at i = 700 (`abrupt doubling of sedimentation rate') are clearly visible. The correlation dimension of the model SINUS was set equal to 1. Additionally, increased correlation dimensions of 1.25 occur periodically in the 250 data intervals used. This increase is a result of insucient window size, which provides a slight arti®cial trend for intervals with dx=di max dx=di cos 2pi=250 in this model. The variations of the correlation dimension of the model MODSR result from the same transformation. In the model SINUS the smallest correlation dimension (2.5) occurs at i = 900. . . 1100
1014
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
for the highest sedimentation rate and, therefore, the smallest number of external periodic forcing in the model data. Intervals with the highest change of sedimentation rate have the highest correlation dimension (3.5). These intervals occur in the transition zone from condensed to diluted sedimentation. The sliding window correlation dimension of the single nonlinear process model FRAC39 varies between 1.2 and 1.6, providing a rather inaccurate estimate compared to the single window correlation dimension of 1 (Fowler and Roach, 1993). In contrast, the random noise model RANDOM provides a correlation dimension of approximately 7 with an embedding dimension of 9. This means that the saturation limit can not be reached because of the stochastic noise. Many of the preceding results were already obtained by Fowler and Roach (1993). However, in our case, samples with time-averaged information are studied. For example, a single nonlinear process controls sedimentation with annual cycles (e.g. preservation of winter (sand)/summer (clay)-cycles in lake strata), but samples or logging measurements contain time-averaged information for ®ve years periods (e.g. mixtures of sand and clay in single samples). Such primary changes can also be induced by bioturbation. Consequently, the resulting sampling interval of model FRAC39 becomes xz
zN=5 X
1 xj 5 jzN=zÿ5
with z=(1 . . . N/5). The correlation dimension of the time-averaged model FRAC39 rises to 6.5 and remains unsaturated at this value. Consequently, the time-averaged sampled time series shows stronger resemblance to stochastic noise (random data) than to the simple nonlinear process originally applied.
4. Box-counting analysis Box-counting analysis is a fractal method used for detection of self-similar distributions of event data (Mandelbrot, 1983; Feder, 1988). The Cantor set is one of the most used models to explain the characteristics of fractal distribution in sedimentary events (e.g. Plotnick, 1986). Attempts of applications of this technique to one-dimensional data include events such as turbidite deposition (Rothman et al., 1994), bedload transport (Ivanov, 1996), and fracture patterns (Walsh and Watterson, 1993). In practice, continuous time series do not provide event data. Events can, however, be embedded in these time series, such as bentonite layers or fossil horizons in clay/limestone successions.
Fig. 6. Box-counting analysis for RANDOM model with (A) relation between ®lled boxes N(e ) and box sizes e and (B) scaling behavior in power-law plot of fractal dimension d versus box sizes e.
In this paper the box-counting method is used to demonstrate the distribution of arti®cial events in the form of binary data, which exceed distinctive thresholds in the models used. For ease of understanding and recalculation, the range between Xmax and the mean plus one standard deviation is used to create `events'. The box-counting analysis supports detection of a single (unifractal) process and distribution pattern. The method is less useful for the analysis of distribution of events in time series, for which the in¯uence of more than a single process for the creation of a distinctive pattern (e.g. fossil horizons) is known in advance. Mathematical fractals show perfect self-similarity, although in reality, the geometry of fractals is never perfectly self-similar at all scales. Real life fractals are considered self-similar in a statistical sense. The idea of box-counting dimension is that the stratigraphic section of interest is covered with boxes of a ®xed size, and the number of boxes with data (part of the attractor) is recorded against box size; the process being repeated for dierent size boxes, and Ne Aeÿd where N = the number of boxes of a certain size (side length e ) which have point data occurrences in them; e=the length of box side; d = box counting dimension.
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
1015
Fig. 7. Box-counting analysis for time series models with relations between ®lled boxes N(e ) and box size e.
In this study, the number of event occurrences N(e ) in intervals (boxes) of 19 dierent sizes e has been calculated for all ®ve model datas. The box sizes range from 1 (original interval size) to the whole data set. Usually, the relationship between N(e ) and e is characterized by three stages: 1. Saturation occurs when events exist in all boxes (intervals) and an increase of the box size (which is inversely proportional to a decrease in the number of boxes) provides the same linear increase in the number of events. 2. Nonlinear functionality occurs when the number of boxes that contain events increases slower than the number of boxes according to a functional dependence. This scaling interval may be characterized by fractal scale invariant growth. 3. Depopulation occurs when the number of events increases much slower than the number of boxes or becomes stagnant. Results of box-counting analysis for the model RANDOM are shown in Fig. 6. The following three parameters are involved in the characterization of the scaling behavior of the data. 1. The relationship between N(e ) and box size e provides the number of `®lled' boxes. 2. The relationship of the ratio N(e )/Ng(e ) versus e gives the probability of getting a `®lled' box at a dis-
tinctive interval size e where Ng(e ) represents total number of boxes of size e. If N(e )/Ng(e )=1 there is complete saturation with events occurring in all boxes at least once. 3. The slope of D log(N(e )/log(e ) versus e gives the box-size range with power-law behavior; it means that scale-invariance is obtained for this box range. The interval used is relatively small in this example with box sizes e=4.5 . . . 8.6 and with a box-counting dimension d=ÿ0.7, selected to characterize stochastic random time series. Whereas the determination of the box size of saturation is exact, the determination of the boundary between depopulation and nonlinear function is subject to determination by the researcher. Fig. 7 shows power-law plots and growth probability of N(e ) versus e for events greater than mean plus one standard deviation in the time series model data. The nonlinear model FRAC39 can be detected only in box sizes between 3 and 6 with d=ÿ0.7 indicating a weak power-law relationship for these event-like data less than one scaling interval. The models SINUS, MODSR and UNCONF are characterized by stepwise increase of the probability function N(e )/Ng(e ), where each step may represent one of the cycle periodicities in the model data. The saturation limits of models SINUS and UNCONF are for 250 and 50 data, re-
1016
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Fig. 8. Multifractal spectra f(a ) versus a for ®ve continuous and event-like time series models, respectively.
spectively, representing major periodicities in the time series. Power-law relationships are suitable for describing the cycle-dominated models SINUS and UNCONF for distinctive box size ranges with d=ÿ0.93 and d=ÿ0.82, respectively (Fig.7). The cause of such scaling behavior is not mathematically described here and the reason why a periodic function can produce scale invariance is not fully understood.
5. Multifractal analysis In contrast to the fractal box-counting technique used for the study of binary data sets, multifractals analyze self-similar patterns using measures of variables such as those provided by continuous time series. Originally, multifractal analysis was introduced to describe growth probability distribution in kinetic aggregation processes where the growth probability for the perimeter sites can be regarded as a measure associated with each site (Amitrano et al., 1986). In this study, the sum of X(i ) was scaled to 1 by setting X(i )=X(i )/SX. The measures (data X(i )) in these intervals were treated as probabilities p of occurrence. Two possible multifractal patterns may occur in discretized continuous sedimentary time series: 1. One-dimensional features of spatially intertwined fractals coexisting in the same time series, such as
Cantor bars with dierent sizes (masses) (cf. Feder, 1988, pp. 64±65). 2. Nonstationarities of fractal (or multifractal) patterns in a time series caused by nonstationary nonlinear processes (e.g. volcanic or storm events). The following method can be used to determine a multifractal spectrum and ascertain if the data are best described by a fractal or multifractal distribution (Feder, 1988): An estimate of the amount of speci®c variable (e.g. logging value) per sample interval for a number of boxes that cover the time series of interest is obtained. All samples X(i ) within a box (interval) are used to obtain an average value for all boxes used. The measure for X(i ) per box can be set equal to this average value multiplied by e. Given a function wq(e ) for any real number q and box size e: wq e
N e X mqi e: i1
The measure mi(e ) satis®es a multifractal model if wq(e ) forms a power law relationship with e for any value of q, or wq(e )0 e t(q ) where t(q ) represents an auxiliary function. The HoÈlder exponent a satis®es a=d{t(q )}/dq, and can be estimated by using the Legendre-transformation (Feder, 1988) using a ®ve-point window. This procedure allows the multifractal spectrum f(a ) to be plotted according to f(a )=qaÿt(q ). Because we are
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
1017
interested in the scaling properties of the multifractal, it is desirable to know how wq scales with respect to e for a continuous range of q values with wq eAe qÿ1D q (Gould and Tobochnik, 1990). D(q ) is the generalized dimension of order q or dimension function introduced by Grassberger and Procaccia (1983). It can be estimated as D(q )=t(q )/(1ÿq ), q $ 1 (Feder, 1988). An unifractal distribution exists, if D(q ) versus q gives a straight line. Feder (1988) discussed that amin=ln p0/ln e0 and amax=ln p1/e1, with maximum probability or mass p0=Xmax/SX and minimum probability p0=Xmin/SX for associated intervals l0 and l1, respectively. For sedimentary successions, e1 and e0 are given by the thickness of layers with speci®c measurements such as grain size, not provided by the previous computer models in this paper. Relatively small HoÈlder exponents a support large values of X(i ), which tend to be clustered, whereas high a supports small values of X(i ) dispersed over the time series. For q = 0, all measurements become Xq(i )=1, indicating equal probability for all data in each scaling interval e. Then the multifractal spectrum reduces to fmax=d where d is the boxcounting dimension. For each calculation, six scaling intervals (segments) I for q=ÿ21, . . .,+21 were used, with I1=N/2, and Ii + 1=Ii/2, until I6=N/26, and N representing the size of the time series; e.g. 1000 values were used in model SINUS. Therefore, the number of segments ranged from 2 to 64. The distributions of measurements in very small boxes were not taken into account; Xmin and Xmax do not apply to the original data range, but represent averaged maxima and minima for I6. The multifractal spectra f(a ) versus a for all models are shown in Fig. 8. In general, the multifractal spectrum of the `upper threshold' time series considers more HoÈlder exponents a than the complete and continuous time series model. In each model fmax=1, which is a result of the saturation of the values for small box size e, as previously shown in Fig. 6. The spectrum f(a ) for model FRAC39 is close to a single spike characteristic of a unifractal distribution with fmax=1. The same applies to model RANDOM. Both models provide a nearly uniform distribution for each class of measurements (see Fig. 3). These distributions imply averaging of Xmin and Xmax at small e, resulting in small variances in the singularities a of the measurements. Relatively high f(a ) for small a is seen for model FRAC39 (upper threshold) whereas the distribution of f(a )/a is almost symmetrical in the model RANDOM (upper threshold). This dierence is caused by the relatively large number of large measurements X(i ) in model
Fig. 9. Dimension function plots D(q ) versus q for (A) three continuous and (B) ®ve event-like model data time series models.
FRAC39 in comparison with model RANDOM in which the measurements are uniformly distributed (see Fig. 3). In contrast, the models UNCONF, MODSR and SINUS show broad multifractal spectra for complete as well as `upper-threshold' data. There is no comparable multifractal spectrum for model SINUS because its completely periodic distribution does not provide a slope for a(q )=ÿdt(q )/dq. In model MODSR (upper part) at both ends its multifractal spectrum f(a ) is small (0.05) and the growth probability for the measurements with increasing box numbers I is small as well. The dimension function D(q )/q plot for the model UNCONF shows a large value of D(q ) for small q. This means that measurements X(i ) from just above the threshold value are widely dispersed, whereas high values, such as the single event at X(300)=11 results in a small dimension D(q ) (Fig. 9). The model RANDOM shows constant dimension both in the complete time series as well as in the `upper-threshold' time series. This is a result of the randomness of the process underlying this model. The model FRAC39 also shows approximately equal values
1018
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Fig. 11. Distribution plot of gamma-ray data of Kirchrode 1/ 91 well.
Fig. 10. Correlation dimension analysis for (A) gamma-ray data of Kirchrode 1/91 well (Upper Albian of North Germany); gamma-ray readings in API; (B) sliding windows correlation dimension; note the small variation of correlation dimension at embedding dimensions between 3 and 4; (C) D correlation dimension marks stability of underlying 3- to 4dimensional system; x marks most obvious sedimentary discontinuity (fossiliferous mud-turbidite horizon) with high D correlation dimension.
for all D(q ), which may be a result of time-averaging within the box sizes used for analysis. Computer programs in C for UNIX (TIMEMF.C for multifractal analysis, TIMESWCD.C for sliding window correlation dimension analysis) or TurboPascal for PC (FRACZEIT.PAS for box-counting analysis) are available on the IAMG server (ftp.iamg.org). Input and output ®les are generally ASCII-text ®les. Choice of input parameters such as window size, data range, box size range and box numbers is controlled interactively. 6. Case study: gamma-ray log, Kirchrode 1/91 (Cretaceous, North Germany) A gamma-ray log consisting of 2275 continuous measurements of marine marlstones was used for this case study (Fig. 10, Table 1). The section covers about 2.5 Ma representing almost the entire Late Albian (Cretaceous) in the North German Basin (BCCP Group, 1994). The cycle pattern in time and frequency was previously studied by wavelet analysis showing well-pronounced 12 m cyclicity in the upper part (40±165 m) and 18 m cyclicity in the lower part (178±216 m) of the section (Prokoph and Barthelmes, 1996). Mud-turbidites at the base and top of the section probably mark short-time hiatusses
within the corresponding intervals (28±8 m, 235± 210 m). The sliding window correlation dimension (200 data per window=20 m) of the gamma-ray data varies between 3.2 at 180 m and 4.1 at 80 m, possibly indicating a complex control of sedimentation by as many as three to four dierent factors in the section (Fig. 10(B),(C)). The saturation of the correlation function is shown by its ®rst derivative D, which reaches its maximum of 0.3 at about 220 m, where an hiatus is most likely to occur (BCCP Group, 1994; Prokoph, 1994). This high ®rst derivative indicates that chaotic randomness is most likely to occur in the interval 210± 230 m. In contrast, the interval from about 40 m to 170 m has D correlation dimension of less than 0.15. The Milankovitch cyclicity found at depths from 160 to 40 m (Prokoph and Barthelmes, 1996) probably was the predominant control of the sedimentation. The data are nearly normally (Gaussian) distributed with most data concentrated in a class around the average (70±80 API) as shown in Fig. 11. Both data distribution and correlation dimension are similar to those of model MODSR suggesting that nonlinear changes of sedimentation rate were at least temporarily superimposed by Milankovitch cyclicity. The box-counting technique shows a power-law relationship in the scaling interval e from 5 to 50 (=0.5± 5 m) with a box-counting dimension of d=ÿ0.5807 and saturation limit N(e )/Ng(e )=1 at box size 300 (30 m) (Fig. 12). However, this estimate can not be extrapolated to narrow data intervals. For example, 429 ®lled boxes predicted for box size 1 (=10 cm) exceed the 373 events counted by about 13% (Fig. 12(A)). This discrepancy is a result of depopulation. The multifractal spectrum of the gamma-ray data shows a spike-like curve of f(a )/a with f(a ) between 0.7 and 1 (Fig. 12(C)). Therefore, all classes of measurements X(i ) exhibit almost the same scaling
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
1019
Fig. 12. (A), (B) results of box-counting analysis, and (C) multifractal spectra of gamma-ray data of Kirchrode 1/91 well; r mark continuous data, mark time series data below mean-standard deviation, r mark time series data above mean+standard deviation.
behavior. On the other hand, the event-like data from the upper and lower threshold shows relatively high f(a ) > 0.2 for small a < 0.8. This may be a result of power-law relationships involving high gamma-ray values, probably exceeding 100 API. In contrast, gamma-ray values just above the threshold (about 85 API) are dispersed in the series and provide a > 2. 7. Concluding remarks Sliding window correlation dimension provides a method for the detection of temporary changes of the number of process controlling parameters through time (depth) and, therefore, allows the detection of sedimentary discontinuities such as facies changes, events and unconformities better than spectral and wavelet-analysis. Multifractal analysis is shown to be a tool for discrimination between unifractal, random and multiple time-related data patterns and for interpretation of the scaling behavior of classes of measurements in continuous time series. However, a parameterization of processes that controlled the data distribution shown in the multifractal spectrum and the dimension function plot is probably only possible for very high (Xmax) and very low values (Xmin) in the time series.
Box-counting techniques permit the discrimination of saturation and depopulation limits (sample intervals) of event-like data and the calculation of the scaling behavior between the corresponding two box size (interval length) limits. It can also give the box-counting dimension for power-law functionality between event clustering (`®lled' boxes) and the sampling interval size of observation (box sizes), results similar to those provided by lacunarity analysis (Plotnick et al., 1996). In this way it is possible to predict the maximum interval between events and probably the cluster characteristics of these events both in sedimentary sections and lateral distribution. Therefore, the ideas of Plotnick (1986) about the fractal distribution of stratigraphic hiatuses and Sadler's (1981) prediction of power-law relationship between sedimentation rate and observation time may also be related to scaling intervals between depopulation and saturation. Log±log regression with errors of less than 0.2% do not permit the determination of signi®cant scale invariance and subsequent prediction of sedimentary events for dierent sampling intervals. The time-event distribution patterns are likely to be expressed in their spatial distribution (i.e. sedimentary basin) as discussed by Yamaji (1988). Fig. 13 shows a model of the importance of saturation and depopulation limit as well as nonlinear func-
1020
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Fig. 13. Two-dimensional basin model for recognition of events due to dierent sedimentation rates and scaling (observation) intervals. For further explanation see text.
tionality in the distribution of discrete events of the same type (e.g. bentonite layers, fossil horizons, sequence boundaries) due to observation limits (measurement interval). For example, a sedimentary succession with four unevenly distributed volcanic ash layers (bentonites) in time and depth covers a time interval of 1 million years. The maximum resolution of a certain measurement method (e.g. gamma-ray logging) permits a resolution of 1 m-intervals. Therefore, in case the succession deposited during the 1 million years interval is only 1 or 2 m thick, every measurement detects a bentontite, i.e. the measurement is above the saturation limit, and the probability of detection of a bentonite layer is equal to 1. With increasing non-bentonite background sedimentation (e.g. marlstones) towards the basin centre, the number of measurements without bentonite signals (e.g. positive gamma-ray value) decreases and, ®nally, measurements represent only a single bentontite layer each. The depopulation limit is reached and the probability of ®nding more bentonite layers approaches 0 for increasing sedimentation rate. For analysis of nonlinear and nonperiodic continuous sedimentary time/depth series the following algorithms are suggested: 1. Sliding window correlation dimension calculation to
detect temporal discontinuities and the approximate number of underlying processes that probably controlled the dierent populations of the data groups (e.g. sedimentary pattern such as lamination and turbidite layers). 2. Truncation of time series at distinct discontinuities. 3. Multifractal analysis of segment cuts from the original time series by using threshold values to distinguish between unifractal, multifractal and nonfractal (linear, periodic) distribution. Estimation of growth probabilities at amin, a0 and amax. 4. Box-counting analysis of event-like data from the cut time series to estimate depopulation and saturation limits as well as nonlinear relationships in the restricted scale range. Together with trend and wavelet analysis, these techniques are tools to improve the knowledge about the temporary and spatial distribution of sedimentary patterns and processes. Computer programs of time-series analysis for box-counting analysis (FRACZEIT.PAS) in TurboPascal for PC, multifractal analysis (TIMEMF.C) and sliding window correlation dimension analysis (TIMESWCD.C) in C for UNIX are available on the Editor's Home Page (http://www.iamg.org/CGEditor/ index.htm).
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Acknowledgements I am grateful for discussions and reviewing dierent versions of this article by Walter Schwarzacher, Karl Stattegger, Bradley Sim and Frits Agterberg. The German Research Council (DFG) supports my research in Grant Pr 545/1-1/2. References Amitrano, C., Coniglio, A., di Liberto, F., 1986. Growth probability distribution in kinetic aggregation process. Physical Review Letters 57 (8), 1016±1019. Bak, P., Tang, C., Wiesenfeld, K., 1988. Self-organized criticality. Physical Review A38, 364±374. Baker, G.L., Gollub, J.P., 1990. Chaotic Dynamics: an Introduction. Cambridge University Press, Cambridge, p. 182. BCCP Group, 1994. The Upper Albian of northern Germany: results from the Kirchrode 1/91 borehole, Boreal Cretaceous Cycles Project (BCCP). Zentralblatt fuÈr Geologie und PalaÈontologie 1 (7/8), 809±822. Brescia, M., D'Argenio, B., Ferreri, V., Longo, G., Pelosi, N., Rampone, S., Tagliaferri, R., 1996. Neural net aided detection of astronomical periodicities in geological records. Earth & Planetary Science Letters 139, 33±45. Davis, J.C., 1986. Statistics and Data Analysis in Geology. John Wiley & Sons, New York, p. 646. Essex, C., 1991. Correlation dimension and data sample size. In: Schertzer, D., Lovejoy, S. (Eds.), Non-linear Variablility in Geophysics, Scaling and Fractals. Kluwer, Dordrecht, pp. 93±98. Feder, J., 1988. Fractals. Plenum Press, New York, p. 283. Fowler, A.D., Roach, D.E., 1993. Dimensionality analysis of time series data: nonlinear methods. Computers & Geosciences 19 (1), 41±52. Gould, H., Tobochnik, J., 1990. More on fractal and chaos: multifractals. Computers in Physics 3/4, 202±207. Grassberger, P., Procaccia, I., 1983. Measurement of the strangeness of strange attractors. Physica 9D, 189±208. Ivanov, S., 1996. Variability of sedimentary sequence: numerical modeling of the deposition-erosion process. Geologische Rundschau 85, 12±18. Kowalewski, M., 1996. Time-averaging, overcompleteness, and the geological record. Journal of Geology 104, 317± 326. Mandelbrot, B., 1983. The Fractal Geometry of Nature. W. H. Freeman, Freeman, p. 468. May, R.M., 1976. Simple mathematical models with very complicated dynamics. Nature 261 (5560), 459±467. Mudelsee, M., Stattegger, K., 1994. Application of the
1021
Grassberger±Procaccia algorithm to the d18O Record from ODP Site 659: selected methodical aspects. In: Kruhl, J.H. (Ed.), Fractals and Dynamic Systems in Geoscience. Springer, Berlin, pp. 399±421. Nicolis, C., Nicolis, G., 1984. Is there a climate attractor? Nature 311 (5986), 529±532. Nicolis, G., Nicolis, C., 1991. Nonlinear dynamic systems in the geosciences. In: Fransen, E.K., Watney, W.L., Kendall, C.G.StC, Ross, W. (Eds.), Sedimentary Modeling: Computer Simulations and Methods for Improved Parameter De®nition, vol. 233. Kansas Geological Survey Bulletin, pp. 33±42. Plotnick, R.E., 1986. A fractal model for the distribution of stratigraphic hiatuses. Journal of Geology 94 (6), 885±890. Plotnick, R.E., Gardner, R.H., Hargrove, W.W., Prestegaard, K., Perlmutter, M., 1996. Lacunarity analysis: a general technique for the analysis of spatial pattern. Physical Review E53 (5), 1±8. Prokoph, A., 1994. Zyklische Sedimentation im Oberalb des Norddeutschen Beckens. Tuebinger Geowissenschaftliche Arbeiten A19, 163. Prokoph, A. 1997. New mathematical tools and their applications in quantitative stratigraphy: Detection of cycles, chaos and unconformities by wavelet-and fractal analysis. In: Proceedings of the IAMG'97 Conference, Barcelona, vol 2, pp. 507±512. Prokoph, A., Barthelmes, F., 1996. Detection of nonstationarities in geological time series: Wavelet transform of chaotic and cyclic sequences. Computers & Geosciences 22 (10), 1097±1108. Rothman, D.A., Grotzinger, J.P., Flemings, P., 1994. Scaling in turbidite deposition. Journal Sedimentary Research A64 (1), 59±67. Sadler, P.M., 1981. Sediment accumulation rates and the completeness of stratigraphic sections. Journal of Geology 89, 569±584. Sanchez, A., Vos, H., 1994. From solar variability to climate proxy data Ð some ideas about signal transformation and the interpretation of time series. Terra Nostra 1/94, 19±24. Schwarzacher, W., 1993. Cyclostratigraphy and Milankovitch Theory. In: Developments in Sedimentology, vol. 52. Elsevier, Amsterdam, p. 225. Walsh, J.J., Watterson, J., 1993. Fractal analysis of fracture patterns using the standard box-counting technique: valid and invalid methodologies. Journal of Structural Geology 15 (12), 1509±1512. Wiesenfeld, K., Moss, F., 1995. Stochastic resonance and the bene®ts of noise: from ice ages to cray®sh and SQUIDs. Nature 373, 33±36. Yamaji, A., 1988. A fractal model for the distribution of stratigraphic hiatuses: a discussion. Journal of Geology 96 (6), 101±102.
Fractal, multifractal and sliding window correlation p dimension analysis of sedimentary time series A. Prokoph* University of Ottawa, Department of Geology, P.O. Box 450, Stn. A, Ottawa, Ontario, Canada, K1N 6N5 Received 20 April 1998; received in revised form 22 August 1998; accepted 18 September 1998
Abstract Multifractal analysis, the box counting technique and sliding window correlation dimension analysis are methods for the detection and description of nonlinear variations in sedimentary time/depth series measurements such as well-log data. Combinations of these methods allow the detection of sedimentary discontinuities and fractal clustering of events such as bentonite layers or fossils horizons in sedimentary sections. The methods are used to analyze linear, periodic, fractal and random patterns with dierent sampling intervals and after setting distinctive measurement thresholds. Application of the methods is illustrated by means of computer simulation models and on the gamma-ray log from a Mid-Cretaceous sedimentary time series. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: Time series analysis; Sedimentation models; Nonlinear; Fractal; Multifractal
1. Introduction Sedimentary time series lacking good lithological background are common in sedimentology and stratigraphy (e.g. logging data). Measurement of distinct sediment layers can be particularly time-consuming or impossible in weathered sections or in boreholes. Classical time series analysis methods, such as use of Markov chains, trend analysis and spectral analysis (Fourier or Walsh analysis) are widely used and explained in detail in Davis (1986) and Schwarzacher (1993). In recent years, wavelet analysis which permits the detection of signals in the frequency spectrum and
Code available at http://www.iamg.org/CGEditor/ index.htm * Tel.: +1-613-562-5800/6799; fax: +1-613-562-5192. E-mail address: [email protected] (A. Prokoph) p
in depth, (e.g. Prokoph and Barthelmes, 1996), and neural network analysis (Brescia et al., 1996) have been applied to geological time series. However, none of these methods are useful for the characterization of nonlinear systems. Sedimentary time series dier from physical measurements and experimental time series in several respects. 1. Only the preserved signals can be measured, not the original forcing due to physical processes. 2. The signal intervals can not be interpreted as timeintervals, because of ¯uctuating and often unknown sedimentation rates of the embedding rock (or ice). 3. Sedimentary processes are often nonstationary and discontinuous. This means that they can begin and end abruptly, e.g. by volcanic activity or subsidence in a sedimentary basin as illustrated in Fig. 1. Origin and preservation of the signals can also be in¯uenced by stochastic processes (e.g. Nicolis and Nicolis, 1991; Sanchez and Vos, 1994; Wiesenfeld
0098-3004/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 0 9 8 - 3 0 0 4 ( 9 9 ) 0 0 0 6 3 - 1
1010
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Fig. 1. Example of variations of simultaneous nonstationary processes with thresholds and dierent impacts on sedimentary successions. Signals of process 1 above threshold remain measurable but stops at certain time (e.g. vanishing volcanic activity); process 2 is measureable and increases its intensity linearily up to certain time and vanishes there abruptly; the intensity of process 3 declines continuously and it not measureable below threshold, but still exists.
and Moss, 1995) or postsedimentary erosion of the embedding rocks. 4. Data used for analysis are not equally spaced with respect to time in the original process, but set on the basis of economic restrictions and a priori considerations (e.g. about sedimentation rate). Therefore, time series data often lack sucient information for geological interpretation. 5. Data used may be completely or partially non-representative in regard to the underlying processes and patterns. Some of the nonlinear time series analysis methods are suitable for experimental time series with well-known parameters; for example, a nonlinear model can be advantageous for correlation dimension function analysis (Fowler and Roach, 1993). The objective of this paper is focused on the question whether or not nonlinear analysis techniques can detect, localize and parameterize deterministic patterns in complex sedimentary time series. Dierent methods of nonlinear analysis (multifractal spectrum, boxcounting and correlation dimension) are used to describe and detect characteristics of nonlinear, sinusoidal, random, and other nonstationary and discon-
Fig. 2. Experimental data sets from sedimentary time series models.
tinuous time series without relying on a priori assumptions about their dynamics. 2. Time series computer simulation models Generally, a process can be modeled by the dierential equation: dx=dt ÿdU=dx A sin ot F t where U is the potential bar, Asin(ot ) is the signal, and F is the noise (Wiesenfeld and Moss, 1995). Partially, the methods used allow for bar potentials as thresholds for binary (event) data in box-counting and multifractal analysis. The models including signal Asin(ot ) are without noise F, and those with noise
Table 1 Statistics of data sets generated by time series models and by gamma-ray log (BCCP group, 1994) time series Discontinuity Sinoidal function Random data Logistic function Gradual change Gamma-ray
UNCONF SINUS RANDOM FRAC39 MODSR GAMMA
number of data
mean
standard deviation
range
upper threshold value
1024 1000 1000 1024 1996 2275
0.7 0 0.5 0.5 0.01 74.79
20.902 20.708 20.287 20.296 21.06 29.36
ÿ1. . .+11 ÿ1. . .+1 0. . .+1 0. . .+1 ÿ3.72. . .+3.54 44.71±107.18
1.6 0.708 0.787 0.796 1.07 84.15
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
1011
Fig. 3. Distribution plots of modeled time series data.
without signal, to provide an easier understanding of the results. Five models were used to generate metrical as well as binary time series by selection of values greater than the mean plus one standard deviation in each case. Basic statistics are shown in Table 1. None of these models contain a superimposed linear trend. This is because trends generally can be removed by trend analysis (Davis, 1986). All model data are illustrated in Fig. 2. 2.1. SINUS model This model supports continuous periodic forcing of sedimentation at constant accumulation rates due to periodic ¯uctuations of the Earth's orbital parameters involving annual cycles and Milankovitch cycles. The model data simply are de®ned as X(i )=sin(2pt/250), i = 1 . . . 1000, with four sinusoidal curves. 2.2. RANDOM model Random processes (e.g. white noise) are common in sedimentary time series because of stochastic measurement errors and stochastic processes of sediment accumulation and composition. Results for this model are produced using the random generator in Microsoft EXCEL 5.0, with i = 1 . . .1000. Consequently, the values X(i ) are uniformly distributed throughout the entire data range (Fig. 3(A)).
2.3. FRAC39 model This model serves to illustrate the preservation of feedback and `self-organized criticality' processes in sedimentary sequences, which are related to the nonlinear models of May (1976) or Bak et al. (1988). The equation used is de®ned by the logistic map x(i + 1)=mx(i )(1ÿx(i )) with i = 1 . . . 1024. For values of m > 3.57, the process model becomes (deterministically) chaotic. The value of m=3.9 was used for application and the initial value was set equal to x(0)=0.6. The logistic model has been studied extensively in experiments with other nonlinear analysis methods (Fowler and Roach, 1993; Baker and Gollub, 1990). Consequently, most results obtained here are not new, but they are helpful in explaining the relative advantages and disadvantages of dierent analysis techniques. The single feedback process of the logistic map provides an almost uniform distribution throughout the entire data range with only a slight over-presentation of the very high (>0.9) and very low (<0.2) values (Fig.3(B)). 2.4. UNCONF model This model data may be used to describe abrupt temporal events in the sedimentary succession without in¯uence on later sedimentation (e.g. bentonite layers), abrupt increases/decreases in the sedimentation rate, which may be caused by abrupt changes in subsidence or sediment in¯ux and by abrupt facies changes (e.g.
1012
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Fig. 4. Correlation dimension analysis: power-law relationship between correlation C(r ) and eective data range r of (A) FRAC39 model and (B) MODSR model; saturation and depopulation limits used in estimation are marked by dotted lines; (C) correlation dimensions d due to embedding dimension (E ) for FRAC39 (q), MODSR (r) models and white noise (r).
from limestone/marl-cycles to sandstone/siltstonecycles), which may be caused by changes in the composition of terrigenous in¯ux or carbonate productivity. The frequency distribution of the data is characterized predominantly by four data classes at X(i )=(ÿ1,ÿ0.8), (0,0.2), (0.8,1) and (1.8,2) as shown in Fig. 3(C). The model used is described by
original data information. This eect has been described by Kowalewski (1996) as `time-averaging'. In the model, three high-frequency, equal-amplitude, sinusoidal signals (period=20, 40 and 100) are transformed with respect to a low-frequency sinusoidal transformation of accumulation (sedimentation) rate (period=1000 data). Therefore
x i sin 2pi=100
x t sin 2pt=20 sin 2pt=40 sin 2pt=100
x i sin 2pi=50
i 700 . . . 1024 i 300 . . . 699
x i sin 2pi=50 10 i 500 x i sin 2pi=50 ÿ 0:5 i 1, . . . 300: 2.5. MODSR model The complex model describes the preservation of signals transformed from time-dependency (dt ) to depthdependency (ds ). In its simplest form, y(s )=ax(t )dt. The resulting time series y(s ) was transformed to become equidistant with Ds = 1 by means of a simple weighted average calculation. Therefore, some of the new sample intervals contain an increased number of original data, whereas other intervals completely lack
and
y s Sx t cos 2pt=500 dt with t = 1 . . . 2000 and s = 1 . . . 1996. For the analysis, the variable s is replaced by discrete values i. After transformation, the data are quasi-normally (Gaussian) distributed (Fig. 3(D)). Because of the cosine function used, this model may describe gradual changes in sedimentary sections, such as thinning-upward or thickening-upward successions. The models UNCONF and MODSR have been previously used by Prokoph (1997) to demonstrate the preservation of periodic cycles and discontinuities by wavelet analysis. 3. Sliding window correlation dimension To evaluate the sensitivity of variables at dierent stages of the depositional process, the correlation
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
1013
dimension function plot (Grassberger and Procaccia, 1983; Nicholis and Nicholis, 1984) was used to estimate the dimensionality of the time series (or degrees of freedom of the system). The analytical procedure of correlation dimension analysis, including the transformation of a time series x(t ) into an array of E-dimensional vectors Xi, was already described by Nicolis and Nicolis (1984, p. 531), and Mudelsee and Stattegger (1994). The correlation function of the attractor C(r ) is estimated by
C r
N 1 X y r ÿ jXi ÿ Xj j 2 N i, j1
where y(x )=0 if x < 0 and y(x )=1 if x > 0. Relatively small C(r ) varies approximately as C(r )=rd implying that the dimensionality d of the attractor is given by the slope of log C(r ) versus log r within the restricted range of values of r where log C(r ) 0 dvlog rv. If d reaches its saturation limit beyond some relatively small embedding dimension E, then the value of E provides the minimum number of variables necessary to model the behavior represented by the attractor. For this study, the interval of distance r was selected empirically and set equal to a value between rmin greater than the probable measurement error and rmax equal to half of the data range R/2. Twenty distance values of r were chosen according to a logarithmical interval and the C(r )-interval from 0.001 . . .0.1 was automatically selected for calculation of d. This is because the C(r )-interval is relatively far from the saturation limit as well as from the `depopulation zone' without measurable power-law relationship (Fowler and Roach, 1993). Fig. 4 gives the dependence of C(r ) versus r in a log±log plot, and straight lines in the plot of model FRAC39 indicate a relatively strong power law relationship for scale invariance. Curved and discontinuous placements of C(r ) for model MODSR characterize absent or limited-range power law relationships between C(r ) and distance r. Correlation dimensions become saturated at d = 1.06 in FRAC39 and at d = 2.7 in MODSR, respectively. Notice that, initially, a single process controls the model FRAC39 and three processes control the model MODSR. These correlation dimension estimations are well established and allow the detection of the underlying nonlinear stationary processes. For the detection of nonstationary processes, a sliding window algorithm is applied where a distinctive analysis window size o < N is shifted in discrete analysis steps s over the whole data set. The ®nal sliding window correlation is de®ned as
Fig. 5. Sliding window correlation dimension plots for six sedimentary time series models.
C r, s
so=2 X 1 y r ÿ jXio ÿ Xjo j 2 o i, j1sÿo=2
The length of the time series analyzed is much less than the entire time series and the accuracy of the estimation is limited (cf. Essex, 1991). The analysis windows contain 200 data points for every calculation. Fig. 5 shows the temporary variations of the correlation dimension analysis of ®ve dierent time-series models. Notice that the discontinuities in the model UNCONF at i = 300 (`facies change') and at i = 700 (`abrupt doubling of sedimentation rate') are clearly visible. The correlation dimension of the model SINUS was set equal to 1. Additionally, increased correlation dimensions of 1.25 occur periodically in the 250 data intervals used. This increase is a result of insucient window size, which provides a slight arti®cial trend for intervals with dx=di max dx=di cos 2pi=250 in this model. The variations of the correlation dimension of the model MODSR result from the same transformation. In the model SINUS the smallest correlation dimension (2.5) occurs at i = 900. . . 1100
1014
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
for the highest sedimentation rate and, therefore, the smallest number of external periodic forcing in the model data. Intervals with the highest change of sedimentation rate have the highest correlation dimension (3.5). These intervals occur in the transition zone from condensed to diluted sedimentation. The sliding window correlation dimension of the single nonlinear process model FRAC39 varies between 1.2 and 1.6, providing a rather inaccurate estimate compared to the single window correlation dimension of 1 (Fowler and Roach, 1993). In contrast, the random noise model RANDOM provides a correlation dimension of approximately 7 with an embedding dimension of 9. This means that the saturation limit can not be reached because of the stochastic noise. Many of the preceding results were already obtained by Fowler and Roach (1993). However, in our case, samples with time-averaged information are studied. For example, a single nonlinear process controls sedimentation with annual cycles (e.g. preservation of winter (sand)/summer (clay)-cycles in lake strata), but samples or logging measurements contain time-averaged information for ®ve years periods (e.g. mixtures of sand and clay in single samples). Such primary changes can also be induced by bioturbation. Consequently, the resulting sampling interval of model FRAC39 becomes xz
zN=5 X
1 xj 5 jzN=zÿ5
with z=(1 . . . N/5). The correlation dimension of the time-averaged model FRAC39 rises to 6.5 and remains unsaturated at this value. Consequently, the time-averaged sampled time series shows stronger resemblance to stochastic noise (random data) than to the simple nonlinear process originally applied.
4. Box-counting analysis Box-counting analysis is a fractal method used for detection of self-similar distributions of event data (Mandelbrot, 1983; Feder, 1988). The Cantor set is one of the most used models to explain the characteristics of fractal distribution in sedimentary events (e.g. Plotnick, 1986). Attempts of applications of this technique to one-dimensional data include events such as turbidite deposition (Rothman et al., 1994), bedload transport (Ivanov, 1996), and fracture patterns (Walsh and Watterson, 1993). In practice, continuous time series do not provide event data. Events can, however, be embedded in these time series, such as bentonite layers or fossil horizons in clay/limestone successions.
Fig. 6. Box-counting analysis for RANDOM model with (A) relation between ®lled boxes N(e ) and box sizes e and (B) scaling behavior in power-law plot of fractal dimension d versus box sizes e.
In this paper the box-counting method is used to demonstrate the distribution of arti®cial events in the form of binary data, which exceed distinctive thresholds in the models used. For ease of understanding and recalculation, the range between Xmax and the mean plus one standard deviation is used to create `events'. The box-counting analysis supports detection of a single (unifractal) process and distribution pattern. The method is less useful for the analysis of distribution of events in time series, for which the in¯uence of more than a single process for the creation of a distinctive pattern (e.g. fossil horizons) is known in advance. Mathematical fractals show perfect self-similarity, although in reality, the geometry of fractals is never perfectly self-similar at all scales. Real life fractals are considered self-similar in a statistical sense. The idea of box-counting dimension is that the stratigraphic section of interest is covered with boxes of a ®xed size, and the number of boxes with data (part of the attractor) is recorded against box size; the process being repeated for dierent size boxes, and Ne Aeÿd where N = the number of boxes of a certain size (side length e ) which have point data occurrences in them; e=the length of box side; d = box counting dimension.
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
1015
Fig. 7. Box-counting analysis for time series models with relations between ®lled boxes N(e ) and box size e.
In this study, the number of event occurrences N(e ) in intervals (boxes) of 19 dierent sizes e has been calculated for all ®ve model datas. The box sizes range from 1 (original interval size) to the whole data set. Usually, the relationship between N(e ) and e is characterized by three stages: 1. Saturation occurs when events exist in all boxes (intervals) and an increase of the box size (which is inversely proportional to a decrease in the number of boxes) provides the same linear increase in the number of events. 2. Nonlinear functionality occurs when the number of boxes that contain events increases slower than the number of boxes according to a functional dependence. This scaling interval may be characterized by fractal scale invariant growth. 3. Depopulation occurs when the number of events increases much slower than the number of boxes or becomes stagnant. Results of box-counting analysis for the model RANDOM are shown in Fig. 6. The following three parameters are involved in the characterization of the scaling behavior of the data. 1. The relationship between N(e ) and box size e provides the number of `®lled' boxes. 2. The relationship of the ratio N(e )/Ng(e ) versus e gives the probability of getting a `®lled' box at a dis-
tinctive interval size e where Ng(e ) represents total number of boxes of size e. If N(e )/Ng(e )=1 there is complete saturation with events occurring in all boxes at least once. 3. The slope of D log(N(e )/log(e ) versus e gives the box-size range with power-law behavior; it means that scale-invariance is obtained for this box range. The interval used is relatively small in this example with box sizes e=4.5 . . . 8.6 and with a box-counting dimension d=ÿ0.7, selected to characterize stochastic random time series. Whereas the determination of the box size of saturation is exact, the determination of the boundary between depopulation and nonlinear function is subject to determination by the researcher. Fig. 7 shows power-law plots and growth probability of N(e ) versus e for events greater than mean plus one standard deviation in the time series model data. The nonlinear model FRAC39 can be detected only in box sizes between 3 and 6 with d=ÿ0.7 indicating a weak power-law relationship for these event-like data less than one scaling interval. The models SINUS, MODSR and UNCONF are characterized by stepwise increase of the probability function N(e )/Ng(e ), where each step may represent one of the cycle periodicities in the model data. The saturation limits of models SINUS and UNCONF are for 250 and 50 data, re-
1016
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Fig. 8. Multifractal spectra f(a ) versus a for ®ve continuous and event-like time series models, respectively.
spectively, representing major periodicities in the time series. Power-law relationships are suitable for describing the cycle-dominated models SINUS and UNCONF for distinctive box size ranges with d=ÿ0.93 and d=ÿ0.82, respectively (Fig.7). The cause of such scaling behavior is not mathematically described here and the reason why a periodic function can produce scale invariance is not fully understood.
5. Multifractal analysis In contrast to the fractal box-counting technique used for the study of binary data sets, multifractals analyze self-similar patterns using measures of variables such as those provided by continuous time series. Originally, multifractal analysis was introduced to describe growth probability distribution in kinetic aggregation processes where the growth probability for the perimeter sites can be regarded as a measure associated with each site (Amitrano et al., 1986). In this study, the sum of X(i ) was scaled to 1 by setting X(i )=X(i )/SX. The measures (data X(i )) in these intervals were treated as probabilities p of occurrence. Two possible multifractal patterns may occur in discretized continuous sedimentary time series: 1. One-dimensional features of spatially intertwined fractals coexisting in the same time series, such as
Cantor bars with dierent sizes (masses) (cf. Feder, 1988, pp. 64±65). 2. Nonstationarities of fractal (or multifractal) patterns in a time series caused by nonstationary nonlinear processes (e.g. volcanic or storm events). The following method can be used to determine a multifractal spectrum and ascertain if the data are best described by a fractal or multifractal distribution (Feder, 1988): An estimate of the amount of speci®c variable (e.g. logging value) per sample interval for a number of boxes that cover the time series of interest is obtained. All samples X(i ) within a box (interval) are used to obtain an average value for all boxes used. The measure for X(i ) per box can be set equal to this average value multiplied by e. Given a function wq(e ) for any real number q and box size e: wq e
N e X mqi e: i1
The measure mi(e ) satis®es a multifractal model if wq(e ) forms a power law relationship with e for any value of q, or wq(e )0 e t(q ) where t(q ) represents an auxiliary function. The HoÈlder exponent a satis®es a=d{t(q )}/dq, and can be estimated by using the Legendre-transformation (Feder, 1988) using a ®ve-point window. This procedure allows the multifractal spectrum f(a ) to be plotted according to f(a )=qaÿt(q ). Because we are
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
1017
interested in the scaling properties of the multifractal, it is desirable to know how wq scales with respect to e for a continuous range of q values with wq eAe qÿ1D q (Gould and Tobochnik, 1990). D(q ) is the generalized dimension of order q or dimension function introduced by Grassberger and Procaccia (1983). It can be estimated as D(q )=t(q )/(1ÿq ), q $ 1 (Feder, 1988). An unifractal distribution exists, if D(q ) versus q gives a straight line. Feder (1988) discussed that amin=ln p0/ln e0 and amax=ln p1/e1, with maximum probability or mass p0=Xmax/SX and minimum probability p0=Xmin/SX for associated intervals l0 and l1, respectively. For sedimentary successions, e1 and e0 are given by the thickness of layers with speci®c measurements such as grain size, not provided by the previous computer models in this paper. Relatively small HoÈlder exponents a support large values of X(i ), which tend to be clustered, whereas high a supports small values of X(i ) dispersed over the time series. For q = 0, all measurements become Xq(i )=1, indicating equal probability for all data in each scaling interval e. Then the multifractal spectrum reduces to fmax=d where d is the boxcounting dimension. For each calculation, six scaling intervals (segments) I for q=ÿ21, . . .,+21 were used, with I1=N/2, and Ii + 1=Ii/2, until I6=N/26, and N representing the size of the time series; e.g. 1000 values were used in model SINUS. Therefore, the number of segments ranged from 2 to 64. The distributions of measurements in very small boxes were not taken into account; Xmin and Xmax do not apply to the original data range, but represent averaged maxima and minima for I6. The multifractal spectra f(a ) versus a for all models are shown in Fig. 8. In general, the multifractal spectrum of the `upper threshold' time series considers more HoÈlder exponents a than the complete and continuous time series model. In each model fmax=1, which is a result of the saturation of the values for small box size e, as previously shown in Fig. 6. The spectrum f(a ) for model FRAC39 is close to a single spike characteristic of a unifractal distribution with fmax=1. The same applies to model RANDOM. Both models provide a nearly uniform distribution for each class of measurements (see Fig. 3). These distributions imply averaging of Xmin and Xmax at small e, resulting in small variances in the singularities a of the measurements. Relatively high f(a ) for small a is seen for model FRAC39 (upper threshold) whereas the distribution of f(a )/a is almost symmetrical in the model RANDOM (upper threshold). This dierence is caused by the relatively large number of large measurements X(i ) in model
Fig. 9. Dimension function plots D(q ) versus q for (A) three continuous and (B) ®ve event-like model data time series models.
FRAC39 in comparison with model RANDOM in which the measurements are uniformly distributed (see Fig. 3). In contrast, the models UNCONF, MODSR and SINUS show broad multifractal spectra for complete as well as `upper-threshold' data. There is no comparable multifractal spectrum for model SINUS because its completely periodic distribution does not provide a slope for a(q )=ÿdt(q )/dq. In model MODSR (upper part) at both ends its multifractal spectrum f(a ) is small (0.05) and the growth probability for the measurements with increasing box numbers I is small as well. The dimension function D(q )/q plot for the model UNCONF shows a large value of D(q ) for small q. This means that measurements X(i ) from just above the threshold value are widely dispersed, whereas high values, such as the single event at X(300)=11 results in a small dimension D(q ) (Fig. 9). The model RANDOM shows constant dimension both in the complete time series as well as in the `upper-threshold' time series. This is a result of the randomness of the process underlying this model. The model FRAC39 also shows approximately equal values
1018
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Fig. 11. Distribution plot of gamma-ray data of Kirchrode 1/ 91 well.
Fig. 10. Correlation dimension analysis for (A) gamma-ray data of Kirchrode 1/91 well (Upper Albian of North Germany); gamma-ray readings in API; (B) sliding windows correlation dimension; note the small variation of correlation dimension at embedding dimensions between 3 and 4; (C) D correlation dimension marks stability of underlying 3- to 4dimensional system; x marks most obvious sedimentary discontinuity (fossiliferous mud-turbidite horizon) with high D correlation dimension.
for all D(q ), which may be a result of time-averaging within the box sizes used for analysis. Computer programs in C for UNIX (TIMEMF.C for multifractal analysis, TIMESWCD.C for sliding window correlation dimension analysis) or TurboPascal for PC (FRACZEIT.PAS for box-counting analysis) are available on the IAMG server (ftp.iamg.org). Input and output ®les are generally ASCII-text ®les. Choice of input parameters such as window size, data range, box size range and box numbers is controlled interactively. 6. Case study: gamma-ray log, Kirchrode 1/91 (Cretaceous, North Germany) A gamma-ray log consisting of 2275 continuous measurements of marine marlstones was used for this case study (Fig. 10, Table 1). The section covers about 2.5 Ma representing almost the entire Late Albian (Cretaceous) in the North German Basin (BCCP Group, 1994). The cycle pattern in time and frequency was previously studied by wavelet analysis showing well-pronounced 12 m cyclicity in the upper part (40±165 m) and 18 m cyclicity in the lower part (178±216 m) of the section (Prokoph and Barthelmes, 1996). Mud-turbidites at the base and top of the section probably mark short-time hiatusses
within the corresponding intervals (28±8 m, 235± 210 m). The sliding window correlation dimension (200 data per window=20 m) of the gamma-ray data varies between 3.2 at 180 m and 4.1 at 80 m, possibly indicating a complex control of sedimentation by as many as three to four dierent factors in the section (Fig. 10(B),(C)). The saturation of the correlation function is shown by its ®rst derivative D, which reaches its maximum of 0.3 at about 220 m, where an hiatus is most likely to occur (BCCP Group, 1994; Prokoph, 1994). This high ®rst derivative indicates that chaotic randomness is most likely to occur in the interval 210± 230 m. In contrast, the interval from about 40 m to 170 m has D correlation dimension of less than 0.15. The Milankovitch cyclicity found at depths from 160 to 40 m (Prokoph and Barthelmes, 1996) probably was the predominant control of the sedimentation. The data are nearly normally (Gaussian) distributed with most data concentrated in a class around the average (70±80 API) as shown in Fig. 11. Both data distribution and correlation dimension are similar to those of model MODSR suggesting that nonlinear changes of sedimentation rate were at least temporarily superimposed by Milankovitch cyclicity. The box-counting technique shows a power-law relationship in the scaling interval e from 5 to 50 (=0.5± 5 m) with a box-counting dimension of d=ÿ0.5807 and saturation limit N(e )/Ng(e )=1 at box size 300 (30 m) (Fig. 12). However, this estimate can not be extrapolated to narrow data intervals. For example, 429 ®lled boxes predicted for box size 1 (=10 cm) exceed the 373 events counted by about 13% (Fig. 12(A)). This discrepancy is a result of depopulation. The multifractal spectrum of the gamma-ray data shows a spike-like curve of f(a )/a with f(a ) between 0.7 and 1 (Fig. 12(C)). Therefore, all classes of measurements X(i ) exhibit almost the same scaling
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
1019
Fig. 12. (A), (B) results of box-counting analysis, and (C) multifractal spectra of gamma-ray data of Kirchrode 1/91 well; r mark continuous data, mark time series data below mean-standard deviation, r mark time series data above mean+standard deviation.
behavior. On the other hand, the event-like data from the upper and lower threshold shows relatively high f(a ) > 0.2 for small a < 0.8. This may be a result of power-law relationships involving high gamma-ray values, probably exceeding 100 API. In contrast, gamma-ray values just above the threshold (about 85 API) are dispersed in the series and provide a > 2. 7. Concluding remarks Sliding window correlation dimension provides a method for the detection of temporary changes of the number of process controlling parameters through time (depth) and, therefore, allows the detection of sedimentary discontinuities such as facies changes, events and unconformities better than spectral and wavelet-analysis. Multifractal analysis is shown to be a tool for discrimination between unifractal, random and multiple time-related data patterns and for interpretation of the scaling behavior of classes of measurements in continuous time series. However, a parameterization of processes that controlled the data distribution shown in the multifractal spectrum and the dimension function plot is probably only possible for very high (Xmax) and very low values (Xmin) in the time series.
Box-counting techniques permit the discrimination of saturation and depopulation limits (sample intervals) of event-like data and the calculation of the scaling behavior between the corresponding two box size (interval length) limits. It can also give the box-counting dimension for power-law functionality between event clustering (`®lled' boxes) and the sampling interval size of observation (box sizes), results similar to those provided by lacunarity analysis (Plotnick et al., 1996). In this way it is possible to predict the maximum interval between events and probably the cluster characteristics of these events both in sedimentary sections and lateral distribution. Therefore, the ideas of Plotnick (1986) about the fractal distribution of stratigraphic hiatuses and Sadler's (1981) prediction of power-law relationship between sedimentation rate and observation time may also be related to scaling intervals between depopulation and saturation. Log±log regression with errors of less than 0.2% do not permit the determination of signi®cant scale invariance and subsequent prediction of sedimentary events for dierent sampling intervals. The time-event distribution patterns are likely to be expressed in their spatial distribution (i.e. sedimentary basin) as discussed by Yamaji (1988). Fig. 13 shows a model of the importance of saturation and depopulation limit as well as nonlinear func-
1020
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Fig. 13. Two-dimensional basin model for recognition of events due to dierent sedimentation rates and scaling (observation) intervals. For further explanation see text.
tionality in the distribution of discrete events of the same type (e.g. bentonite layers, fossil horizons, sequence boundaries) due to observation limits (measurement interval). For example, a sedimentary succession with four unevenly distributed volcanic ash layers (bentonites) in time and depth covers a time interval of 1 million years. The maximum resolution of a certain measurement method (e.g. gamma-ray logging) permits a resolution of 1 m-intervals. Therefore, in case the succession deposited during the 1 million years interval is only 1 or 2 m thick, every measurement detects a bentontite, i.e. the measurement is above the saturation limit, and the probability of detection of a bentonite layer is equal to 1. With increasing non-bentonite background sedimentation (e.g. marlstones) towards the basin centre, the number of measurements without bentonite signals (e.g. positive gamma-ray value) decreases and, ®nally, measurements represent only a single bentontite layer each. The depopulation limit is reached and the probability of ®nding more bentonite layers approaches 0 for increasing sedimentation rate. For analysis of nonlinear and nonperiodic continuous sedimentary time/depth series the following algorithms are suggested: 1. Sliding window correlation dimension calculation to
detect temporal discontinuities and the approximate number of underlying processes that probably controlled the dierent populations of the data groups (e.g. sedimentary pattern such as lamination and turbidite layers). 2. Truncation of time series at distinct discontinuities. 3. Multifractal analysis of segment cuts from the original time series by using threshold values to distinguish between unifractal, multifractal and nonfractal (linear, periodic) distribution. Estimation of growth probabilities at amin, a0 and amax. 4. Box-counting analysis of event-like data from the cut time series to estimate depopulation and saturation limits as well as nonlinear relationships in the restricted scale range. Together with trend and wavelet analysis, these techniques are tools to improve the knowledge about the temporary and spatial distribution of sedimentary patterns and processes. Computer programs of time-series analysis for box-counting analysis (FRACZEIT.PAS) in TurboPascal for PC, multifractal analysis (TIMEMF.C) and sliding window correlation dimension analysis (TIMESWCD.C) in C for UNIX are available on the Editor's Home Page (http://www.iamg.org/CGEditor/ index.htm).
A. Prokoph / Computers & Geosciences 25 (1999) 1009±1021
Acknowledgements I am grateful for discussions and reviewing dierent versions of this article by Walter Schwarzacher, Karl Stattegger, Bradley Sim and Frits Agterberg. The German Research Council (DFG) supports my research in Grant Pr 545/1-1/2. References Amitrano, C., Coniglio, A., di Liberto, F., 1986. Growth probability distribution in kinetic aggregation process. Physical Review Letters 57 (8), 1016±1019. Bak, P., Tang, C., Wiesenfeld, K., 1988. Self-organized criticality. Physical Review A38, 364±374. Baker, G.L., Gollub, J.P., 1990. Chaotic Dynamics: an Introduction. Cambridge University Press, Cambridge, p. 182. BCCP Group, 1994. The Upper Albian of northern Germany: results from the Kirchrode 1/91 borehole, Boreal Cretaceous Cycles Project (BCCP). Zentralblatt fuÈr Geologie und PalaÈontologie 1 (7/8), 809±822. Brescia, M., D'Argenio, B., Ferreri, V., Longo, G., Pelosi, N., Rampone, S., Tagliaferri, R., 1996. Neural net aided detection of astronomical periodicities in geological records. Earth & Planetary Science Letters 139, 33±45. Davis, J.C., 1986. Statistics and Data Analysis in Geology. John Wiley & Sons, New York, p. 646. Essex, C., 1991. Correlation dimension and data sample size. In: Schertzer, D., Lovejoy, S. (Eds.), Non-linear Variablility in Geophysics, Scaling and Fractals. Kluwer, Dordrecht, pp. 93±98. Feder, J., 1988. Fractals. Plenum Press, New York, p. 283. Fowler, A.D., Roach, D.E., 1993. Dimensionality analysis of time series data: nonlinear methods. Computers & Geosciences 19 (1), 41±52. Gould, H., Tobochnik, J., 1990. More on fractal and chaos: multifractals. Computers in Physics 3/4, 202±207. Grassberger, P., Procaccia, I., 1983. Measurement of the strangeness of strange attractors. Physica 9D, 189±208. Ivanov, S., 1996. Variability of sedimentary sequence: numerical modeling of the deposition-erosion process. Geologische Rundschau 85, 12±18. Kowalewski, M., 1996. Time-averaging, overcompleteness, and the geological record. Journal of Geology 104, 317± 326. Mandelbrot, B., 1983. The Fractal Geometry of Nature. W. H. Freeman, Freeman, p. 468. May, R.M., 1976. Simple mathematical models with very complicated dynamics. Nature 261 (5560), 459±467. Mudelsee, M., Stattegger, K., 1994. Application of the
1021
Grassberger±Procaccia algorithm to the d18O Record from ODP Site 659: selected methodical aspects. In: Kruhl, J.H. (Ed.), Fractals and Dynamic Systems in Geoscience. Springer, Berlin, pp. 399±421. Nicolis, C., Nicolis, G., 1984. Is there a climate attractor? Nature 311 (5986), 529±532. Nicolis, G., Nicolis, C., 1991. Nonlinear dynamic systems in the geosciences. In: Fransen, E.K., Watney, W.L., Kendall, C.G.StC, Ross, W. (Eds.), Sedimentary Modeling: Computer Simulations and Methods for Improved Parameter De®nition, vol. 233. Kansas Geological Survey Bulletin, pp. 33±42. Plotnick, R.E., 1986. A fractal model for the distribution of stratigraphic hiatuses. Journal of Geology 94 (6), 885±890. Plotnick, R.E., Gardner, R.H., Hargrove, W.W., Prestegaard, K., Perlmutter, M., 1996. Lacunarity analysis: a general technique for the analysis of spatial pattern. Physical Review E53 (5), 1±8. Prokoph, A., 1994. Zyklische Sedimentation im Oberalb des Norddeutschen Beckens. Tuebinger Geowissenschaftliche Arbeiten A19, 163. Prokoph, A. 1997. New mathematical tools and their applications in quantitative stratigraphy: Detection of cycles, chaos and unconformities by wavelet-and fractal analysis. In: Proceedings of the IAMG'97 Conference, Barcelona, vol 2, pp. 507±512. Prokoph, A., Barthelmes, F., 1996. Detection of nonstationarities in geological time series: Wavelet transform of chaotic and cyclic sequences. Computers & Geosciences 22 (10), 1097±1108. Rothman, D.A., Grotzinger, J.P., Flemings, P., 1994. Scaling in turbidite deposition. Journal Sedimentary Research A64 (1), 59±67. Sadler, P.M., 1981. Sediment accumulation rates and the completeness of stratigraphic sections. Journal of Geology 89, 569±584. Sanchez, A., Vos, H., 1994. From solar variability to climate proxy data Ð some ideas about signal transformation and the interpretation of time series. Terra Nostra 1/94, 19±24. Schwarzacher, W., 1993. Cyclostratigraphy and Milankovitch Theory. In: Developments in Sedimentology, vol. 52. Elsevier, Amsterdam, p. 225. Walsh, J.J., Watterson, J., 1993. Fractal analysis of fracture patterns using the standard box-counting technique: valid and invalid methodologies. Journal of Structural Geology 15 (12), 1509±1512. Wiesenfeld, K., Moss, F., 1995. Stochastic resonance and the bene®ts of noise: from ice ages to cray®sh and SQUIDs. Nature 373, 33±36. Yamaji, A., 1988. A fractal model for the distribution of stratigraphic hiatuses: a discussion. Journal of Geology 96 (6), 101±102.