- Email: [email protected]
Tracer dynamics in ocean sediments and the deciphering of past climates
Mathl.
Compvt.
Modelling Vol. 21, No. 6, pp. 27-38, 1995 Copyright@1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0895-7177/95 $9.50 + 0.00
Pergamon
08957177(95)00021-6
Tracer Dynamics in Ocean Sediments the Deciphering of Past Climates
and
C. NICOLIS Institut
Royal
3, Avenue
Mktkorologique
Circulaire,
de Belgique
1180 Bruxelies,
Belgium
(Received and accepted April 1 L&X$) Abstract-An analysis of tracer dynamics in ocean sediments taking into account the disordered, fractal, character of geological materials is performed. The complexity of the medium is modeled by augmenting the equation of motion of the tracer through the addition of a stochastic forcing. In the one-dimensional version of the model, it is found both analytically and numerically that the dynamics is characterized by a large dispersion of arrival times at a given depth of the sediment as deduced from the first passage time probability distribution of the tracer. Moreover, the trajectories of two initially close tracer particles deviate at later times according to a power law, thereby raising the question of reliability of climatic trends deduced from geological data. Numerical simulations in a more realistic, two-dimensional setting show that this uncertainty becomes even more pronounced entailing the formation of a fractal precipitation front whose width increases with time. The implications of the results on the interpretation of geological records are discussed.
Keywords-Ractals,
Glaciations,
Predictability,
Random
walk
1. INTRODUCTION Much of our understanding of the earth’s past environmental and climatic conditions rests on our ability to decipher geological data, particularly those pertaining to ocean sediment formation. Now, reconstructing a time sequence of events from the spatial distribution of a representative tracer within the sediment involves, perforce, a number of highly nontrivial intermediate steps. Specifically, since as a rule only a limited number of points can be dated with a good accuracy, a model of sediment formation has to be adopted before the “space axis” along the sediment can be converted into a “time axis” along which past history begins to unfold. A classical approach to tracer distribution within sediments makes use of advection-diffusion models, sometimes augmented by the inclusion of adsorption or production-destruction terms. The important issue of the mixing of biological origin is usually addressed through a variety of phenomenological assumptions relative to the depth dependence of the diffusion coefficients involved in the theory [l-4]. Typically, all these models adopt a deterministic description limited to a one-dimensional view, in which attention is focused on the tracer distribution along the direction of advection. Our principal objective in the present work is to develop an enlarged description of tracer dynamics within a sediment, which incorporates the influence of stochastic effects. The principal motivation behind this idea is the fact that geological materials are complex media characterized This work is supported, in part, by the Commission Scientific, Technical and Cultural Affairs (O.S.T.C.).
of the European
Communities
and by the Federal Office for
Typeset Ia4 21:6-C
27
by dn/ls-QX
C. NICOLIS
28
by a pronounced disorder. Within such a medium, the movement of a tracer is expected to become, effectively, a random process since the instantaneous direction and speed are likely to vary continuously depending on the locally prevailing conditions. This should entail, in turn, that the position of the tracer after a sufficiently long lapse of time will be subjected to a considerable uncertainty. We want to characterize this uncertainty in a quantitative manner and to explore its consequences in the interpretation
of the geological record, especially the one pertaining
to
climatic evolution. In Section 2, we discuss in some detail the origin of the randomness of geological materials and formulate the basic equations governing the dynamics in such a material.
A one-dimensional
version of these equations is analyzed in Section 3. In particular, we evaluate the probability density of the position of the tracer as well as the first passage time distribution determining the times of arrival at a given depth. Both distributions a large uncertainty
display a considerable dispersion, entailing
on the position of the tracer, in agreement with the comments made above.
A complementary view of this uncertainty is explored in Section 4, where we study the mean deviation of the trajectories of two initially nearby tracers. As it turns out, this deviation increases steadily in time according to a power law. This illustrates further the problems inherent in the deciphering of climatic trends from the geological record. In Section 5, the analysis is extended to two dimensions. It is found that the presence of a second dimension gives rise to a new mechanism of variability in the form of a fractal diffusion front whose average width increases with time. The implications of the results are discussed in Section 6, in connection with the glacial-interglacial climatic transitions.
2. GENERAL
FORMULATION
Consider a tracer moving in a geological medium. Depending on the circumstances, it may represent a chemical substance dissolved in water, a macroscopic particle (like, for instance, debris of a certain biological species) containing a radioactive isotope, etc. In either case, owing to the disorder of the medium-including the one arising from the action of living organisms-we stipulate that, fundamentally, this movement is to be regarded as a random walk. Denoting by r the instantaneous position of the tracer and by t the time, we therefore write
dr - = v(r) + qg(r).F(t). dt Here v denotes the advection (drift) velocity; F a random force accounting for the action of the medium on the particle; g a coupling term; and q the strength of the coupling. The complexity of the medium can be further manifested through the fact that v and g may have a complicated r-dependence, which may even be viewed most adequately as a random function of r. We do not allow for an explicit time dependence of v and g as this is taken into account by the random force F(t). We shall make on F(t)
the simplest possible assumption of a (normalized)
Gaussian white
noise:
Pi(t)) = 0 (Fi(t)Fj(t'))
=
6;; b(t -
(2)
t’)
We emphasize that these assumptions can be relaxed, but in the sequel, we shall not be concerned with this generalization. It is well known that equations (1) and (2) are equivalent to a Fokker-Planck equation for the probability density P(r, t) of finding the tracer around point r at time t (see [5]):
at=-c& w(r)P(r, t) + $ C dP
i
2
ijk
d gik & dr,,
gjk
P(r, t)
(3)
Deciphering
of Past
Climates
29
Notice the formal similarity between equation (3) and the advection-diffusion equations describing the transport of material in sediments: Here the role of advection is played by the drift velocity V, whereas
the term g coupling
diffusion
whose magnitude
the tracer depends
movement
on the strength
to the random parameter
force gives rise to an effective
q.
Let us comment in somewhat more detail, on the connection between the complexity of the medium and the random character of the dynamics. It is well known that geological media present a very pronounced irregularity and fragmentation, which can hardly be viewed as a small, quantitative modification of some regular “reference” crystal lattice. To capture the essence of this irregularity, we assimilate the material to a fractal [6,7]. The specific choice of fractal is a major question whose answer depends on the nature of the material intervening in a given problem. Here we consider an abstract-but typical-model in which the medium is regarded as a triadic Cantor bar, as depicted in Figure 1 (see [S]). Specifically, it is assumed that all bar sizes between two extreme values r,in and T,~,, are present with a probability P(r) induced by the rules leading to the construction of the Cantor set.
-
m
m
-
I
II
II
II
II
II
II
II
II
II
II
II
II
II II II II
II II II II
II II II II
II II II II
II II II II
II1
IIII
II II
IIII
II1
IIII
II II
II II
Figure
II II
II II II II 1. A Cantorkm
triadic
II1
II II II II
bar.
Our objective is to evaluate the spatial variability of a relevant quantity affecting the dynamics in such a medium. Consider as an example porosity, 4, a typical property affecting transport. Assimilating, for instance, the pore space to the bars of the Cantor set (see Figure l), we want to estimate the relative importance of the standard deviation of 4 with respect to its mean. Because a pore of size r contributes to porosity (in a one-dimensional medium) a term of order r/L (L being the total length of the medium), the contribution of such pores to porosity is
where C is a proportionality factor and Df is the fractal dimension (Df = In 2/ In 3 N 0.63 for the Cantor bar). To estimate the variability of $, we need to know explicitly the pore size distribution, P(r). To this end, we notice that in the nth iteration length r and number N of the bars are given by
leading
to the Cantor
set the
1 r = - 1’0, 3” N = 2nNo, NO, rc being, we then have
respectively,
the initial
number
and length
N=Noen11121No
or, finally,
P(r) = Z-‘rMDf,
(5a) (5b) of the segments. Di
From equations
(5),
C. NICOLIS
30
wherein the normalization factor Z depends on the range of variation of r. Combining equation (4) and equation (7), we may now compute the mean porosity (4) and its variance (642) through
(4) =
Jrma’ dr r1-2DJ,
&
(84
TInin
(G2)= (412) - w2,
WI
with
Performing the integrals and setting r,in = 0, r,,,
(4) = g
and
= L = 1, we obtain:
(tM2)1’2= &.
In other words, the variability of porosity is comparable to the mean. This justifies the statistical view adopted here and suggests that a typical noise strength
(factor q in equation
(l)),
would be in the macroscopic range. One might also interpret this effect as some sort of effective “renormalized” diffusion whose rate is by orders of magnitude larger than molecular diffusion. Granted now that, fundamentally, the tracer dynamics is a random process described by equation (1) or (3), the following questions may naturally be raised [5,9]: - Starting at t = 0 on some horizontal surface r = Rs, what is the distribution of arrival times t of the tracers at a lower lying surface r = R ? Clearly, the dispersion of these first passage times will tell us whether or not the natural association between time and a given depth along a sediment, which is at the basis of the interpretation of the geological record, is a legitimate one. - Starting at t = 0 on the surface r = Ro, what is the distribution of the position r of the tracers at a later time t? Obviously, the properties of this front will give us a complementary view of whether the space and time courses follow a one-to-one correspondence or a more subtle relationship. - Starting at t = 0 with two tracers in nearby positions, what is the mean deviation of their trajectories for subsequent times ? Depending on its value, this deviation will pose limits on our ability to predict the instantaneous position of the tracer, and hence, to infer climate trends from the geological data. In the sequel, we provide an answer to these questions using a hierarchy of models of increasing complexity.
3. ONE-DIMENSIONAL
MODEL
In this section, we consider a one-dimensional version of equation (1) with w = volz, where z is directed along the vertical. Furthermore, for illustrative purposes, we limit ourselves to the simplest nontrivial case in which ~0 as well as the coupling term g are space-independent. Taking g = 1 without loss of generality and denoting q2/2 by D, the Fokker-Planck equation, equation (3), becomes d2P ---= t) D(10) at &2 In order to get a feeling about the behavior to be expected, we first consider an infinite system subjected to natural boundary conditions (P = 0, g = 0 at z = *co). Multiplying both sides of equation (10) successively by z and z2 and integrating over Z, one finds
wz,
-uo;P+
d (4 - dt = ‘UO+ (z) = d (6~~)
___
dt
where 6~ = t - (z).
.zo +
vat,
= 20 -+ (Sz2) = (6~“)~ + 2Dt,
(114 (lib)
31
Deciphering of Past Climates
In other words, as the tracer front advances, there is an increasingly large dispersion (proportional to the elapsed time) around the front. This suggests that there is an inherent mechanism of unpredictability in this system. Our next objective is to characterize it more sharply, using a more satisfactory model than the infinite medium model used above. We consider a column of vertical depth e. It is assumed that particles reaching the top z = 0 are trapped
within
the system,
whereas
particles
the lower inactive part of the sediment. following boundary conditions
reaching
This implies
the bottom
that equation
z = e are absorbed
by
(10) is to be solved with the
(P),d = 0. absorbing
reflecting To discuss
the first passage
time
problem
raised
(12)
in the preceding
section,
we need
instead
of P(z,t) the probability density g(z,t) that starting initially at z(z 2 0), the particle reaches the boundary z = e at time t. Notice that P(z, t) is associated to the stochastic variable z, whereas in g(z,t) the role of stochastic variable is played by t. One can show (see [5]) that g = -g, where G is the probability to be at time t still in the interval (O$) starting initially at z. This function obeys to the adjoint of equation (10): dG(z,
at together
with the boundary
conditions
t)
dG
d2G
=w()-++dZ
(13)
az2
(12) and the initial
condition
G(z,O) = 1.
Its solution can be found straightforwardly time probability is obtained in the form
dz,t) = g exp [-g
(z+
+0t)]
by Fourier
{exp (-g)
(14
series methods
+exp C nexp n odd
from which the first, passage
(g)} [-D
(z)2
t] sin “r(“,,
‘).
(15)
The details of the derivation of equation (15) starting from (13) are given in the Appendix. To estimate the parameters wg and D, we make use of the analogy, pointed out, in Section 2, between the Fokker-Planck equation (equation (10)) and the advection-diffusion equations describing transport of material in sediments. On this basis, we regard wo as the analogue of sedimentation velocity (considered as a constant in the approximation of the present section), and D as the analogue of an effective diffusion coefficient accounting for dispersion and mixing. A dimensionless quantity frequently used to compare the relative role of these two types of phenomena is the Peclet number P, = &o/D, where e is the system’s size. Ordinarily this parameter varies considerably depending on the type of the natural environment (see for instance, [lo, Table 21). For abyssal sediments, we choose for illustrative purposes ~0 = 1 cm Kyr-‘, yielding values of P, from 4 to 1. Figure 2a depicts the probability density g(.z, t) for reaching the boundary starting at z = 0, for three different values of D. We observe a dispersion of passage times measured by the variance (6t2), which is comparable to the mean and gets increasingly large as D (or q2) gets of D. small. However, the normalized dispersion (St2) U2 / (t) , is less sensitive to the variations An interesting manifestation of the role of stochastic effects is in the fact that the probability density attains its maximum at a value of time that is less than the deterministic time t = e/vo.
32
C. NICOLIS
P
0.2 0.15 0.1 0.05
0
10
5
15
t(Kyd
15
t(Kyr)
(a) zIo = 1 cmKyr_‘.
g 0.2 0.15 0.1 0.05
0
5
10 (b) o. - 0.
Figure 2. Probability density g(z, t) as computed from equation (15) for reaching the depth of the sediment e = 1Ocm starting at z = 0 for three different values of D in cm2 Kyr-’
I
g 0.12
0.08
0.04
0
5
10
15
20
twyr)
Figure 3. Probability density g(z, t) obtained by direct numerical simulation of the stochastic dynamics equation (1) using 10,000 particles. Parameters as in Figure 2a, with D = 2 cm* Kyr-l.
This result suggests that one should be extremely careful before adopting the usual (linear) correspondence between the elapsed time and the depth of a sediment. Notice that for fixed values of D and wo the shift of the maximum relative to t becomes less pronounced as the size of the system increases. Figure 2b depicts the probability density g(z, t) in the case where the drift velocity is negligible. We observe that for large values of D the difference with the preceding case is rather small. Figure 3 represents the form of g(z,t) obtained by direct numerical simulation of the stochastic dynamics of the tracer, equation (1) under the same conditions as in Figure 2a with
Deciphering of Past Climates
33
1 0.3 0.2 0.1
Figure 4. Probability density P(z, t) of arrival depths at t = 10 Kyr as deduced from numerical simulation of equation (1) involving 10,000 particles starting at z = 0. Parameter values as in Figure 2a with D = 5 cm2 Kyr-‘.
D = 2.5cm’Kyr-‘. To obtain a good statistic a number of N = 10,000 realizations of the process have been carried out. The overall agreement with the analytical results is quite good. Let us now turn to the distribution of arrival depths at a given time. Figure 4 depicts the results of a numerical simulation of the stochastic dynamics involving N = 10,000 realizations, each starting on the surface z = 0 and running over 10 Kyrs. We again obtain an appreciable dispersion since, for one thing, in the absence of fluctuations and for ZIO= 1 cm Kyr-’ all realizations should be at depth z = -10 cm.
4. SENSITIVITY
TO INITIAL
CONDITIONS
We now turn to a complementary view in which attention is focused on the realizations of the stochastic process (1) rather than on its probability distribution. The specific question we raise is, how two stochastic trajectories starting with initial conditions differing by a small amount E will deviate-if at all-at subsequent times. In other words, to what extent can a reliable prediction of the tracer position be made on the basis of the knowledge of its initial position up to an error margin E? Clearly, the answer to this question will condition the ability to assess reliably climatic trends from geological data. Under the simplifying assumptions spelled out in Section 3, equation (1) becomes dz - = -uJ + qF(t), dt where F(t) is a Gaussian white noise satisfying equations (2). Integrating equation (16) using two different initial conditions, we obtain
zl(t)
= zo -uot
+ q
t J t J
dt’Fl(t’),
0
(17)
dt’Fz(t’),
zz(t) = ze + E - vat + Q
0
where it is understood that Fl, Fz are two different realizations of the white noise process. Subtracting the two equations (17), we obtain the mean instantaneous error uE(t) as z&(t) = Izg(t) - z1(t)l = IE+ q(W2(t) where
Wz(t) =Jt
dt’Fi(t’)
0
wl(t))l,
(18b)
34
C. NICOLIS
is a Wiener process. As expected, equation (Isa) fluctuates considerably over time. To obtain a systematic behavior, we average over the probability distribution of the process
P(W)=le
-w2/4q=t
(19)
&q
The result is (see [11,12])
(&(t))
= 2
J-
$
e-E2/4q2t + Eerf
& (
We notice that beyond a characteristic
(20) )
time t* of the order of t*
M
”
(214
492 ’
the instantaneous
error increases monotonically (l&(t))
according to
= g
@lb)
P/2.
The initial error E is thus amplified in time by a power law. Although this is a milder form of amplification than the exponential sensitivity to initial conditions characterizing deterministic chaos [13], it does imply that, time going on, the uncertainty on the position of the tracer is steadily increasing. For times T of the order of (or beyond) r M .rr2/4q2 the uncertainty is of the order of unity, which is about the distance traveled by the tracer in a geological time unit. Notice that r may be quite small in this scale-for instance, using the same numerical values as in Figure 4 7 M 0.25 Kyr. Figure 5 depicts the time evolution of the mean error as obtained from equation (20) (solid line) and from direct numerical simulation of equation (16) with q” = 5/4cm2 Kyr-l and ]E] = 1 cm averaged over 2,000 samples (open circles). We observe that the initial regime merges with the one given by (21b) at a time t* of about 0.2Kyr in agreement with the theoretical estimate.
1.2
0.8 i
0
0.5
1
Q/2
Figure 5. Time evolution of the mean error (uE) as obtained from equation (20) (solid line) and from direct numerical simulation of equation (16) averaged over 2,000 samples (open circles), with q2 = 1.25cm2 Kyr-’ and 1~1= Icm.
In reality, it is unlikely that errors will grow indefinitely. As deeper sediments become more compact, the tracer will be less mobile. Such a “saturation” cannot be taken into account in the setting of the present model. It would certainly be interesting to elaborate to that effect an enlarged description accounting in a more realistic manner for the depth dependence of the structural properties of the medium.
Deciphering of Past Climates
5. TWO-DIMENSIONAL Despite
the variability
the preceding of random
sections,
in the distribution
walk remain
smooth.
In the present
from equation
qx,
z, t)
at
and arrival times of the tracer
system considered section,
quite different in a two-dimensional medium. Let x denote the horizontal direction. Taking, gZ = gZ = 1, we obtain
MODEL
of positions
in the one-dimensional
35
therein,
we show that
as in Section
the statistical the situation
reported
in
properties may become
3, v = uolz,D = q2/2 and
(3) dP = -Qz+T
q2 (
@P p+=
d2P
>
(22)
.
We want to explore the new features arising from the presence of the second (horizontal) direction. To get as transparent a view as possible, we begin by a simulation of random walk with ve N 0 on a discrete The simulation
two-dimensional
is carried
moves by one lattice be visited is decided
lattice
out in discrete
whose continuous
steps, t, = nAt,
limit corresponds
to equation
(10).
as follows. At each time unit the tracer
distance toward one of its first neighbors. Which particular neighbor will by using a random number generator with uniform distribution. Figure 6
depicts the results. We observe that, as time grows, the tracer visits the lattice in an intermittent fashion [8,9]. Specifically, it first remains in a certain neighborhood for an appreciable amount of time covering it in a more or less dense manner and subsequently transits rather quickly toward another neighborhood in which it exhibits a similar behavior. Given an arbitrarily large amount of time, all regions of the lattice will be visited with probability one; however in a finite amount of time as in a typical real-world problem the intermittent character of the two-dimensional dynamics will entail a markedly nonuniform distribution of the tracer. We suggest that this feature may be related to the hiatuses characterizing certain geological deposits [14].
-20 -30 -40
Figure 6. Stochastic trajectory of a particle performing random walk on a 2-d lattice, starting at z = y = 0 and moving during 2,000 time steps.
Figures 7a to 7c report an alternative manifestation of the new features arising from the existence of a second dimension. The figures register the position, at times t = 500 (Figure 7a), t = 1,500 (Figure 7b) and t = 6,000 steps (Figure 7c) of N = 20,000 particles performing random walk and starting initially at z = 0. The random walk is modeled as in Figure 6 taking into account the boundary conditions, equations (12). We observe the appearance of a fractal diffusion front [15,16], separating a region of densely occupied sites from a region of sparse occupation. This separation is not clear-cut: the front has a complex, irregular form in space and, furthermore, broadens as time advances. Two types of inhomogeneities are formed during this process: First,, islands of unoccupied sites surrounded by regions of occupied ones and vice versa, and second, connected clusters of unoccupied sites invading occupied regions, as in classical percolation problems.
36
C. NICOLIS
-80
I -20
20
I
x
(b) t = I:000
-20
20 (c)
t =
x
GPO00
Figure 7. Spatial distribution of 20,000 particles performing random walk in a 2-d lattice initially at .z = 0 and -40 5 2 5 40 after the appropriate time steps. Notice the high degree of fragmentation of the front and the increase of its width with time.
In summary, it becomes increasingly difficult to predict, in an unambiguous manner, which characteristic depth will be reached in a given lapse of time: tracer deposits in certain parts of two horizontal layers at different depths could correspond to a chronological order which is different from, or even opposite to the one suggested by sheer inspection of the depths.
6. DISCUSSION In this paper, we developed a stochastic model of tracer dynamics in sediments accounting for the complex, disordered character of geological materials. We have worked out explicitly, both analytically and by direct numerical simulation, the simplest version of the model. We found that the complexity of the medium entails a large dispersion of arrival times of the tracer at a given depth as well as of the depth along the sediment at a given time. This phenomenon becomes especially pronounced in two dimensions where the precipitation front displays a considerable degree of fragmentation. As a corollary of the above, in regions of the sediment in which no absolute dating is possible, the usually admitted linear relationship between the elapsed time and the depth of a sediment appears to be seriously compromised. In view of the potential importance of this conclusion more elaborate models accounting realistically for the fractal properties of the geological medium would be worth developing in the future. One of the most challenging conclusions of our analysis pertains to the interpretation of deepsea core data. It is well known that the oxygen isotope composition inferred from these data serves as a most valuable indicator of paleoclimates. Furthermore, it is believed that a number of
Deciphering
climatic
episodes
time scale.
have not been gradual
Clearly,
if, as suggested
of Past Climates
but have, rather,
in the present
the position of the tracer along the sediment series data becomes necessary.
37
occurred
abruptly
on a relatively
work, there exists an inherent
on a similar
short
uncertainty
time scale, a reassessment
on
of the time
APPENDIX Derivation
of Equation
(15)
In view of the relation
g = -2
probability (equation
density
pointed
g it is sufficient
out in Section
to find the solution,
G(z,t)
the first passage
of the adjoint
time
Fokker-Planck
(13)) dG(z,t)
to the following
initial
aG
d2G
az
az2 1
=vo---+D-
at subject
3, to compute
and boundary
conditions
(A.1)
(eqs. (14) and (12)):
G(z,O) = 1,
= az > (Wz,t)
(A4
0
(A.3a)
7
z=o
(G(zz, t)),+ = 0. It will be convenient to extend the spatial domain satisfying to the symmetric boundary conditions G(4,
Furthermore, transformation
as an intermediate
Inserting
(A.5) into (A.l),
one easily finds that
and boundary
f(z,t).
f(z, t) satisfies
of the
(A.5) the diffusion
equation
az2
(A.6)
,
(A.7)
= 0.
(A.8)
= exp (2)
f(-l,O)
= f(&O)
f in Fourier
series,
f(z, t) = c
A,(t)
sin Ku
n Inserting into (A.6) and using the orthogonality (-e, e) one easily finds . A,(t)
by means
conditions f(z,O)
Owing to (A.8), we expand
(A.4)
the drift term in (A.l)
(z+tuot)]
[-5
at
with the initial
seek for solutions
t) = G(C, t) = 0.
df&Y
together
(0, e) to (-e, e) and
step, we will eliminate
G=exp
(A.3b)
= A,(O)exp
properties
[-D
($)’
“‘I. of sin[(rn(L! - ~))/(2[)]
t] .
(A.9
in the interval
(A.lO)
C. NIC~LIS
38
To determine in (A.9)
A,(O),
we make use of the initial condition
equation
(A.7).
Inversion
of the series
yields
=-
Combining
(A.9)
G(z,t)
212 (s)‘: (5)’
to (A.ll)
= &exp
and (A.5),
[ -20‘O
{exp (-g) {exp
(z+iu0t)]
Differentiating
with respect
cgj”:
(-2)
+exp
($)‘t]
(A.ll)
(g)}
we arrive at the following explicit
exp [-D x 2
+exp
form for G(z,t)
(s)}
sin nT’~~
‘).
(A’12)
($j2
to t, one finally obtains
equation
(15).
REFERENCES in the modeling of early diagnosis, Earth Planet. Sci. Lett. 29, 33-340, (1976). N.L. Guinasso and D.R. Schink, Quantitative estimates of biological mixing rates in abyssal sediments, J. Geophys. Rev. 80, 3032-3043, (1975). J.B. Fisher, W.J. Lick, P.L. McCall and J.A. Robbins, Vertical mixing of lake sediments by turbificid oligochsetes, J. Geophys. Res. 85, 3997-4006, (1980). B.P. Boudreau and D.M. Imboden, Mathematics of tracer mixing in sediments: III. The theory of nonlocal mixing within sediments, Am. Jour. Sci. 287, 693-719, (1987). C.W. Gardiner, Handbook of Stochastic Methods, Springer Verlag, Berlin, (1983). P.-Z. Wong, The statistical physics of sedimentary rock, Physics Today, 24-32, (December 1988). P. Friedlingstein, Sedimentation dans un milieu complexe: Modelisation mathematique et application iL l’interpretation des don&es geologiques et climatiques. Memoire de Licence, Faculte des Sciences Appliqu&s, Universite Libre de Bruxelles. (1989). B. Mandelbrot, fiactals, Form, Chance and Dimension, Freeman, San Francisco, (1977). W. Feller, An Introductzon to Probability Theory, Vol. I, Wiley, New York, (1968). B.P. Boudreau, Mathematics of tracer mixing in sediments I. Spatially-dependent, diffusive mixing, Am. Jour. Sci. 286, 161-198, (1986). K. Seki, V. Balakrishnan and G. Nicolis, Sensitivity to initial conditions in stochastic systems, Phys. Rev. E47, 155-163, (1993). G. Nicolis, C. Nicolis and K. Seki, Sensitivity to initial conditions in complex systems, In Statistical Physics and Thermodynamics of Nonlinear Nonequilibrium Systems, (Edited by W. Ebeling and W. Muschik), pp. 23-46, (1993). H. Schuster, Determznzstic Chaos, V.C.H. Verlag, Weinheim, (1988). R.E. Plotnik, A fractal model for the distribution of stratigraphic hiatuses, Jour. Geology 94, 885-890, (1986). B. Sapoval, M. Rosso and J.F. Gouyet, The fractal nature of a diffusing front and the relation to percolation, J. Phys. Lett. 48, L149-L156, (1985). J. Feder, Fractals, Plenum, New York. (1988). E.R. Christensen and P.K. Bhunia, Modeling radiotracers in sediments: Comparison with observations in lakes Huron and Michigan, J. Geophys. Res. 91, 855998571, (1986).
1. R.A. Berner, Inclusion of adsorption 2.
3. 4. 5. 6. 7.
8. 9. 10. 11. 12.
13. 14. 15. 16. 17.
Compvt.
Modelling Vol. 21, No. 6, pp. 27-38, 1995 Copyright@1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0895-7177/95 $9.50 + 0.00
Pergamon
08957177(95)00021-6
Tracer Dynamics in Ocean Sediments the Deciphering of Past Climates
and
C. NICOLIS Institut
Royal
3, Avenue
Mktkorologique
Circulaire,
de Belgique
1180 Bruxelies,
Belgium
(Received and accepted April 1 L&X$) Abstract-An analysis of tracer dynamics in ocean sediments taking into account the disordered, fractal, character of geological materials is performed. The complexity of the medium is modeled by augmenting the equation of motion of the tracer through the addition of a stochastic forcing. In the one-dimensional version of the model, it is found both analytically and numerically that the dynamics is characterized by a large dispersion of arrival times at a given depth of the sediment as deduced from the first passage time probability distribution of the tracer. Moreover, the trajectories of two initially close tracer particles deviate at later times according to a power law, thereby raising the question of reliability of climatic trends deduced from geological data. Numerical simulations in a more realistic, two-dimensional setting show that this uncertainty becomes even more pronounced entailing the formation of a fractal precipitation front whose width increases with time. The implications of the results on the interpretation of geological records are discussed.
Keywords-Ractals,
Glaciations,
Predictability,
Random
walk
1. INTRODUCTION Much of our understanding of the earth’s past environmental and climatic conditions rests on our ability to decipher geological data, particularly those pertaining to ocean sediment formation. Now, reconstructing a time sequence of events from the spatial distribution of a representative tracer within the sediment involves, perforce, a number of highly nontrivial intermediate steps. Specifically, since as a rule only a limited number of points can be dated with a good accuracy, a model of sediment formation has to be adopted before the “space axis” along the sediment can be converted into a “time axis” along which past history begins to unfold. A classical approach to tracer distribution within sediments makes use of advection-diffusion models, sometimes augmented by the inclusion of adsorption or production-destruction terms. The important issue of the mixing of biological origin is usually addressed through a variety of phenomenological assumptions relative to the depth dependence of the diffusion coefficients involved in the theory [l-4]. Typically, all these models adopt a deterministic description limited to a one-dimensional view, in which attention is focused on the tracer distribution along the direction of advection. Our principal objective in the present work is to develop an enlarged description of tracer dynamics within a sediment, which incorporates the influence of stochastic effects. The principal motivation behind this idea is the fact that geological materials are complex media characterized This work is supported, in part, by the Commission Scientific, Technical and Cultural Affairs (O.S.T.C.).
of the European
Communities
and by the Federal Office for
Typeset Ia4 21:6-C
27
by dn/ls-QX
C. NICOLIS
28
by a pronounced disorder. Within such a medium, the movement of a tracer is expected to become, effectively, a random process since the instantaneous direction and speed are likely to vary continuously depending on the locally prevailing conditions. This should entail, in turn, that the position of the tracer after a sufficiently long lapse of time will be subjected to a considerable uncertainty. We want to characterize this uncertainty in a quantitative manner and to explore its consequences in the interpretation
of the geological record, especially the one pertaining
to
climatic evolution. In Section 2, we discuss in some detail the origin of the randomness of geological materials and formulate the basic equations governing the dynamics in such a material.
A one-dimensional
version of these equations is analyzed in Section 3. In particular, we evaluate the probability density of the position of the tracer as well as the first passage time distribution determining the times of arrival at a given depth. Both distributions a large uncertainty
display a considerable dispersion, entailing
on the position of the tracer, in agreement with the comments made above.
A complementary view of this uncertainty is explored in Section 4, where we study the mean deviation of the trajectories of two initially nearby tracers. As it turns out, this deviation increases steadily in time according to a power law. This illustrates further the problems inherent in the deciphering of climatic trends from the geological record. In Section 5, the analysis is extended to two dimensions. It is found that the presence of a second dimension gives rise to a new mechanism of variability in the form of a fractal diffusion front whose average width increases with time. The implications of the results are discussed in Section 6, in connection with the glacial-interglacial climatic transitions.
2. GENERAL
FORMULATION
Consider a tracer moving in a geological medium. Depending on the circumstances, it may represent a chemical substance dissolved in water, a macroscopic particle (like, for instance, debris of a certain biological species) containing a radioactive isotope, etc. In either case, owing to the disorder of the medium-including the one arising from the action of living organisms-we stipulate that, fundamentally, this movement is to be regarded as a random walk. Denoting by r the instantaneous position of the tracer and by t the time, we therefore write
dr - = v(r) + qg(r).F(t). dt Here v denotes the advection (drift) velocity; F a random force accounting for the action of the medium on the particle; g a coupling term; and q the strength of the coupling. The complexity of the medium can be further manifested through the fact that v and g may have a complicated r-dependence, which may even be viewed most adequately as a random function of r. We do not allow for an explicit time dependence of v and g as this is taken into account by the random force F(t). We shall make on F(t)
the simplest possible assumption of a (normalized)
Gaussian white
noise:
Pi(t)) = 0 (Fi(t)Fj(t'))
=
6;; b(t -
(2)
t’)
We emphasize that these assumptions can be relaxed, but in the sequel, we shall not be concerned with this generalization. It is well known that equations (1) and (2) are equivalent to a Fokker-Planck equation for the probability density P(r, t) of finding the tracer around point r at time t (see [5]):
at=-c& w(r)P(r, t) + $ C dP
i
2
ijk
d gik & dr,,
gjk
P(r, t)
(3)
Deciphering
of Past
Climates
29
Notice the formal similarity between equation (3) and the advection-diffusion equations describing the transport of material in sediments: Here the role of advection is played by the drift velocity V, whereas
the term g coupling
diffusion
whose magnitude
the tracer depends
movement
on the strength
to the random parameter
force gives rise to an effective
q.
Let us comment in somewhat more detail, on the connection between the complexity of the medium and the random character of the dynamics. It is well known that geological media present a very pronounced irregularity and fragmentation, which can hardly be viewed as a small, quantitative modification of some regular “reference” crystal lattice. To capture the essence of this irregularity, we assimilate the material to a fractal [6,7]. The specific choice of fractal is a major question whose answer depends on the nature of the material intervening in a given problem. Here we consider an abstract-but typical-model in which the medium is regarded as a triadic Cantor bar, as depicted in Figure 1 (see [S]). Specifically, it is assumed that all bar sizes between two extreme values r,in and T,~,, are present with a probability P(r) induced by the rules leading to the construction of the Cantor set.
-
m
m
-
I
II
II
II
II
II
II
II
II
II
II
II
II
II II II II
II II II II
II II II II
II II II II
II II II II
II1
IIII
II II
IIII
II1
IIII
II II
II II
Figure
II II
II II II II 1. A Cantorkm
triadic
II1
II II II II
bar.
Our objective is to evaluate the spatial variability of a relevant quantity affecting the dynamics in such a medium. Consider as an example porosity, 4, a typical property affecting transport. Assimilating, for instance, the pore space to the bars of the Cantor set (see Figure l), we want to estimate the relative importance of the standard deviation of 4 with respect to its mean. Because a pore of size r contributes to porosity (in a one-dimensional medium) a term of order r/L (L being the total length of the medium), the contribution of such pores to porosity is
where C is a proportionality factor and Df is the fractal dimension (Df = In 2/ In 3 N 0.63 for the Cantor bar). To estimate the variability of $, we need to know explicitly the pore size distribution, P(r). To this end, we notice that in the nth iteration length r and number N of the bars are given by
leading
to the Cantor
set the
1 r = - 1’0, 3” N = 2nNo, NO, rc being, we then have
respectively,
the initial
number
and length
N=Noen11121No
or, finally,
P(r) = Z-‘rMDf,
(5a) (5b) of the segments. Di
From equations
(5),
C. NICOLIS
30
wherein the normalization factor Z depends on the range of variation of r. Combining equation (4) and equation (7), we may now compute the mean porosity (4) and its variance (642) through
(4) =
Jrma’ dr r1-2DJ,
&
(84
TInin
(G2)= (412) - w2,
WI
with
Performing the integrals and setting r,in = 0, r,,,
(4) = g
and
= L = 1, we obtain:
(tM2)1’2= &.
In other words, the variability of porosity is comparable to the mean. This justifies the statistical view adopted here and suggests that a typical noise strength
(factor q in equation
(l)),
would be in the macroscopic range. One might also interpret this effect as some sort of effective “renormalized” diffusion whose rate is by orders of magnitude larger than molecular diffusion. Granted now that, fundamentally, the tracer dynamics is a random process described by equation (1) or (3), the following questions may naturally be raised [5,9]: - Starting at t = 0 on some horizontal surface r = Rs, what is the distribution of arrival times t of the tracers at a lower lying surface r = R ? Clearly, the dispersion of these first passage times will tell us whether or not the natural association between time and a given depth along a sediment, which is at the basis of the interpretation of the geological record, is a legitimate one. - Starting at t = 0 on the surface r = Ro, what is the distribution of the position r of the tracers at a later time t? Obviously, the properties of this front will give us a complementary view of whether the space and time courses follow a one-to-one correspondence or a more subtle relationship. - Starting at t = 0 with two tracers in nearby positions, what is the mean deviation of their trajectories for subsequent times ? Depending on its value, this deviation will pose limits on our ability to predict the instantaneous position of the tracer, and hence, to infer climate trends from the geological data. In the sequel, we provide an answer to these questions using a hierarchy of models of increasing complexity.
3. ONE-DIMENSIONAL
MODEL
In this section, we consider a one-dimensional version of equation (1) with w = volz, where z is directed along the vertical. Furthermore, for illustrative purposes, we limit ourselves to the simplest nontrivial case in which ~0 as well as the coupling term g are space-independent. Taking g = 1 without loss of generality and denoting q2/2 by D, the Fokker-Planck equation, equation (3), becomes d2P ---= t) D(10) at &2 In order to get a feeling about the behavior to be expected, we first consider an infinite system subjected to natural boundary conditions (P = 0, g = 0 at z = *co). Multiplying both sides of equation (10) successively by z and z2 and integrating over Z, one finds
wz,
-uo;P+
d (4 - dt = ‘UO+ (z) = d (6~~)
___
dt
where 6~ = t - (z).
.zo +
vat,
= 20 -+ (Sz2) = (6~“)~ + 2Dt,
(114 (lib)
31
Deciphering of Past Climates
In other words, as the tracer front advances, there is an increasingly large dispersion (proportional to the elapsed time) around the front. This suggests that there is an inherent mechanism of unpredictability in this system. Our next objective is to characterize it more sharply, using a more satisfactory model than the infinite medium model used above. We consider a column of vertical depth e. It is assumed that particles reaching the top z = 0 are trapped
within
the system,
whereas
particles
the lower inactive part of the sediment. following boundary conditions
reaching
This implies
the bottom
that equation
z = e are absorbed
by
(10) is to be solved with the
(P),d = 0. absorbing
reflecting To discuss
the first passage
time
problem
raised
(12)
in the preceding
section,
we need
instead
of P(z,t) the probability density g(z,t) that starting initially at z(z 2 0), the particle reaches the boundary z = e at time t. Notice that P(z, t) is associated to the stochastic variable z, whereas in g(z,t) the role of stochastic variable is played by t. One can show (see [5]) that g = -g, where G is the probability to be at time t still in the interval (O$) starting initially at z. This function obeys to the adjoint of equation (10): dG(z,
at together
with the boundary
conditions
t)
dG
d2G
=w()-++dZ
(13)
az2
(12) and the initial
condition
G(z,O) = 1.
Its solution can be found straightforwardly time probability is obtained in the form
dz,t) = g exp [-g
(z+
+0t)]
by Fourier
{exp (-g)
(14
series methods
+exp C nexp n odd
from which the first, passage
(g)} [-D
(z)2
t] sin “r(“,,
‘).
(15)
The details of the derivation of equation (15) starting from (13) are given in the Appendix. To estimate the parameters wg and D, we make use of the analogy, pointed out, in Section 2, between the Fokker-Planck equation (equation (10)) and the advection-diffusion equations describing transport of material in sediments. On this basis, we regard wo as the analogue of sedimentation velocity (considered as a constant in the approximation of the present section), and D as the analogue of an effective diffusion coefficient accounting for dispersion and mixing. A dimensionless quantity frequently used to compare the relative role of these two types of phenomena is the Peclet number P, = &o/D, where e is the system’s size. Ordinarily this parameter varies considerably depending on the type of the natural environment (see for instance, [lo, Table 21). For abyssal sediments, we choose for illustrative purposes ~0 = 1 cm Kyr-‘, yielding values of P, from 4 to 1. Figure 2a depicts the probability density g(.z, t) for reaching the boundary starting at z = 0, for three different values of D. We observe a dispersion of passage times measured by the variance (6t2), which is comparable to the mean and gets increasingly large as D (or q2) gets of D. small. However, the normalized dispersion (St2) U2 / (t) , is less sensitive to the variations An interesting manifestation of the role of stochastic effects is in the fact that the probability density attains its maximum at a value of time that is less than the deterministic time t = e/vo.
32
C. NICOLIS
P
0.2 0.15 0.1 0.05
0
10
5
15
t(Kyd
15
t(Kyr)
(a) zIo = 1 cmKyr_‘.
g 0.2 0.15 0.1 0.05
0
5
10 (b) o. - 0.
Figure 2. Probability density g(z, t) as computed from equation (15) for reaching the depth of the sediment e = 1Ocm starting at z = 0 for three different values of D in cm2 Kyr-’
I
g 0.12
0.08
0.04
0
5
10
15
20
twyr)
Figure 3. Probability density g(z, t) obtained by direct numerical simulation of the stochastic dynamics equation (1) using 10,000 particles. Parameters as in Figure 2a, with D = 2 cm* Kyr-l.
This result suggests that one should be extremely careful before adopting the usual (linear) correspondence between the elapsed time and the depth of a sediment. Notice that for fixed values of D and wo the shift of the maximum relative to t becomes less pronounced as the size of the system increases. Figure 2b depicts the probability density g(z, t) in the case where the drift velocity is negligible. We observe that for large values of D the difference with the preceding case is rather small. Figure 3 represents the form of g(z,t) obtained by direct numerical simulation of the stochastic dynamics of the tracer, equation (1) under the same conditions as in Figure 2a with
Deciphering of Past Climates
33
1 0.3 0.2 0.1
Figure 4. Probability density P(z, t) of arrival depths at t = 10 Kyr as deduced from numerical simulation of equation (1) involving 10,000 particles starting at z = 0. Parameter values as in Figure 2a with D = 5 cm2 Kyr-‘.
D = 2.5cm’Kyr-‘. To obtain a good statistic a number of N = 10,000 realizations of the process have been carried out. The overall agreement with the analytical results is quite good. Let us now turn to the distribution of arrival depths at a given time. Figure 4 depicts the results of a numerical simulation of the stochastic dynamics involving N = 10,000 realizations, each starting on the surface z = 0 and running over 10 Kyrs. We again obtain an appreciable dispersion since, for one thing, in the absence of fluctuations and for ZIO= 1 cm Kyr-’ all realizations should be at depth z = -10 cm.
4. SENSITIVITY
TO INITIAL
CONDITIONS
We now turn to a complementary view in which attention is focused on the realizations of the stochastic process (1) rather than on its probability distribution. The specific question we raise is, how two stochastic trajectories starting with initial conditions differing by a small amount E will deviate-if at all-at subsequent times. In other words, to what extent can a reliable prediction of the tracer position be made on the basis of the knowledge of its initial position up to an error margin E? Clearly, the answer to this question will condition the ability to assess reliably climatic trends from geological data. Under the simplifying assumptions spelled out in Section 3, equation (1) becomes dz - = -uJ + qF(t), dt where F(t) is a Gaussian white noise satisfying equations (2). Integrating equation (16) using two different initial conditions, we obtain
zl(t)
= zo -uot
+ q
t J t J
dt’Fl(t’),
0
(17)
dt’Fz(t’),
zz(t) = ze + E - vat + Q
0
where it is understood that Fl, Fz are two different realizations of the white noise process. Subtracting the two equations (17), we obtain the mean instantaneous error uE(t) as z&(t) = Izg(t) - z1(t)l = IE+ q(W2(t) where
Wz(t) =Jt
dt’Fi(t’)
0
wl(t))l,
(18b)
34
C. NICOLIS
is a Wiener process. As expected, equation (Isa) fluctuates considerably over time. To obtain a systematic behavior, we average over the probability distribution of the process
P(W)=le
-w2/4q=t
(19)
&q
The result is (see [11,12])
(&(t))
= 2
J-
$
e-E2/4q2t + Eerf
& (
We notice that beyond a characteristic
(20) )
time t* of the order of t*
M
”
(214
492 ’
the instantaneous
error increases monotonically (l&(t))
according to
= g
@lb)
P/2.
The initial error E is thus amplified in time by a power law. Although this is a milder form of amplification than the exponential sensitivity to initial conditions characterizing deterministic chaos [13], it does imply that, time going on, the uncertainty on the position of the tracer is steadily increasing. For times T of the order of (or beyond) r M .rr2/4q2 the uncertainty is of the order of unity, which is about the distance traveled by the tracer in a geological time unit. Notice that r may be quite small in this scale-for instance, using the same numerical values as in Figure 4 7 M 0.25 Kyr. Figure 5 depicts the time evolution of the mean error as obtained from equation (20) (solid line) and from direct numerical simulation of equation (16) with q” = 5/4cm2 Kyr-l and ]E] = 1 cm averaged over 2,000 samples (open circles). We observe that the initial regime merges with the one given by (21b) at a time t* of about 0.2Kyr in agreement with the theoretical estimate.
1.2
0.8 i
0
0.5
1
Q/2
Figure 5. Time evolution of the mean error (uE) as obtained from equation (20) (solid line) and from direct numerical simulation of equation (16) averaged over 2,000 samples (open circles), with q2 = 1.25cm2 Kyr-’ and 1~1= Icm.
In reality, it is unlikely that errors will grow indefinitely. As deeper sediments become more compact, the tracer will be less mobile. Such a “saturation” cannot be taken into account in the setting of the present model. It would certainly be interesting to elaborate to that effect an enlarged description accounting in a more realistic manner for the depth dependence of the structural properties of the medium.
Deciphering of Past Climates
5. TWO-DIMENSIONAL Despite
the variability
the preceding of random
sections,
in the distribution
walk remain
smooth.
In the present
from equation
qx,
z, t)
at
and arrival times of the tracer
system considered section,
quite different in a two-dimensional medium. Let x denote the horizontal direction. Taking, gZ = gZ = 1, we obtain
MODEL
of positions
in the one-dimensional
35
therein,
we show that
as in Section
the statistical the situation
reported
in
properties may become
3, v = uolz,D = q2/2 and
(3) dP = -Qz+T
q2 (
@P p+=
d2P
>
(22)
.
We want to explore the new features arising from the presence of the second (horizontal) direction. To get as transparent a view as possible, we begin by a simulation of random walk with ve N 0 on a discrete The simulation
two-dimensional
is carried
moves by one lattice be visited is decided
lattice
out in discrete
whose continuous
steps, t, = nAt,
limit corresponds
to equation
(10).
as follows. At each time unit the tracer
distance toward one of its first neighbors. Which particular neighbor will by using a random number generator with uniform distribution. Figure 6
depicts the results. We observe that, as time grows, the tracer visits the lattice in an intermittent fashion [8,9]. Specifically, it first remains in a certain neighborhood for an appreciable amount of time covering it in a more or less dense manner and subsequently transits rather quickly toward another neighborhood in which it exhibits a similar behavior. Given an arbitrarily large amount of time, all regions of the lattice will be visited with probability one; however in a finite amount of time as in a typical real-world problem the intermittent character of the two-dimensional dynamics will entail a markedly nonuniform distribution of the tracer. We suggest that this feature may be related to the hiatuses characterizing certain geological deposits [14].
-20 -30 -40
Figure 6. Stochastic trajectory of a particle performing random walk on a 2-d lattice, starting at z = y = 0 and moving during 2,000 time steps.
Figures 7a to 7c report an alternative manifestation of the new features arising from the existence of a second dimension. The figures register the position, at times t = 500 (Figure 7a), t = 1,500 (Figure 7b) and t = 6,000 steps (Figure 7c) of N = 20,000 particles performing random walk and starting initially at z = 0. The random walk is modeled as in Figure 6 taking into account the boundary conditions, equations (12). We observe the appearance of a fractal diffusion front [15,16], separating a region of densely occupied sites from a region of sparse occupation. This separation is not clear-cut: the front has a complex, irregular form in space and, furthermore, broadens as time advances. Two types of inhomogeneities are formed during this process: First,, islands of unoccupied sites surrounded by regions of occupied ones and vice versa, and second, connected clusters of unoccupied sites invading occupied regions, as in classical percolation problems.
36
C. NICOLIS
-80
I -20
20
I
x
(b) t = I:000
-20
20 (c)
t =
x
GPO00
Figure 7. Spatial distribution of 20,000 particles performing random walk in a 2-d lattice initially at .z = 0 and -40 5 2 5 40 after the appropriate time steps. Notice the high degree of fragmentation of the front and the increase of its width with time.
In summary, it becomes increasingly difficult to predict, in an unambiguous manner, which characteristic depth will be reached in a given lapse of time: tracer deposits in certain parts of two horizontal layers at different depths could correspond to a chronological order which is different from, or even opposite to the one suggested by sheer inspection of the depths.
6. DISCUSSION In this paper, we developed a stochastic model of tracer dynamics in sediments accounting for the complex, disordered character of geological materials. We have worked out explicitly, both analytically and by direct numerical simulation, the simplest version of the model. We found that the complexity of the medium entails a large dispersion of arrival times of the tracer at a given depth as well as of the depth along the sediment at a given time. This phenomenon becomes especially pronounced in two dimensions where the precipitation front displays a considerable degree of fragmentation. As a corollary of the above, in regions of the sediment in which no absolute dating is possible, the usually admitted linear relationship between the elapsed time and the depth of a sediment appears to be seriously compromised. In view of the potential importance of this conclusion more elaborate models accounting realistically for the fractal properties of the geological medium would be worth developing in the future. One of the most challenging conclusions of our analysis pertains to the interpretation of deepsea core data. It is well known that the oxygen isotope composition inferred from these data serves as a most valuable indicator of paleoclimates. Furthermore, it is believed that a number of
Deciphering
climatic
episodes
time scale.
have not been gradual
Clearly,
if, as suggested
of Past Climates
but have, rather,
in the present
the position of the tracer along the sediment series data becomes necessary.
37
occurred
abruptly
on a relatively
work, there exists an inherent
on a similar
short
uncertainty
time scale, a reassessment
on
of the time
APPENDIX Derivation
of Equation
(15)
In view of the relation
g = -2
probability (equation
density
pointed
g it is sufficient
out in Section
to find the solution,
G(z,t)
the first passage
of the adjoint
time
Fokker-Planck
(13)) dG(z,t)
to the following
initial
aG
d2G
az
az2 1
=vo---+D-
at subject
3, to compute
and boundary
conditions
(A.1)
(eqs. (14) and (12)):
G(z,O) = 1,
= az > (Wz,t)
(A4
0
(A.3a)
7
z=o
(G(zz, t)),+ = 0. It will be convenient to extend the spatial domain satisfying to the symmetric boundary conditions G(4,
Furthermore, transformation
as an intermediate
Inserting
(A.5) into (A.l),
one easily finds that
and boundary
f(z,t).
f(z, t) satisfies
of the
(A.5) the diffusion
equation
az2
(A.6)
,
(A.7)
= 0.
(A.8)
= exp (2)
f(-l,O)
= f(&O)
f in Fourier
series,
f(z, t) = c
A,(t)
sin Ku
n Inserting into (A.6) and using the orthogonality (-e, e) one easily finds . A,(t)
by means
conditions f(z,O)
Owing to (A.8), we expand
(A.4)
the drift term in (A.l)
(z+tuot)]
[-5
at
with the initial
seek for solutions
t) = G(C, t) = 0.
df&Y
together
(0, e) to (-e, e) and
step, we will eliminate
G=exp
(A.3b)
= A,(O)exp
properties
[-D
($)’
“‘I. of sin[(rn(L! - ~))/(2[)]
t] .
(A.9
in the interval
(A.lO)
C. NIC~LIS
38
To determine in (A.9)
A,(O),
we make use of the initial condition
equation
(A.7).
Inversion
of the series
yields
=-
Combining
(A.9)
G(z,t)
212 (s)‘: (5)’
to (A.ll)
= &exp
and (A.5),
[ -20‘O
{exp (-g) {exp
(z+iu0t)]
Differentiating
with respect
cgj”:
(-2)
+exp
($)‘t]
(A.ll)
(g)}
we arrive at the following explicit
exp [-D x 2
+exp
form for G(z,t)
(s)}
sin nT’~~
‘).
(A’12)
($j2
to t, one finally obtains
equation
(15).
REFERENCES in the modeling of early diagnosis, Earth Planet. Sci. Lett. 29, 33-340, (1976). N.L. Guinasso and D.R. Schink, Quantitative estimates of biological mixing rates in abyssal sediments, J. Geophys. Rev. 80, 3032-3043, (1975). J.B. Fisher, W.J. Lick, P.L. McCall and J.A. Robbins, Vertical mixing of lake sediments by turbificid oligochsetes, J. Geophys. Res. 85, 3997-4006, (1980). B.P. Boudreau and D.M. Imboden, Mathematics of tracer mixing in sediments: III. The theory of nonlocal mixing within sediments, Am. Jour. Sci. 287, 693-719, (1987). C.W. Gardiner, Handbook of Stochastic Methods, Springer Verlag, Berlin, (1983). P.-Z. Wong, The statistical physics of sedimentary rock, Physics Today, 24-32, (December 1988). P. Friedlingstein, Sedimentation dans un milieu complexe: Modelisation mathematique et application iL l’interpretation des don&es geologiques et climatiques. Memoire de Licence, Faculte des Sciences Appliqu&s, Universite Libre de Bruxelles. (1989). B. Mandelbrot, fiactals, Form, Chance and Dimension, Freeman, San Francisco, (1977). W. Feller, An Introductzon to Probability Theory, Vol. I, Wiley, New York, (1968). B.P. Boudreau, Mathematics of tracer mixing in sediments I. Spatially-dependent, diffusive mixing, Am. Jour. Sci. 286, 161-198, (1986). K. Seki, V. Balakrishnan and G. Nicolis, Sensitivity to initial conditions in stochastic systems, Phys. Rev. E47, 155-163, (1993). G. Nicolis, C. Nicolis and K. Seki, Sensitivity to initial conditions in complex systems, In Statistical Physics and Thermodynamics of Nonlinear Nonequilibrium Systems, (Edited by W. Ebeling and W. Muschik), pp. 23-46, (1993). H. Schuster, Determznzstic Chaos, V.C.H. Verlag, Weinheim, (1988). R.E. Plotnik, A fractal model for the distribution of stratigraphic hiatuses, Jour. Geology 94, 885-890, (1986). B. Sapoval, M. Rosso and J.F. Gouyet, The fractal nature of a diffusing front and the relation to percolation, J. Phys. Lett. 48, L149-L156, (1985). J. Feder, Fractals, Plenum, New York. (1988). E.R. Christensen and P.K. Bhunia, Modeling radiotracers in sediments: Comparison with observations in lakes Huron and Michigan, J. Geophys. Res. 91, 855998571, (1986).
1. R.A. Berner, Inclusion of adsorption 2.
3. 4. 5. 6. 7.
8. 9. 10. 11. 12.
13. 14. 15. 16. 17.