Analysis for temporal patterns in Geoscience

Lecture index:


Reading:

Chapt. 6 Data Through Time - Swan & Sandilands, 1995, Introduction to Geological Data Analysis, Blackwell Science. 213-264. Read 213-220 in detail. Just skim the rest.


Introductory thoughts

Humans have a fascination with history. Naturally we especially treasure human history as our own, but for context for that history geologic history is only surpassed by cosmic history, and geologic history has the added advantage of being more accessible. Through understanding history we come to understand process and guess at destiny. You may remember quotes along the lines of those who do not know history are doomed to repeat it. A basic question is as to the pattern or character of history. Does it repeat itself or in what manner does it repeat itself? This is captured in the discussion on time's arrow and time's cycle. Irreversible system evolution is time's arrow. The seasons epitomize time's cyle. The distinction between and the necessity for both is nicely explored by Stephen J. Gould in his book Time's Cycle, Time's Arrow, and is encapsulated in the qupte below.

"time's arrow and time's cycle is, if you will, a "great" dichotomy because each of its poles captures, by its essence, a theme so central to intellectual (and practical) life that Western people who hope to understand history must wrestle intimately with both - for time's arrow is the intelligibility of distinct and irreversible events, while time's cycle is the intelligibility of timeless order and lawlike structure. We must have both. " Stephen Gould, p. 15-16, 1987.

One doesn't often think of math as provided a metaphor. However, one mathematically one can represent time's arrow with a linear relationship with time as the x variable. Given a value of y there is only one unique value of t that can work. A solitary sin or cosine function on the other hand can represent time's cycle. For a give value of y there are many possible times that it could correspond to. You know where you are in the cycle, but not where you are in the history. Put together they provide both times arrow and times cycle. There is not need that it be a linear function. An exponential decay function also provides times arrow (it just isn't "straight"). For realism one could add some random noise, and these curves could look even more realistic. History can be thought of as some mix of the directional, the cyclic, and the unpredictable.

Types of temporal patterns in geology:

What are external forcing agents that could induce cyclicity in the geologic record:

What geo-systems might have cyclic behavior driven by internal dynamics:

The eruptive history of Mt. St. Helens. From USGS site: http://vulcan.wr.usgs.gov/Volcanoes/MSH/Publications/MSHPPF/MSH_past_present_future.html

Eruptive history for the Cascades. From USGS site, PROJECT: Geologic Mapping in Southern Washington Cascades - http://vulcan.wr.usgs.gov/Projects/Mapping/SoWA_Cascades/summary.html

Eruption history for Mauna Loa from USGS site: http://pubs.usgs.gov/gip/hawaii/page14.html

A complicated approximation to the truth and how to capture it: Much of the above deals with simplified end-member types of behaviors. In reality, as mentioned above one can expect mixtures. Cyclic sea level changes driven by glacial-eustatic coupling and orbital forcing are likely superimposed on long term tectonic changes, with a possible random noise component, due to other factors. An equation that might describe such behavior could look something like:

y = k + m * t**p + RND(seeded on t)+ a * sin ( w*t) + b *cos ( u*t),

where y is sea level, and k, m, a, b, p, u and w are constants, and t represents time. The first two terms would represent a linear evolution if p=1 and a nonlinear powerlay relationship if t is not, and the last would represent a cyclic component with the constants a and b, and w and u, controlling the amplitudes and wavelengths respectively. By adding multiple sin and cos functions we can introduce interacting cycle components and some very interesting behavior. The relative wavelengths and whether the cycles are in phase or out of phase are important considerations in the behavior. As suggested one could also build an exponential component with ease.

This lab will introduce you a bit into analysis of time series data for temporal patterns. As a perusal of the latter part of your reading suggests this can be a very complicated process and it can be. However, initial entry is not difficult and if you take the time to think about the analyses they make sense. In addition to something called Markov chain analysis, you will also explore some simple plots that can identify patterns in time-series event data.


Vostok ice core data

If you drill at the right place on a glacier, as you go down you are sampling older and older snow that has been turned into ice. Looking at a range of impurities, at trapped air bubbles, and at the isotopic character of the H2O, there is a tremendous historical record. As snow falls on to a glacier, seasonal differences often result in a layering in snow/ice texture so that annual increments can be recognized. There are also other ways to provide time constraints. This has resulted in a lot of ice cores being drilled from glaciers around the world, and the longest and most spectacular of these may be the Vostok ice core data.

Antarctic Vostok ice core data source:

Graph of T changes recorded in Vostok core as measured by hydrogen isotopes. Look at the pattern qualitatively. What do you notice about it? Here the pattern is strongly developed - obvious. Often it is not.

Other ice core data from Antartica showing some of the variables measured with depth. From NOAA site Eight glacial cycles from an Antarctic ice core.Image source: http://www.ncdc.noaa.gov/paleo/pubs/epica2004/epica2004.html


Below are two more in-depth examples of how one might analyze for temporal patterns: the runs test and Markov chain analysis. The focus in this course used to be on the Markov chain analysis, but has switched to the runs test simply because it seems the runs test may be more useful. Data to apply the Markov chain analysis on is harder to come by and there is more complexity/ambiguity when it is applied to stratigraphic data in the transformation of the data to a sequence of states. So you may want to familiarize yourself with both, but in terms of completing the assignment please focus on the runs test.


Using the runs test.

This exercise is based on p. 261-264 in the Swan and Sandilands textbook we have referred to before. We will introduce it in the context of fining-upwards or coarsening upwards trends in sedimentary sequences, but it can definitely be generalized to other uses.

In sedimentology you will learn that one of the hallmark traits of a deltaic progradation sequencing is a coarsening upwards trend and that of a fluvial channel lateral migration is a fining upwards trend (in sediment grain size). In turbidite fans similar coarsening upwards and fining upwards trend are expected and analyzed for, and these can reflect sediment source and alluvial fan dynamics. Similar trends can be expected and sought in bedding thickness. However, the signal can be noisy, and given our human tendency to assign pattern when there is none, or the challenge of recognizing a pattern in a very noisy system, one can use a variety of statistical analyses and tests. A robust and simple test is the runs test.

What is a run? It is simply a consistent sequence. In this case we will consider 2 simple events given passage through a temporal (stratigraphic) sequence; either a decrease or an increase with time. In a perfect coarsening upwards sequence one would see one long positive run where grain size (or bedding thickness) would consistently increase. If there was a random progression of events this would be very unlikely. However, just by chance a random sequence would have some runs. The longer the random sequence, the more runs, and the longer the longest run would be expected. A perfect alternation of increases and decreases would also be a clear non-random sequence. Why might such a sequence exist? What if it is not a perfect coarsening upwards or a fining upwards sequence? Can we test for how unlikely the sequence was produced randomly? This is what the runs tests does.

The Z statistic for the runs tests basically measures the difference between the observed number of runs and the expected number of runs for a random distribution. This means that the null hypothesis is that the data resulted from randomness, and the alternate hypothesis is that it did not.

Z = ( (observed # of runs) – (expected number of runs) ) / (expected variance)

The expected number of runs is a function of the length of your record, and the number of observations of each state (the number of positives = n1, and the number of negatives = n2) and is:

expected number of runs= 1+ (2 * n1 * n2)/ (n1+n2)

The expected variance is also a function of n1 and n2 and is equal to the following:

= ((((2 * n1 * n2) * ((2 * n1 * n2)– n1 – n2)) / ((n1 + n2)^2 * (n1 + n2 – 1)))^.5

How do you interpret the z statistic itself? Using the standard threshold of 5% (your willingness to conclude it is non-random when it is not – or on the flip side of the coin, a 95% level of confidence), results in Z values of 1.96 and -1.96. If your observed number of runs is much higher than expected (a high positive Z value) then it suggests some type of alternating or oscillating system. There could be good sedimentologic reasons you would expect this. If your number or runs is lower than expected (longer, more continuous trends), a more negative Z value would result. So the steps you can use to apply the runs test in Excel are as follows.

  1. Get your stratigraphic grain size, bed thickness or other data, in to a column (label it of course).
  2. Use the logical IF statement in Excel to assign an increase from one cell to the next cell a state designator (+1 is handy), and a decrease from one cell to the other a different state designator (-1 works well).
  3. Count the number of observed runs. This can either be done manually (tedious for a long sequence, and possibly prone to error), or it can be done in Excel. You can use the IF function so that if the two adjacent values in the column where you have your state values (your +1 and -1 values) is the same you are still in a run (and a value of 0 can be assigned) versus if they are different and the end of a run has been reached (and a value of 1 can be assigned).
  4. Using two COUNTIF statements you have used before you can have Excel compute n1 and n2 from the values in the column that has the values calculated in step 2.
  5. Using the SUM function and the values in the column calculated in step 3 will give the observed number of runs.
  6. Use the formulas to compute the expected standard error and the expected number of runs (using n1 and n2 values).
  7. You can now compute the Z statistic, and then start thinking about what it means, and why it might the way it is.

You should clearly label all your columns and cells were calculations of n1, n2, etc.. are made. This is typically only the start of an analysis, unless it tests as random, and then you might be done. If you conclude it is non-random you can think of other steps that might be taken to look for patterns.


Markov Chain analysis of stratigraphic successions

Markov chain analysis is used to detect repetitive patterns in how sedimentary units are stacked. One of the powerful aspects of it is that it can be used on nominal data such as mineral associations in a layered mafic complex, where rhythmic patterns are suspected.

Another place that Markov chain analysis has been used is in the analysis for cyclothems in sedimentary sequences. So what is a cyclothem?

Cyclothems: A cyclothem is a model for a pattern of repeated facies migrations that leaves a distinctive reptititve pattern in the vertical stacking of lithologies (remember Walther's law). They are particularly well known from the Carboniferous, and also from the Cretaceous (Brenner & McHargue, 1988). They are developed for and can bear the name of a region and/or geologic period. The cause for the repetition can be internal factors (such as lobe switching in a deltaic environment) or external forcing agents such as eustatic changes. Better known examples include the Illinois cyclothem, which is thought to be a result of delta shifting, and the Pennsylvanian Kansas cylothem (see figure to left) which is attributed to glacial transgression-regression cycles on a shallow carbonate-shale shelf. Cyclothems related to faulting have also been explored. You might think of other geologic situations that might produce a repetitive vertical sequence or rocks.


This strat column shows an idealized section of a Pennsylvanian cyclothem and is based on work mainly in Kansas, but similar rocks outcrop in SE Nebraska also. Diagram from http://capp.water.usgs.gov/gwa/ch_d/D-summary.html.

Remember that humans are biased to seeing patterns even where they do not exist. This is why Markov chain and other types of even more rigorous statistical analysis are needed to establish that there is actually some sort of repetitive pattern and what the pattern is. Field work and a qualitative/intuitive assessment may lead one to develop a cyclothem model. Markov analysis can be used as a test of the hypothesis.

You can think of a sequence of lithologies as a 'chain' of letters and ask the question as to whether there is a pattern in this chain. Naturally, the longer your chain the more confidence you will have in the results of the analysis ("Please, sir, can I have some more data."). Markov Chain analysis is based on the simple question of whether a given lithology is independent or not on the lithology below. The greater the dependence the more likely that transition from one bed to another is part of a pattern of behavior. Technically this is known as a first-order Markov analysis. You can also ask to what degree a lithology is dependent on a sequence of lithologies below. If you look at the two underlying lithologies then this would be a second-order Markov analysis. The simple comparison is to what would be expected randomly (equivalent to independence). If 25% of our lithologies as represented in the chain, are limestone, then 25 percent of the time we would expect limestone to lie on top of any other lithology (including limestone) if stacking were random. Using a matrix we can look at all the observed transitions and see how they compare to a random stacking pattern. We will look only at first order Markov Chains.

Your reading does a very good job of describing the basis of first order Markov chain analysis. The key is constructing transition frequency and probability matrixes. You can model your analysis directly after the example in the book.

We are always interested in finding new places to use our skills and tools to profit. Can you think of other non-geologic applications where this type of analysis could provide insight.


An example - peritidal carbonate sequences

Source of information and data: Journal of Sedimentary Research; November 1996; v. 66; no. 6; p. 1065-1078 © 1996 SEPM Society for Sedimentary Geology Facies successions in peritidal carbonate sequences, Bruce H. Wilkinson, Nathaniel W. Diedrich, and Carl N. Drummond

One might expect shallow carbonate platforms to be very sensitive to sea level changes. The table below is part of such a lithologic sequence, measured bed by bed, and the assignment of each lithology to one of four categories, and then the transition that is seen. For example, for the second lithology, TK-BEDDED MUD_SILT, which is arbitrarily assigned the letter D, it is follow above by FINE RIBBON ROCK which has been assigned an A, and so the transition is that A follows D. The table below that is the Markov matrix analysis. Your reading goes through the steps in detail.

FINE RIBBON ROCK A
TK-BEDDED MUD-SILT D AD
COARSE SAND C DC
MASSIVE SILT D CD
MEDIUM SAND C DC
FINE RIBBON ROCK A CA
CRYPTALGAL LAMINITE A AA
MASSIVE SILT D AD
TK-BEDDED MUD-SILT D DD
FLAT PEBB CGL B DB
TK-BEDDED MUD-SILT D BD
SHALE D DD
TK-BEDDED MUD-SILT D DD
MASSIVE MUD D DD
FINE SAND C DC
COARSE SAND C CC
CRSE RIBBON ROCK A CA
COARSE SAND C AC
FINE RIBBON ROCK A CA
COARSE SAND C AC
FLAT PEBB CGL B CB
CRYPTALGAL LAMINITE A BA
FINE SAND C AC
MEDIUM SAND C CC
MASSIVE SILT D CD
EXPOSURE BRECCIA D
MASSIVE MUD D D
MEDIUM SAND C DC
CRSE RIBBON ROCK A CA
CRYPTALGAL LAMINITE B AB
TK-BEDDED MUD-SILT D BD
FINE RIBBON ROCK A DA
TK-BEDDED MUD-SILT D AD
FINE RIBBON ROCK A DA
MEDIUM SAND C AC
CRSE RIBBON ROCK A CA

A B C D OBSERVED TRANSITION MATRIX
A 16 10 30 58 114
B 11 0 5 14 30
C 26 5 7 64 102
D 60 15 60 139 274
520
A B C D OBSERVED TRANSITION P MATRIX
A 0.1403 0.0877 0.2631 0.5087
B 0.3666 0 0.1666 0.4666
C 0.2549 0.0490 0.0686 0.6274
D 0.2189 0.0547 0.2189 0.5072
A B C D RANDOM TRANSITION P MATRIX
A 0.2192 0.0576 0.1961 0.5269
B 0.2192 0.0576 0.1961 0.5269
C 0.2192 0.0576 0.1961 0.5269
D 0.2192 0.057 0.1961 0.5269
A B C D DIFFERENCE P MATRIX
A -0.0788 0.0300 0.0670 -0.0181
B 0.14743 -0.0576 -0.0294 -0.0602
C 0.03567 -0.0086 -0.1275 0.10052
D -0.0002 -0.0029 0.0228 -0.0196
A B C D RANDOM COUNT MATRIX
A 24.9923 6.5769 22.3615 60.0692
B 6.57692 1.7307 5.88461 15.8076
C 22.3615 5.8846 20.0076 53.7461
D 60.0692 15.8076 53.7461 144.376
A B C D DIFFERENCE COUNT MATRIX
A -8.9923 3.42307 7.6384 -2.0692
B 4.4230 -1.7307 -0.8846 -1.8076
C 3.6384 -0.8846 -13.0076 10.2538
D -0.0692 -0.8076 6.2538 -5.37692
A B C D CHI-SQUARED MATRIX
A 3.2354 1.7816 2.6092 0.071279 7.6975
B 2.9745 1.7307 0.1329 0.206719 5.0450
C 0.5920 0.1329 8.4567 1.956258 11.138
D 0.00007 0.0412 0.7276 0.200248 0.9692
24.849 CHI-SQUARED STATISTIC
16.92 REJECTION LEVEL .O5 AND N.F. = 9

The conclusion is that this sequence is nonrandom. Note that the Chi-square test says nothing about the validity of any proposed cyclothem model - only that the sequence as analyzed is unlikely to be random, opening up the possibility that there is a rhythm present. The researcher then can use the more common transitions to propose a repeated sequence, and then a geologic model to try and explain why that sequence is common.


Exploratory plots for time series analysis

There are a number of simple plots that can be relatively quickly done for 'event' data that can detect patterns and guide further analysis. We will explore four such plots using geyser eruption data from Old Faithful. This geyser is famous for its regularity, and so it is optimal for learning purposes (and indeed this data set is used in some statistic courses - see reference below). Yet, the temporal behavior is richer than you might think, and it has changed through the decades. The source of the data is: http://www.stat.cmu.edu/~larry/all-of-statistics/=data/faithful.dat . One can get annual data for eruptions from 200-2011 from The Geyser Observation and Study Association site - http://www.geyserstudy.org/geyser.aspx?pGeyserNo=OLDFAITHFUL if you want to do more sophisticated analysis.

Azzalini, A. and Bowman, A. W. (1990). A look at some data on the Old Faithful geyser. Applied Statistics 39, 357-365.

Old Faithful eruption in 2012 to the right, and people waiting for the next eruption below (by the time it occurred it was a wall of people 10 or 15 deep around a good bit of the geyser).

A simple first plot to do is a histogram of the intervals between eruptions and a histogram of the Old Faithful data referenced above is below. n is 272 in this data set.

The plot suggest that there is at least two modes, one positioned at roughly 52-54 minutes and one at roughly 76-80 minutes. One might be tempted to conclude that there are more peaks, but there is insufficient data in this data set to really tell if this is the case.

A second plot one can do is a serial correlation plot, where for each event the time between it and the preceding event is x and the time between it and the following event is y. The resulting plot for the Old Faithful data is below.

In these plots clusters indicate a periodicity to events. Note here that we have three clusters. Compare the position of the clusters to the peaks in the histogram - do they match. Which potential cluster is missing or under represented in this plot, and what does this tell you?

Another type of plot is the cumulative frequency plot where x is time, and y is the number of events up to that point in the history. The resulting plot for the Old Faithful data is below.

This plot shows a very nicely defined line, which can be interpreted as constancy of activity over the history of events - Old Faithful is particularly well named according to this graph. Breaks or steps in the line would be useful for detecting some distinct change in activity, perhaps some threshold change in the system. There is some small scale fluctuation in this line (think of the two modes) that is hidden by the size of the symbols used to plot each point.

The finaly type of plot we will consider is an empirical Survivor Plot where x is eruption interval and y is the proportion of the intervals of that duration or greater. With some thought it makes sense that a negatively sloping line should result starting at the left with 100% of the points to the left of the minimum value and 0% of the points right of the maximum value in the data set.

This type of plot is perhaps more difficult to interpret thatn the previous three. The plot will have an exponential form if events occur randomly in time. To see departures more easily a log-log plot will turn an exponential form into a straight line and then departures from that line can be seen. Here there are clearly some steps and the data is clearly non-random. Think of how a distinct mode would appear in a plot like this.

What might be driving such regular geyser eruption behavior? What might be driving the more irregular portion of the behavior (look at how well the clusters are defined) or not?

Image from the USGS showing an idealized view of the plumbing for a hot-spring - geyser system. Source - http://pubs.usgs.gov/fs/2002/fs101-02/ .


Exercise 13

1) Analyzing for runs in turbidite sequences.

Turbidites often build out as parts of submarine fan sequences, and it could make sense that there are runs that are the result of fan progradation, or turbidite source evolution (perhaps driven by changes in sea level).

Photograph of part of a turbidite sequence from up in Alaska, Bay of Pillars Formation. Source: Reconnaissance Geologic Map of the Duncan Canal-Zarembo Island Area, Southeastern Alaska By Susan M. Karl, Peter J. Haeussler, and Anne McCafferty, U.S. Geological Survey Open-File Report 99-168 Version 1.0 1999 http://pubs.usgs.gov/of/1999/of99-168/of99-168-pamphlet.html

In a 1997 Master's thesis by C. Chen, in the appendices, are tables of bedding thickness and grain size data that is very amenable to this type of analysis. Indeed, the title of the thesis is Statistical Analysis of Turbidite Cycles in Submarine Fan Successions. Chose one of the sections in Appendix 1 and analze either the bed thickness or the base grain size for runs. The steps for doing this analysis are described above.

Your final product should include the following:

Probably the best approach is to hand it in as your Excel file, with text boxes providing information and interpretation where needed.

1(Optional) - Analyzing for cyclothems in the Billefjorden basin of Spitsbergen.

You have (or will be given) a copy of a portion of a stratigraphic fence diagram from a Carboniferous basin in Spitsbergen (taken from Johannessen & Steel, 1992, Mid-Carboniferous extension and rift-infill sequences in the Billefjorden Trough, Svalbard). Authors claim cyclicity. These are sediments deposited next to a fault in a half graben. Cycles could be related to episodes of fault movement and sedimentary readjustment or to orbital or glacial/eustatic forcing or some type of combination. You will test for it. You should basically replicate the Markov Chain analysis as described in the first part of your reading for this diagram.

The lithologies you will be working with are a) limestone/dolomite, b) evaporites, c) fluvial sandstone and shallow marine sandstone and conglomerates (lumped together), d) and alluvial fanglomerates. You should thus have a 4 by 4 matrix and 16 possible transitions (similar to what is in your reading).

Your final product should include:

View of Trikolorfjellet where the Carboniferous basin fill is well exposed on Spitsbergen.

2) Graphical analysis of a series of events - The Aso eruption case history:

Much of this is taken from Davis, 1986, Statistics and Data Analysis in Geology, 2nd edition. Here they consider a long history of eruptions for the Volcano Aso on Japan (Kyushu). The data is given below.

1229, 1239, 1240, 1265, 1269, 1270, 1272, 1273, 1274, 1281, 1286, 1305, 1324, 1331, 1335, 1340, 1346, 1369, 1375, 1376, 1377, 1387, 1388, 1434, 1438, 1473, 1485, 1505, 1506, 1522, 1533, 1542, 1558, 1562, 1563, 1564, 1576, 1582, 1583, 1584, 1587, 1598, 1611, 1612, 1613, 1620, 1631, 1637, 1649, 1668, 1675, 1683, 1691, 1708, 1709, 1765, 1772, 1780, 1804, 1806, 1814, 1815, 1826, 1827, 1828, 1829, 1830, 1854, 1872, 1874, 1884, 1894, 1897, 1906, 1916, 1920, 1927, 1928, 1929, 1931, 1932, 1933, 1934, 1935, 1938, 1949, 1950, 1951, 1953, 1954, 1955, 1956, 1957, 1958, 1962

Data in form easily loaded into Excel.

You may also use the Vostok ice core data if you want to explore that data set. I also encourage you to look for other eruption histories to analyze. Just make sure you have enough data points to do something with. A lot of eruption history data can be found at http://www.volcano.si.edu/world/ . You also need to provide the source for your data.

First, look at this data and see if you detect any pattern as you peruse it.

The following plots are simple to do in Excel and provide a first look for any pattern in the series.

A cumulative frequency plot: This plots the total number of events that occured at or before time t versus t. The slope between any two points on the plot gives the average number of events per unit time. A straight line would indicate constancy of occurrence. Distinct breaks in the plot would suggest a change in behavior throughout the history. Sorting your data will help make creating this plot fairly easily.

Histogram plot: Basically create a bar graph showing the number of events within time bins of constant interval. If a well developed trend or pattern is developed it may show up here. Remember that histograms are sensitive to the interval position and width.

Empirical survivor plot: In this plot x is the length of a time interval between events, and y is the proportion of such event intervals in the record longer than x (you might want to reread that). This plot will have an exponential form if events occur randomly in time.

Serial correlation (first-order correlation) plot: this is plot of x as the interval before a given data point versus y the interval after a given data point. A plot with a lot of scatter and a higher concentration of points near the origin is typical of a random series of events. Think about what a plot would look like if intervals were spaced about every 10 then 20 year cycles. Two clusters symmetrically disposed should be evident on the plot.

Hand in these 4 plots (labeled) and describe your results in 250 words or less. Can you think of any problems with the data set that might obscure a pattern?

Note by the way that even if eruptions happen randomly or close to randomly through time, potentially they can still be predicted using precursors and a monitoring system.


References:

Copyright by Harmon D. Maher Jr.. This material may be used for non-profit educational purposes if proper attribution is given. Otherwise please contact Harmon D. Maher Jr.