Chapter 13 Some problems of three-dimensional stratigraphic analysis

Chapter 13

Some problems of three-dimensional stratigraphic analysis

13.1 INTRODUCTION

A complete stratigraphical analysis must involve the study of three-dimensional rock bodies within a reference system of time, as was emphasized in the preceding chapter. Given ideal mapping conditions, it ought t o be possible to obtain the shape and distribution of the lithostratigraphic units, at least. A lithological parameter, let us say r, can then be expressed as a function of the three Cartesian coordinates x , y , and z and it may be written as r = r (x,y,z). Matheron (1965)calls the variable r a “regionalized” variable. This is intended to be a neutral term which indicates that the value of r at any given point, is the outcome of natural phenomena (Matheron, 1967). The regionalized variable can be generated by some stochastic process but it does not necessarily have to be so. Functions of the type r ( x , y , z ) are usually very complicated and erratic in geology and they are consequently very difficult t o describe by straightforward mathematical expressions. For this reason, one must resort to sampling techniques and these latter invariably involve a random element when the function r is determined. Quite apart from the “sampling error”, there is a distinct possibility that the actual shape of an observed rock body is the outcome of a random process. The model concept can be employed once again, in order to simplify three-dimensional stratigraphic analysis and it is clear that two types of models are required. Firstly, there must be a “process model” which deals with the dynamic system of the environmental factors which will determine eventually the type of sediment which is produced. Secondly, it may be necessary to have in addition, some simplification of the observed configuration of the lithostratigraphic units and this can be achieved by substituting a model for the geometrical shape of the rock body. It follows that the “process model” can be either some mathematical expression, or it can be any set of rules which portray the actual process. This model is used to generate the “geometrical model” which itself may be a simplified version of a rock distribution that has been observed in the field. The term “process model” was introduced by Sloss (1962)and was contrasted with the term “response model” which refers to the resulting configuration of sediments. A response model might be thought to be almost the same as a geometrical model but the geometrical model is always independent of time, whereas the sedimenta-

350

PROBLEMS OF THREE-DIMENSIONAL STRATIGRAPHIC ANALYSIS

ry response model is a feature that develops with time. Thus time-dependent response models are feasible. Indeed, the process model is time-dependent by definition and it often interacts with the response model and feed-back relationships can exist between the two (Krumbein and Sloss, 1963). Three-dimensional stratigraphic analysis can be simplified by thinking of a rock body as if it consisted of a series of two-dimensional sections. This treatment can be applied also t o the different types of models which have just been discussed. For example, the geometrical model can be thought of as being a series of superimposed maps and if each map represents an isochronous surface, this can be made the basis for a dynamic process model. In a similar way, maps of this kind can be used for plotting any sedimentary parameter and a time-dependent response model can be constructed by choosing the vertical spacings between the maps as a time axis. Let us assume that the parameter which is plotted on a system of maps like this, represents the thickness z that was deposited during unit time. Each map then represents a record of the function z t = z ( x , y ) , . Geographical variation can be ignored in dealing with single sections, making it possible t o concentrate on the simpler function z = z ( t ) , in the earlier models. When it is realized that the variable z is a function of both time and areal variation, the analysis becomes considerably more complicated and one of the major problems in the three-dimensional approach is to assess the extent to which the two functional dependencies can be separated. 13.2 THE AREAL VARIATION OF SEDIMENTS

That conditions of sedimentation and consequently the rocks which are formed in an environment change with their geographical distribution, is self-evident. The facies variation is very often systematic in its nature and it is related directly to such geographical factors as the distance from land, the water depth, or a combination of both. Indeed, palaeogeographic maps of an environment can only be reconstructed because these relationships are known. Deterministic models can be developed for some facies variations. For instance, the thickness of an individual bed which has been produced by a unidirectional current pattern is known t o decrease exponentially with the distance from its source. Similar dependencies exist for other lithological parameters like grain size or roundness, for example. Theoretical models which attempt t o explain these observations, can be based on the simple differential equation that was introduced in Section 3.5. The diffusion process that was introduced in the same section, provides an alternative way of obtaining a facies distribution model. One may assume, for example, that the source of the sediment is a topographically high area and that the sediment will spread from the source in accordance with [3.14],which gives the thickness of the sediment for each time increment.

AREAL VARIATION OF SEDIMENTS

351

Sedimentation processes incorporate various stochastic elements, by their very nature and almost any process will be associated with a certain number of random fluctuations, regardless of whether it involves sedimentation from transporting currents or settling or organic precipitation. If these variations occur independently during deposition in different localities, they will generate non-correlated noise in the areal distribution of the sediments. There must be many additional processes which occur during sedimentation and which can lead to correlated random terms. For example, a stream bed may be covered with pebbles which are scattered at random on a sandy surface. The pebbles will cause disturbances in the flow pattern which lead to the well known scour marks being made behind the pebbles and possibly, to increased deposition in front of them. Thus the random disturbances will affect adjacent areas and introduce correlated random deviations into the areal distribution of the sediments. Very similar effects can be expected to arise from wash-outs, sand bars and ripple marks. Indeed, the latter can be regarded as a case of oscillating areal distribution, where sediment transport leads t o geological cycles in a horizontal direction. The formation of correlated random disturbances is by no means restricted to the sediments that are produced from unidirectional currents alone. The ecologist, for example, is quite familiar with the occurrence of high-density patches in the distribution of plants and animals. The clustering together of various plant and animal communities may have a number of causes, and the mode of reproduction is a common one. The distribution of organisms can have a profound effect on sedimentation and can lead to areal variation which, once again, is highly autocorrelated. Unfortunately, no detailed studies have been made yet to investigate the areal correlation structure of sediments or to relate this to the primary processes of deposition. It is obvious that the practical difficulties which arise in making the observations (Section 12.1),have been a major obstacle but if they could be overcome by using suitable exposures, it would afford a very promising field of research. The various types of areal variation which have been discussed in this section, suggest that one could construct models for the regional facies variations which are very similar t o the models that are used for describing the time series properties of single sections. Indeed, Agterberg (1967, 1970) suggests the following general model: data = trend + signal + noise

[13~1

The term “trend” refers to deterministic functions such as polynomials and it is assumed that this part of the areal distribution is caused by one of the deterministic processes which were discussed at the beginning of this section. The signal is assumed t o be stationary in the wide sense and is characterized best by its autocorrelation. The noise is an uncorrelated random disturbance.

352

PROBLEMS OF THREE- DIMENSIONAL STRATIGRAPHIC ANALYSIS

In accepting the Agterberg model ([13.1] ), one might hope that the estimation of areal variation could proceed by using methods which are similar to those used in the estimation of the lithological variation in time-series analysis. Unfortunately, a number of difficulties arise which have not been met with in previous problems and an attempt will be made to show how these can be overcome, at least partially. 13.3 THE ESTIMATION OF AREAL TREND

The study of areal trends can be considered first of all. This is approached traditionally by fitting polynomials to the lithological data, using the methods of regression analysis. Thus the method is similar to some of the trendfitting procedures which have been used previously in time-series analysis, but with one important difference. In time-series analysis, trends like the linear or the exponential trend can be explained by simple physical processes, making it possible to assess the effects which the removal of such trends would have on the data. Although similar, simple processes like the exponential decrease of grain size with increasing distance from the source area, can operate in areal variation, the geographical configuration of depositional basins can lead to very complicated distribution patterns which cannot be approximated by low-order polynomials. For example, if a basin’s shore line is taken to be the source area for the sediments, then this is an essential boundary condition which will determine the distribution of the sediments. Clearly, the sedimentation patterns c m be approximated by low-order polynomials only if the basin shapes are geometrically simple ones like circular or ellipsoidal basins. One has only two options in dealing with the complex facies patterns which occur more usually. One can either abandon the idea that a simple relationship exists between the geographical coordinates and the lithological variation, or one can restrict the investigation to a sub-area which is small enough to be uninfluenced by the boundary conditions of the basin shape. For example, one might study the lithological variation along a series of sections or narrow strips which are taken at right angles to the shore line, in order to find the relationship between sedimentation and the distance from the land. In this case, one has t o assume that the individual sections are independent of each other and that lateral relationships can be ignored. This sitpation can be compared with a study of sediment transport which restricts itself to laminar flow where the particles are allowed to move in definite planes only and where turbulence terms are omitted. Models like these can give valuable information about current velocities and sediment transport but they are unavoidably incomplete. In a similar way, low-order polynomial trends can be used t o bring out the essential facies variation but they cannot always by interpreted in a physical sense.

ESTIMATION OF AREAL TREND

353

An alternative approach to the analysis of areal trend was initiated by Krige (1966) who designed a method which can produce a trend surface map which is not based necessarily on the assumption of a simple functional relationship existing between the lithological variation and its geographical distribution. Indeed, it is implicit in Krige’s approach that the method could be applied equally successfully to data that are generated by a random or by a deterministic process. The method was developed primarily as an aid in the evaluation of ore reserves. The amounts of rare minerals like gold fluctuate widely within a given rock body so that one has to introduce a certain amount of smoothing in order to make predictions or interpolations. This is achieved by making use of a weighted moving average, using methods which are similar to the ones which are used in time-series analysis. The area is subdivided into a number of blocks, measuring 30 m2 for example. An estimate of the mean mineral content yo is made at a point which is centred in one of these blocks. The estimate is based on all the sampling data X1 within this block and those for a number of surrounding blocks ( X i , i = 1,2,...n). A linear prediction equation is used and one can write for the estimate of Xo : -

2 , = f :aixi+6,x

[13.2]

i= 1

The weights ai can be determined in areas for which accurate estimates of X O exist. This done, one can apply [13.2] to areas for which only the Xi’s are is a correction for the regional mean and the known. The term a, constant regression is given by:

x

x

m ...

6,

=1-

c /-. i= 1

ai

[13.3]

An alternative method of finding the weight factors ai can be used when the two-dimensional autocorrelation of an area is known (Agterberg, 1970). Consider a number of points, p l , p 2 ....pn which surround a point p o and for which a measured lithological variable X i with an assumed zero mean is known. The value at p o can be estimated by:

go=

c aixi n

i= 1

The least-square estimates for the coefficients ai lead to the solution:

[ 13.41

3 54

PROBLEMS OF THREE-DIMENSIONAL STRATIGRAPHIC ANALYSIS -1

r2.1

1

r2,3

. [ 13.51

which is the two-dimensional equivalent of the Yule-Walker relation ([7.31] ) as it is used in time analysis. If the autocorrelation function is isotropic in the plane, then a value rij will depend only on the distance tijbetween two data points. Such distances can be determined from the geographical coordinates x,y and:

The distances tijare calculated between each of the surrounding points p l , p2...p,, and rij is determined from the known autocorrelation function rr. The coefficients ai can then be determined from [13.5] and they may be used for smoothing in the manner which is indicated by [ 13.21 . The term “kriging”, has been used by some geologists to describe this type of smoothing operation, although it is not entirely clear why such an exotic terminology should be thought necessary. Equation [ 13.21 constitutes a linear filter which can be applied to the data and as such, it is a very flexible method. It should prove very useful in sedimentological problems where the analytical expressions of the polynomial trend surface are too rigid. Nevertheless, its main application will continue t o be in the study of areal random processes and these will be discussed in the next section. 13.4 STATIONARY RANDOM PROCESSES IN THE AREAL DISTRIBUTION

Agterberg’s general model which has just been introduced, refers t o the component which is generated by a stationary random process as the “signal”. This does not imply necessarily that this component is the most important one in a regional study. Agterberg (1970) makes it quite clear that the three components depend on the sampling scheme that is being used and if a small area is sampled in detail, the “signal” may very well become part of the trend. The situation is rather similar to the one that is found in time-series analysis and perhaps, it is advantageous t o replace the term “signal” by the more neutral term of “oscillating component”.

STATIONARY RANDOM PROCESSES IN AREAL DISTRIBUTION

355

The analogy between geographical variation and time-series models is not a complete one, however. Time is a unidirectional parameter by definition and this means that the state of a variable f ( t o )only can be determined by L 2 . .and . ) that the future values of f ( t l , t 2 . . . ) can have past events f(L1, no physical influence on the value at t o . This is not true necessarily, for the regional variables and for example, f ( x o ) can be determined by f( ...xP2, X-,, x l r x2 ...) that is, it can be determined by points which are in the past as well as in the future, if distance is regarded as being the equivalent of time. There are exceptions, of course. In the study of sediment transport problems, where f ( x ) may represent a variable which is measured along the non-reversible down-stream direction, the distance function becomes equivalent t o a time function. Whittle (1954) calls such relationships “unilateral”, whereas the normal areal variation which may be observed along a straight line is “bilateral”. Agterberg (1970) gave a summary of Whittle’s approach to multidimensional autocorrelation and the following discussion is based essentially on this. One may consider first the simplest bilateral correlation scheme which can operate along a straight line. Assume that a point p o lies half way between points p1 and p 2 . The values X1 and X 2 are known and X o is to be predicted. It seems reasonable to give the values X , and X 2 equal weights, ul = u2 and the prediction for X o is therefore:

i,= &,X, + 62x2

[ 13.71

Using [13.5], one obtains:

C:] [r2 :1‘ [:1

[13.8]

=

where rl is the correlation coefficient for the distance p o p l = p o p 2 and r2 is the correlation coefficient for the distance p o p 2 = p o p l . It follows that: [ 13.91

For a series of discrete data, a value of xk which is half way between X k and X k + can be obtained from the model:

Xk

=

a(&-1

+X,+,)

+ fk

1

[ 13.101

It can be proved that the term E~ is an uncorrelated random value only if the series X , ( h = 1, 2,...) is a first-order Markov chain. Equation [13.10] is thus

3 56

PROBLEMS OF THREE-DIMENSIONAL STRATIGRAPHIC ANALYSIS

the equivalent expression to the unilateral time-series process:

x, = cx,-,

+ Ek

[ 13.111

and since it has the Markov property, one can write its autocorrelation function as:

rr = r!’

[ 13.121

Equation [13.9] can be simplified by using [13.12], to give: [ 13.131

and it may be noted that the equivalent expression for the unilateral process is C = r l , as can be seen from [7.27]. A continuous autocorrelation function can be formulated by letting the sample interval tend towards zero. From [ 13.101, one can derive a stochastic differential equation:

[13.14] where a 2 = (1 - 2a)/a and this leads to the autocorrelation function:

The equivalent expression for the unilateral process has been noted already and it is: [ 13.161

r7 = e-P171

and the two expressionsin [13.15] and r13.161 are clearly not the same. The Markov property of the bilateral process can be defined as:

whereas the unilateral process has been defined as: P(XtIX7;7

< t)=

lim P(X,IXeAt) At-rO

[ 13.181

STATIONARY RANDOM PROCESSES IN AREAL DISTRIBUTION

357

The differential equation which was given in [13.14] can be generalized for the p-dimensional space and Whittle (1954) gives a general solution for the autocorrelation functions which arise from such processes. He considers particular solutions for the two-dimensional as well as for the three-dimensional process. The autocorrelation function of the two-dimensional process turns out to be:

[13.19] where Yl ( a t ) is a modified Hankel function, meaning a Bessel function of the second kind. The autocorrelation for the three-dimensional process gives the familiar exponential:

113.201 Whittle’s results are of considerable importance in the interpretation of geological data. They suggest that if one makes the transition from the discrete regional processes to the continuous random processes, the exponential autocorrelation function (which is a condition for the Markov property) can be found only in three-dimensional processes. Empirical autocorrelation functions which were derived from geological data, have been sampled along straight line transects. For example, samples were taken at regular intervals and analysed for ore content in studies by Krige (1962) and Agterberg (1965, 1967). Such data were then treated exactly like a time series and a correlogram was calculated. It was found that the autocorrelation function of all the examples could be approximated by an exponential expression that is comparable with [13.20]. One must conclude in applying Whittle’s results, that the random processes which were responsible for the mineral distribution must have operated in three dimensions. This is a reasonable assumption to make in the particular case of the mineral distribution. It should be noted that if the autocorrelation structure is exponential in the three-dimensional space, then it will also be exponential in any lower dimension, such as along the direction of transects which were used in these studies, for example. However, if one were to restrict the physical process model to a lower dimension then it could not provide an explanation of the observed distributions, according to Whittle’s results. This restriction must be kept in mind when one proceeds to analyse areal variation by a series of independent transects, as was recommended in Section 13.3. A one-dimensional model could be quite inadequate to explain the observed correlation structure of a lithological variable. Unfortunately, no quantitative work on the areal correlation structure of either sedimentological or stratigraphical data is available yet. This is in spite of it being likely that many of the features which are observed on bedding planes could have originated from processes which were essentially two-

358

PROBLEMS OF THREE-DIMENSIONAL STRATIGRAPHIC ANALYSIS

dimensional. Of the available data, those which come from ecological observations on vegetation patchworks are the most relevant. For example, Pielou (1964) showed that if an area which is populated by a single species that occurs in clusters, is classified into two states involving plants being either present or absent, then a transect which is taken in any direction and records the states that are encountered can be represented by a two-state Markov chain. It is very likely that the process which is responsible for the distribution of the plants is essentially two-dimensional and this appears t o contradict Whittle’s results. However, Switzer (1965) showed that there exists at least one discrete random process in the plane, which has the property that the states which are encountered along a straight-line transect have the properties of a Markov chain. Switzer’s process will be referred to again in the next section. The estimation of the areal autocorrelation is straightforward, providing that the data points are situated on a regular sampling grid. Under such circumstances, correlograms can be calculated along straight-line transects for any direction in the plane. Naturally, if the areal distribution of the correlation structure is isotropic, a single correlogram is all that needs t o be calculated. The correlogram can be represented graphically by contours in an x,y coordinate system, the origin of which represents the zero lag position. The isotropic case, for example, can be shown by circular contours which normally decline in value from the centre of the circle, which has the maximum correlation coefficient of 1.Most anisotropic correlation structures will yield ellipsoidal contour lines. Difficulties arise in the calculation of correlation functions when the data are obtained from irregularly spaced sample points. The problem of estimating the autocorrelation under such circumstances, was discussed by Agterberg (1970) who showed that at least an approximation of the autocorrelation structure can be obtained.

13.5 THE LITHOLOGICAL VARIATION IN STRATIGRAPHIC PROFILES

The preceding sections dealt with problems of three-dimensional rock bodies which were visualized as being a sequence of superimposed maps. In addition to this, a lithosome can be studied in a series of vertical slices which are the familiar stratigraphical profiles, of course. If one were to deal only with the geometrical distribution of rock types, there would be no difference between these two methods. However, the map represents only a short moment in the environmental history which is represented by a single time horizon. The profile, on the other hand, records geographical variation along its horizontal axis and its vertical axis measures the thickness of the sediments and is also a time axis, at least by implication.

LITHOLOGICAL VARIATION IN STRATIGRAPHIC PROFILES

359

Whether the vertical axis should be divided into equal thickness or equal time intervals, is a matter of choice. The profiles which result from these two alternatives are different and they cannot be converted into each other unless both of the scales are given. From a purely practical point of view, it is found that thickness profiles with added time lines are much more easily read than the time profiles with added thickness information. However, the latter are useful in the study of facies variation, where the thickness of the strata can be ignored to some extent and time-scale profiles have been used in simulation studies as well. The type of information that is needed for constructing a thickness profile with added time lines, may be seen from Fig. 13.1. The two graphs on the left-hand side of the figure give firstly, bed by bed, the length of time for which a certain sedimentation process lasted and secondly, the rate with which this process proceeded in each bed. The thickness of the beds is found by multiplying the sedimentation rate by the time taken for the formation of that particular bed. To keep the example as simple as possible, the sedimentation rates in Fig. 13.1 are chosen so that each bed attains a constant thickness. Next, the time lines are constructed by multiplying the index of the time line by the appropriate sedimentation rate for each locality. It may be noted that the second time line enters bed 2 twice, for instance, and it is the sedimentation rates for this particular bed which are relevant to finding the position of the time line. The assumption of absolutely continuous sedimentation has been made in constructing the thickness profile that is shown in the upper right-hand side of Fig. 13.1. This is a most unlikely situation, as has been discussed at some length in Chapter 2, because petrographic studies suggest that a break of deposition occurs after a bed has been formed and

i

+w :- I

I

w-

&I-----

----

Fig. 13.1. The construction of a profile with time lines. Broken lines on the right show time stages.

360

PROBLEMS OF THREE-DIMENSIONAL STRATIGRAPHIC ANALYSIS

discontinuous sedimentation appears to be a more natural model. This more realistic situation is shown in the lower right-hand side of Fig. 13.1. Time lines have been constructed in this figure under the assumption that at least five time units elapse after the completion of bed 1, before the second bed is deposited. Under these circumstances, the time lines merge with the bedding plane which becomes an isochronous surface by this process, eventually. It is clear from the method of constructing the thickness profile and its time lines, using the data that are given on the left-hand side of Fig. 13.1, that this process is not reversible. As in all the previous examples, the thickness of the sediments is given by the product of the sedimentation rate and the period of time over which the sedimentation occurred. It is natural that in constructing a process model, one is interested mainly in the sedimentation rates and in the duration of the sedimentation process, both variables being dependent on the geographical position. Any model attempting to explain an observed profile must make some assumptions about the sedimentation rates, the length of time for which a particular sedimentation process lasted and the regional distribution of these two factors. One can think of various geological situations where either the sedimentation rates or the duration of the sedimentation processes are the main interest. In high-energy fluviatile environments, for example, one might think of sedimentation as consisting of a patchwork of areas where deposition occurs over different and variable periods of time, but with approximately constant rates for each facies. In basinal sedimentation, on the other hand, the time-dependent factors may be constant over large areas but the sedimentation rates may vary regionally either at random, or following some systematic trend. Different models must take such geological peculiarities into account. The following specific example has been chosen to illustrate some other points that arise from the design of models which generate geological profiles. Assume that some random-walk-like mechanism exists t o produce exponentially distributed sedimentation steps. This situation was discussed in Section 6.3 and may arise when sedimentation is steady, apart from irregular and exponentially distributed interruptions which are followed by a certain time of non-deposition. In this situation, the thickness of the sediment that is deposited in unit time can be expressed by a Poisson process with, let us say, density p . Let it be assumed further that at regular time intervals, sedimentation changes occur so that some marker horizons such as bedding planes, are formed. The incidents producing the bedding planes are also distributed exponentially with density X(X < p ) . This sedimentation model was considered in Sections 6.2 and 11.4 and it was seen from [11.11] that the thicknesses of the beds which are measured in a single section are exponentially distributed. The model may be expanded into two dimensions by considering a string of localities that are arranged along a straight line and by assuming that each locality generates its own section. If sedimentation proceeds quite independently in each locality, then both the sedimentation

--

LITHOLOGICAL VARIATION IN STRATIGRAPHIC PROFILES A

361

B

Fig. 13.2. Two artificial profiles based on independent random sedimentation ( A ) and synchronous random sedimentation (B).

steps occurring with density p and the incidents which lead to bedding planes occurring with density X, will happen everywhere at different times so that they have only their distribution laws in common. This would result in the type of profile which is shown in Fig. 13.2A.The thickness distribution of beds which are measured in individual sections, is the distribution of [11.11] and the thickness distribution which is obtained by measuring the same bed in different localities, is identical to the vertical distribution of single sections. This independence of different localities is very unlikely from a geological point of view since it implies that there are no regional dependencies within such an environment. A much more realistic model can be obtained by assuming that the incidents producing the bedding planes occur almost simultaneously over the distance covered by the profile. For example, one may assume that the time events which subdivide the section are caused by the influx of clastic material into a predominantly carbonate environment and that this spreading is relatively rapid, in comparison with the vertical sedimentation rates. Modifying the model in this way does not alter the bed-thickness distribution as it is measured in individual sections, but it does change the thickness distribution of an individual bed that is measured in different localities. One and the same bed now represents a time interval which is constant in every locality and therefore its thickness must be determined by the compound Poisson distribution that was introduced in [6.26] and [6.27].The effect of bedding planes occurring simultaneously, may be seen in the short simulated section which is illustrated in Fig. 13.2B.The difference between the first model (A) and the second model (B) is obvious. When sedimentation is independent (model A), the beds cannot be correlated by matching the thicknesses of beds from different exposures but this can be done for model B. When model A operates, all the mean values which are measured for an identical bed in different localities will be the same. When model B operates, the averages which are derived from the different localities do permit the reconstruction of the function qT (see [11.8])and one can estimate at least, the relative time history of the sedimentation process. Use was made of this argument when discussing the example in Section 11.4.

362

PROBLEMS OF THREE-DIMENSIONAL STRATIGRAPHIC ANALYSIS

It has been pointed out that model A is unrealistic because of its independence of regional variation. However, the same is true for model B to a lesser extent. The latter assumes that the minor fluctuations which can be described by parameter p are completely independent of the geographical location, whereas the incidents causing bedding planes were taken to operate with deterministic certainty wherever they occurred throughout the region. According to what has been discussed already, it seems more likely that sedimentation processes are regionally autocorrelated on a minor scale and that incidents such as the influx of clastic sediments are stochastically disturbed in their regional pattern. Realistic models for generating geological profiles should combine the stochastic processes that are responsible for the vertical sequence with the stochastic processes that are responsible for the lateral variation but double stochastic models lead to difficulties, as was seen in the previous section. A simple example will illustrate this. The Markov-chain model has proved very popular in the description of vertical sections and so one might consider generating stratigraphic profiles by placing a great number of sections having the Markov property, next to each other. It can be seen immediately from Fig. 13.3 that such Markov chains can only exist next to each other if they are completely independent. The first section can be taken as an example. The state which occupies cell ( 2 , l ) should depend only on the content of the cell ( 1 , l ) . Similarly, the state which occupies cell (2,2) should depend only on the state of (1,2). Next, assume that the Markov property exists in the horizontal distribution of the states and that for example, the series (l,l),(1,2), (1,3) ... constitutes a Markov chain. This must lead to the conclusion that state (1,2) depends on ( 1 , l ) and that state (2,2) depends on both state (1,2) and state (2,1),which is a contradiction of the Markov property. Of course, this was brought out by Whittle’s results which were considered in the previous section, where it was mentioned also that the only simple two-dimensional Markov process that is known so far, is the Switzer process. Switzer (1965) showed that if an area is subdivided by a series of random lines to produce a Poisson field with density hdpde, where p and 8 are the coordinates which are used for defining the lines, then a number of convex cells result. If the states occupying the cells are chosen at random and independently, from a two-state system, then the resulting pattern that is shown in Fig. 13.4 has the property that any straight-line transect through this area constitutes a Markov chain. A model of this type is by no means as unrealistic as it may appear to be from the diagram which is given in Fig. 13.4, where the density of the random lines was kept relatively low. If the concept of sedimentary “quanta”, which was used primarily to describe stratigraphic time units, is expanded into two dimensions, then it might be possible to identify such cells as real sedimentological units but more work will be needed before this can be attempted. Although the Markov process provides a simple model for random variation, other schemes such as the higher-order autocorrelated processes, prove

ll,,

MODELS FOR SIMULATING GEOLOGICAL PROFILES

2 , l 2,2 1,l

1,2 1,3

17

,

363

I :.-I.a -

Fig. 13.3. The combination of a vertical and horizontal Markov chain. Fig. 13.4. Jtealization of the Switzer process.

to be important in the study of vertical variation and these must be combined with the lateral variation processes which themselves might be of similar complexity. Not only are there no data for such studies but also there is no complete theory for two- and three-dimensional random processes. It is not surprising that the models which have been developed for generating geological profiles are as yet very simple and have not yielded very exciting results. 13.6 MODELS FOR SIMULATING GEOLOGICAL PROFILES

A model which is particularly instructive in connection with the problems which have just been discussed, was developed by Krumbein (1968).The model is based on the principle of transgression and regression “cycles”, where it is assumed that a marker horizon such as a beach deposit, for example, moves backwards and forwards within the studied profile. The sedimentation process is studied by observing the position of the marker horizon at various fixed localities, such as the sections A , B, C ... in Fig. 13.5. One can say that the system is in state A at a given time, when the marker horizon is recorded in section A . It is assumed that sedimentation is steady and constant so that the thickness of a deposited bed can be equated directly with time. To give a full description of the strand-line movement, one must indicate whether a given state is approached during the trangressive or the regressive phase of the cycle. If it is assumed in Fig. 13.5 for example, that the sea is on the left-hand side and the land is on the right-hand side, then states A , B, ... are all part of the transgressive phase and states A‘, B’ ... are all part of the regressive phase. One can assume further that the system operates at regular time intervals and this enables one to write a transition

PROBLEMS OF THREE-DIMENSIONAL STRATIGRAPHIC ANALYSIS

364 A

B

C

C

0'-

A>-

Fig. 13.5. Transgression and regression cycles simplified after Krumbein (1968).

matrix for the position of the marker horizon. For example, if the transgression-regression cycle is entirely deterministic, meaning that if the strand line migrates from A t o B to C and returns to A , then the following cyclical matrix describes the process:

-

-

A

0

B O P=

c

0

1 0

0

0

0

0

0

1

0

0

0

0

1

0

0

C'O

0

0

0

1

0

B'O

0

0

0

0

1

A '-1

0

0

0

0

0-

[ 13.211

The deterministic model specifies that whenever one state is reached tb e marker horizon must move t o the next state in the following time incren, 9t. A stochastic element can be introduced by allowing a certain probabk P that the strand line remains in the same position, for a variable amount. time. In the deterministic model of [13.21], the full transgression must completed before a regression can occur and a further stochastic eleme. which was introduced by Krumbein, permits a reversal of the transgression cycle at any intermediate state. If only three sections are considered, as in the example above, then the reversals can occur not only at A and C but also at B . This is indicated in the diagram of the possible transitions which is shown on the right-hand side of Fig. 13.5. The process is Markovian and a

MODELS FOR SIMULATING GEOLOGICAL PROFILES

365

typical transition matrix could be written as:

P=

A

B

C

C‘

B‘

A’

A

0.02

0.98

0

0

0

0

B

0

0.02

0.78

0

0.20

0

C

0

0

0.02

0.98

0

0

C’ 0

0

0

0.02

0.98

0

B‘ 0

0.20

0

0

0.02

0.7E

A‘ 0.98

0

0

0

0

0.02 -

-

13.221

The example assumes that the rates of transgression and regression movement are equal but it would be a simple matter t o write the transition matrix in such a way that the transgression is more rapid than the regression. This would lead to more “natural looking” simulation results. The lateral and vertical variation patterns which result from the profiles that are generated by this model, are of interest. The first important observation is about the lateral distribution of the lithologies. Since the model is based on the concept of migrating parallel facies belts, the order of these lithologies cannot be disturbed and consequently, the sequence of rock types dong any time plane is entirely deterministic. For example, if one associates rock type 1 with marine sediments, rock type 2 with the littoral sands and rock type 3 with terrestrial sediments, then the sequence of the rock types always must be 1, 2, 3, in going from the left t o the right-hand side along any horizontal plane. It follows that the entire profile can be constructed if the rock sequence is known for any of the single sections. Section B may be considered as an example and if the system is in state A , then rock type 3 will be deposited in that section. The states can be equated with the rock types in the following way. A = A’ = 3, B = B‘ = 2, C = C‘ = 1. The matrix in [13.22] can be rearranged in order to bring out the rock sequence wliich is generated in section B (see [13.23] on the following page). This is still a Markov matrix, of course, but the 3 X 3 matrix referring to the succession of the lithologies is not Markovian and the stratigraphical section can be described by a Markov chain, only if it can be decided on lithological grounds, whether a sediment is formed during the transgressive or the regressive phase of a cycle. The stratigraphic sequence of undifferentiated lithologies can be modelled by a semi-Markov process and for this to be achieved, it is only necessary to find the distribution of the waiting times for each state. These may be obtained from 113.231. It can be seen that the

PROBLEMS OF THREE-DIMENSIONAL STRATIGRAPHIC ANALYSIS

366

1

2

3

C

C’

C

0.02

0.98

0

0

0

0

C‘

0

0.02

0

0.98

0

0

B

0.78

0

0.02

0.20

0

0

B’

0

0

0.20

0.02

0

0.78

A

0

0

0.98

0

0.02

0

0

0

0

0

0.98

0.02

1

s=2

3 A’

B

B’

A’

A

[13.23]

waiting time for state 1 for instance, is determined by the length of any chain which consists only of state C or C’. However, the probability of such an uninterrupted chain of C’s having a length of n , is given by the first passage probability of state B to B’, since these two states always terminate such chains. It follows from the arguments in Section 5.4, that one can write: [ 13.241

Equivalent expressions can be found for all nine submatrices of [13.23]and after the appropriate calculations, it should be possible to find a matrix P(n) which provides a complete model for any section that is selected in this system. An examination of Krumbein’s model shows that an observed Markov sequence of lithologies is incompatible with the concept of transgressions and regressions, even if the transgression pattern follows a Markov process. The limitations of this model are caused largely by the strictly deterministic relationships existing between the vertical and the lateral variations, which have been discussed already as “Walther’s law” in Section 12.6.If more flexible models are to be achieved, a search must be made for processes which not only permit random variations in the vertical direction of the profile but in its horizontal direction also. Jacod and Joathon (1971,1972) constructed some models which were intended to be the basis for simulation studies which were designed to investigate an actual stratigraphic problem as well as proving process models for some rather complex sedimentation environments. The basic assumptions in the Jacod-Joathon model are simple. Sedimentation is a non-decreasing function Q t ( x ) which depends on the time t and on the geographical location, x. The thickness of the sediment which is deposited during the time interval At

MODELS FOR SIMULATING GEOLOGICAL PROFILES

367

is proportional t o : . [ 13.251

Sedimentation is a function of the water depth W t and the amount of crustal sinking R t and it can be expressed as: [ 13.261

which is equivalent to [3.21]. The function f [ W t ( x ) ] indicates the dependence of sedimentation on the water depth and this could be determined by the amount of sediment dispersion, as mentioned in Section 3.7. It is assumed in addition, that the environmental factors Rt and Qt(x) are random functions. The tectonic component R t is represented by a random walk, whereby both the step size and the intervals between the steps are distributed exponentially with the parameter y. In addition to the discontinuous tectonic movement, there may be an additional constant subsidence which proceeds with speed p. The specification of the variable Qt(x) is rather more difficult. Not only is it dependent on time of course, but also on the geographical location x. In order t o generate such a variable, Jacod and Joathon (1971) employ a stratagem which is based on the following scheme. The two horizontal axes of a Cartesian system are used t o represent the geographical location while the vertical axis represents the time. The horizontal plane is scattered with Poisson-distributed points of a given density. Each point is in the centre of a disc-shaped body which is flat on its base and its upper surface is formed by a parabolic cap. Each time interval At is represented by a new map of Poisson-distributed bodies. When the scheme is viewed in a two-dimensional profile, the disc-shaped bodies appear in section as they are shown in Fig. 13.6 and the distribution of such cross-sections is again Poisson. Each of the randomly distributed bodies is used to define a function G i ( x ) which will determine the sedimentation Q t ( x ) . The following rules are used. Take a point where x = a , in Fig. 13.6, then the perpendicular line which is erected in point a, will represent the time history of this location. It is assumed that when the line does not pass through a disc, dQ,(a) = dt, but when it passes through disc number i, dQ,(a) = Gi(a). A typical realization of such a process is shown on the right-hand side of Fig. 13.6, where it is assumed that Qt carries out a jump which is proportional t o the intercept of the time path passing through any disc. The scheme which has been outlined in the preceding paragraph can be modified in several ways and a typical application of the Jacod-Joathon model is to simulate the lenticular sedimentation of a high-energy environment. For example, one can specify that the functions G i ( x ) are associated with the deposition of sand, whereas clay may be deposited in between the

PROBLEMS OF THREE-DIMENSIONAL STRATIGRAPHIC ANALYSIS

368

I

> X

O,P)

Fig.,,l3.6. The Jacod-Joathon model of sedimentation.

discrete bodies. If an environment like this is subjected to further irregular tectonic movement which affects the sedimentation rates according to [13.26], then one can simulate some very realistic profiles. Indeed, Jacod and Joathon (1972) have shown that by selecting a limited number of parameters and by comparing the simulation results with actual field measurements, the model can produce results which are in very good agreement with the observed structures. Simulated geological structures can be used in the study of practical applications like the gas or the ground-water flow in complex geological bodies and actual field measurements can be used to either confirm or improve the quality of the model. At the same time, it must be admitted that the work of Jacod and Joathon throws very little light on the genesis of random processes which operate regionally. It is assumed that the function Q ( x ) is constant during a short time interval A t and that it is determined only by the water depth of the basin at that moment. This state of average sedimentation is disturbed regionally by a number of randomly distributed centres which give rise to local lens formation. Sedimentation in localities which are close t o each other is not independent because each lens has a finite size but the model does not allow for the interaction of overlapping lenses. Indeed, it is assumed that the maximum diameter of the lenses is smaller than the distance between the sample points. Therefore the model does not give any information which might be used to investigate why and how lenses form and our ignorance of this process cannot be disguised by the arbitrary assumption that the centres of such disturbances are distributed in an independent and random manner. 13.7 THE ENVIRONMENTAL MODEL

The environmental model was introduced originally in Section 1.3 as a system of factors which are linked together by functional relationships. Perhaps it is wishful thinking to attempt to replace this whole complex of physical activity by a system of mathematical expressions which simulate the

THE ENVIRONMENTAL MODEL

369

environmental processes and which will generate a three-dimensional body of sediments, comparable to stratigraphic systems in nature. One brave attempt in this direction was made by Harbaugh (1966) who compiled a computer programme which can simulate the sedimentation processes occurring in a basin which contains various sediment-producing organisms and which receives clastic sediments from various directions. An attempt is made t o show how organic growth, the distribution of sediments and the geometry of the basin interact to provide various sedimentary bodies. Although this model is vastly simplified, it can simulate quite complicated rock distributions. The main lesson t o be learned from Harbaugh’s attempt is possibly twofold. On the one hand, it can be demonstrated that even quite simple models can produce some very realistic three-dimensional sediment bodies. On the other hand, the model highlights once again how little is known about the sedimentation process. Most of the functional relationships, such as the dependence of organic development on external factors and the dependence of sedimentation on the water depth or the distance from the land, as well as other relationships, are known incompletely. All the initial conditions such as the geometry of the basin, the tectonic framework and the influx of the external sediment have t o be guessed a t before the model can operate. If the model incorporates yet more random processes to make it more realistic, then the number of unknown parameters increases correspondingly. The question which arises is whether an attempt should be made to make the models ever more realistic and t o develop well-fitting empirical models, as discussed in Section 3.4,or whether one should concentrate on the theory of environmental processes. It appears that the latter approach can provide the more meaningful explanations, but it is hampered seriously by the complexity of the processes being investigated. However, theoretical models can deal with simplified situations or situations concerned with definite aspects of sedimentation. For example, Krumbein’s work which was considered in Section 13.4 showed that the relationship between lateral and vertical rock variation needs urgent investigation. Very little is known yet, about the random processes which affect the areal variation of sedimentological parameters and even less is known about the compatibility of various stochastic schemes in three dimensions. It is possible that quite simple theoretical models combined with detailed stratigraphical mapping and sedimentological research may help towards a better understanding of the space-time relationships in sediments. It is only after the basic problems of environmental interaction have been understood, that one can hope to approach complex operations like stratigraphic interpretation or stratigraphic correlation, from a quantitative point of view.