CC BY
11
RUSSIAN JOURNAL OF EARTH SCIENCES, English Translation, VOL 1, NO 2, JANUARY 1999
Russian Edition: DECEMBER 1998
Modulated thermoconvective waves in the Earth’s lithosphere
B. I. Birger
Schmidt United Institute of Physics of the Earth, Russian Academy of Sciences, Bol’shaya Gruzinskaya ul. 10, Moscow, 123810 Russia
Abstract. For flows associated with small strains, the rheology of rocks is described by a linear integral (having a memory) law, which reduces to the Andrade law in the case of constant stresses. The continental lithosphere with such a rheology is overstable. Thermoconvective waves propagating through the lithosphere without attenuation have a period of about 200 Ma and a wavelength of the order of 400 km. A pointwise perturbation of the initial temperature in the lithosphere excites amplitude-modulated thermoconvective waves (wave packets). When the initial perturbation occupies a finite area, thermoconvective waves move outside from this area and thermoconvective oscillations (standing waves) are settled within the area. Thermoconvective waves induce oscillations of the Earth’ surface, accompanied by sedimentation and erosion, and can be viewed as a mechanism for the distribution of sediments on continental cratons.
1. Introduction
A power-law non-Newtonian fluid is usually assumed to model slow flows in the mantle and, in particular, convective flows. However,the power-law fluid has no memory in contrast to a real material. A new nonlinear model with a memory was recently proposed recently by Birger [1998]. The proposed model reduces to the power-law fluid model for stationary flows and to the Andrade model for flows associated with small strains.
The steady-state convection beneath continents was studied by Fleitout and Yuen [1984], who used a power-law fluid model and obtained a cold immobile boundary layer (the continental lithosphere). In stability analysis of this layer, the Andrade model must, however, be used. The analysis shows that the lithosphere is overstable, with a period of oscillations of about 200 Ma. These thermoconvective oscillations of the lithosphere are suggested to provide a mechanism for the formation and evolution of sedimentary basins on continental cratons [Birger, 1998]. The vertical crustal movement in sedimentary basins can be respresented as a slow subsidence on which small-amplitude oscillations are superimposed. The longest period of the oscillatory crustal movement is of the same order of magnitude as the pe-
©1998 Russian Journal of Earth Sciences.
Paper No. TJE98006.
Online version of this paper was published on February 15 1999. URL: http://eos.wdcb.rssi.ru/tjes/TJE98006/TJE98006.htm
riod of convective oscillation of the lithosphere found in the stability analysis. Taking into account the difference between the depositional and erosional transport rates we can explain the permanent subsidence of sedimentary basins, as well as their oscillation [Birger, 1998].
Analysis of convective stability for a horizontal layer involves perturbations harmonically depending of the horizontal coordinate. For a layer with the Andrade rheology, the instability is oscillatory and perturbations can take the form of travelling and standing thermoconvective waves. In this study, we adress the problems on the generation of thermoconvective waves under various initial conditions: the linearized equations for thermal convection in a layer with the Andrade rheology are solved for given initial perturbations of temperature.
2. Governing equations
A linear rheological model (having a memory) of the lithosphere is described by the integral relationship
t
2eij = j I<(t - t1)Tij(t1)dt1, (1)
o
where eij and Tij are the deviatoric strain and stress tensors, respectively, t is time, and K(t) is the creep kernel given by
K(t) = t~2!3/?jA, (2)
where A is the Andrade rheological parameter. The
creep kernel (2) is introduced so that, in the case of
constant stress, the strain depends on time as t1!3 (the Andrade law).
The linearized equations of thermal convection in a horizontal layer heated from below are written as
-dp/dx + drxx/dx + drXy/dy + drxz/dz = 0,
-dp/dy + drxy/dx + dTyy/dy + dryz/dz = 0, -dp/dz + drxz/dx + dryz/dy + drzz/dz + Rad = 0, (3)
dvx/dx + dvy/dy + dvz/dz = 0,
dO/dt — d26/dx2 — d26/dy2 — d26/dz2 — vz = 0,
where z is the vertical coordinate, x and y are the horizontal coordinates, v is the velocity, 9 and p are the perturbations of temperature and pressure. The set of equations (3) is written in the dimensionless form. The length scale is layer thickness d and the temperature scale is a temperature drop AT between the hot lower and cold upper surfaces of the layer (both surfaces are supposed to be isothermal). The time scale is d2/n, where k is the thermal diffusivity, and the velocity scale is n/dt. For a Newtonian fluid, the stress (and pressure) scale nr/d2 is usually taken, and the Rayleigh number is Ra = pgaATd3/r/K, where p is the density, a is the thermal expansion coefficient, g is the gravitational acceleration, and T) is the Newtonian viscosity having the dimension of Pa s. For the Andrade rheological medium (the Andrade parameter A has the dimension of Pa s1/3), we introduce a reference viscosity
T]A = A{d2/n)2/3.
Then, the Rayleigh number is defined as
Ra = pgaATd3/t]ak = pgaATd(d2 / k)1^3/A, (4)
and the stress scale is nr/Ad2 = A(d2/k)-1/3.
The lithosphere is characterized by the following depth-averaged values [Birger, 1995]:
d= 2 105 m, AT = 103 °K, a = 4 10-5 °K~1,
p = 3.5 103 kg m-3, g = 10 m s-2, (5)
k = 10-6 m2 s-1, A = 1012 Pa s1/3.
Equations (3) were written in the Boussinesq approximation, which is valid if several dimensionless parameters are small; one of these parameters is a AT. This parameter is estimated for the lithosphere as «AT Pd 0.04.
As the zero-order approximation in the small parameter a AT, the upper free-deformable surface of the layer, where stresses vanish, behaves like a “free” boundary; i.e., the condition of zero normal stress is replaced by a condition of zero vertical velocity on this boundary. Let us suppose the lower boundary of the layer to be also “free”. Then, the boundary conditions on the upper and lower surfaces of the layer are
z = 0, z = 1; vz = txz = Tyz = 6 = 0. (6)
Note that the use (6) permits us to find an exact solution that is not significantly different from the numerical solutions obtained for more realistic boundary conditions [Birger, 1995, 1998].
The set of equations (l)-(6) has a solution in the form of a thermoconvective wave
9 = CE sin ttz, E = exp i(kxx + kyy + u>t),
vx = CE[inkxYi&/(k2 + 7t2)2F(uj)] cos ttz,
k2 = k2x + k2y (7)
vy = CE[iTrkyRsi/(k2 + 7t2)2F(uj)] cos ttz, vz = CE[Ra k2/(k2 + it2)2F(ui)\ sin7rz,
where C is an arbitrary complex factor (the amplitude of temperature), kx and ky are the components of wave vector describing the periodicity in the horizontal directions, k is the wave number, u> is the complex frequency (its imaginary part describes the wave attenuation), and F(u>) is the complex viscosity defined as
F{u) = l/iuK*{iu), K*{iu>) = j K{t)e~iwtdt, (8)
o
K*(iu>) being the Laplace transform of the creep kernel. The complex viscosity is related to the wave number k by the dispersion relation
iojF{bj) (k2 + it2)2 + F{oj) (k2 + 7r2)3 — Ra k2 = 0. (9)
This allows to find such a value of the Rayleigh number Ram (called the minimal critical Rayleigh number) that only a wave with k = km and frequency u> = uim does not attenuate. For the Andrade rheological model, the complex viscosity is
F(w) = (l/3)r(l/3)(iw)-2/3, (l/3)r(l/3) w 1, (10)
where r(a:) is the gamma-function, and Ram, km and uim take on the following values:
km = i/3tt/2 Prf 2.7, cum = iV3(k^ + 7r2) = 7V37t2/4 Pd 30, (11)
Ram = 2 3~1/3{k2m + Tr2)7l3/k2m n 150.
According to the estimated parameters the lithosphere (5), the Rayleigh number Ra for the lithosphere is on the same order of magnitude as Ram. Thus, the lithosphere is in the state close to its instability threshold. If Ra > Ram, the initial perturbations increase with time, and at large t, both the linearized equations of thermal convection (3) and the linear rheological relationship (1) cannot be used. If Ra < Ram, the solution of the set of equations (1)—(6), for a given but not too great initial perturbation of temperature, completely describes the evolution of the perturbations in the layer modeling the lithosphere.
The solution (7)—(11) is obtained in the zero-order approximation in the small parameter a AT. In this approximation, the vertical displacement uz of the upper surface of the layer is equal to zero. In the the same approximation in a A T, we find
uz = aAT[7r(7r2 + 3 k2)/(k2 + 7t2)2]CE, (12)
which does not depend on the rheology of the layer.
If the initial perturbation of temperature at t = 0 is given in the form
90(x, y, z) = 6o exp[i(kxx + kyy)] sin nz, (13)
the factor C in (7) and (12) is defined as C = 9o- However, the initial condition (13) does not completely determine the evolution of the perturbations. Since the wavenumber k, rather than the components kx and ky of the wave vector, enters the dispersion relation (9), in addition to the thermoconvective wave (7), there is a solution with the wave vector {—kx, —ky) describing the wave that runs in the opposite direction. The superposition of the waves, travelling in the opposite directions, forms a standing wave (thermoconvective oscillation). The solutions of governing equations (3)-(6) in the form of both running and standing waves satisfy the initial condition (13). This ambiguity is removed if the initial perturbation (13) is replaced by a more realistic initial perturbation that occupies a finite region and tends to zero for large x and y. Such centered initial perturbations are treated in the next sections of the paper.
In this case, the initial temperature perturbation will be given as
60(x, y, z) = 60(x, y) sin 7TZ,
where 9o(x, y) is not zero only in a limited area on the plane xy, and the solution is sought for in the form
9(x, y, z, t) = 9(x, y, t) sin 7TZ.
Thus, the dependence of the solution on z is fixed, and hence, the problem of three-dimensional distribution of temperature perturbation in the lithosphere is reduced to a two-dimensional problem. In the case, when the initial temperature perturbation does not depend on y and has the form 9o(x), the problem of two- dimensional temperature distribution is reduced to a onedimensional problem.
3. Transient process
The solution in the form of thermoconvective wave is valid for a sufficiently large time elapsing from the moment of the perturbation onset. Under this condition, the stresses harmonic in time induce strains also harmonic in time in the rheological model (1), and the complex analog of viscosity (10) depends only on the frequency, rather than on time.
When the initial temperature perturbation is given in the form (13), we seek a solution of convection equations in the form
9(x, y, z, t) = 9(t) exp[i(kxx + kyy)\ sin nz.
Thus, the coordinate of the temperature and velocity remains in the same form as (7) and the problem reduces to the determination of the time dependence. We seek the solution for k = km and Ra = Ram, whose the values are found in (11).
THe Laplace transformation of the basic equations (1) and (3) and the elimination of all of the physical variables, except the temperature, yields the Laplace transform of the desired function 9(t)
9*(s)
= 90/[s - 2 i-^ikl + n2)1'3^3 + (k2m + 7T2)]. (14) The viscosity analog is now the function
F(s) = s~2!3. (15)
Rearranging the denominator of the expression in the right-hand side of (14),
s - 2 3~ll3{k2m + 7r2)1/3s2/3 + {k2m + 7T2)
= (*m + 7r2)[«l/3 - (31/6/2)(V3 + ¿)] (16)
x [S;/3 - (3i/,6/2)(V3 - *-)][S;/3+3-1/3],
where Sl = s/(fc„ + 7T2).
Since
[(31/6/2)(31/2 ± ¿)]3 = ±*'31/2, (-3“1/3)3 = -1/3,
we may assume that the expression (16) is zero when si = *31 /2, si = —ii1!2 and si = —1/3. However, expressions (10) and (15) for the complex viscosity imply that only the first value of the root, for which the argument of the complex number s satisfies the condition 0 < arg s < 27r, must be used. Then,
(¿31/2)1/3 = (31/®/2)(31/2 + i), (-¿31/2)1/3 = 31/6i,
(—1/3)1/3 = 3“1/3(1 + ¿31/2)/2,
and hence, (16) is zero for si = ii1!2 but not for si = —ii1!2 and si = —1/3. The Laplace transform (14) can be rewritten as
9*{s) = --------------------
[s - lV^{it2m + 7T2)]
(17)
x [s^/3 + (i73)1/3si/3+(i73)2/3]
X [sl/3 - (31/6/2)(V3 - *')](sl/3 + 3-!/3)'
This expression has two singularities: at the branch point s = 0 and at the pole s = ¿31^2{k2rn + 7r2). In the vicinity of point s = 0,
r(s)
(18)
= [<V(*m + ^2)][1 + 2 3~1/3{k2m + tt2)-2/3s2/3 + ...], and in the vicinity of point s = ¿31^2{k2rn + 7r2),
9*(s)
= (3V3/7){V3-2i)90/[s-iV3{k2m+7T2)] +... . (19)
Using the theorem on the asymptotic behavior of Laplace originals [e.g., Von Doetsch, 1967], we find the following asymptotic solution for large t
where u!m = 31l2{k2m + 7r2), which corresponds to (11). The first aperiodic term in the right-hand side of (20) is much smaller than the second, periodic term, even for t equal to the oscillation period 2ir/uim Pd 1/5. Thus, for t > 2ir/uim, the solution takes the form
9{t) = C\9o exp iu>mt, Ci = (3V3/7)(V3 - 2i). (21)
Since ICi I Pd 2, the amplitude of steady-state thermo-convective wave is two times greater than the initial amplitude of temperature perturbation 9o- The argument of complex number C\ defines a phase difference. In the next sections of the paper, the factor C1, which appears in the transient process, is omitted for brevity.
4. One-dimensional problems with initial conditions
Taking ky = 0 and kx = k, we first consider only one-dimensional perturbations independent of the coordinate y. Expanding the complex frequency into the Taylor series in the neighborhood of k — km,
iui = iuim + iV(k — km) — a(k — km)2 + . . ., (22)
where the coefficient V means the group velocity of a packet of thermoconvective waves. Substituting (22) into the dispersion relation (9), we find the values of the coefficients in (22) for the Andrade model
U = 37T, a = (72 + ¿31/2)/7.
The group velocity V of thermoconvective waves in a medium with the Andrade rheology is slightly lower than the phase velocity u}m/km = 77r/2.
The initial temperature perturbation 90 (x) is represented in the form of the Fourier integral
OO
90(x) = J $(&) exp(ikx)dk, (23)
— OO
where $(&) is the Fourier transform of the initial temperature
OO
$(&) = (l/27r) J 9o(x) exp(—ikx)dx. (24)
— OO
The solution satisfying the initial condition is represented as
= [2/31/,3r(—2/3)]6)0/[t(fcm + 7T2)]5/3 + (3V3/7)(V3-2i )90expiuimt, (20)
OO
9(x,t) = J $(fc)exp[i(kx + w(k)t)\dk. (25)
— OO
Substituting the expansion of frequency (22) into (25) and denoting k — km = «, we obtain
9(x,t)=EmJ $(km +u)exp[f(u)t\du, (26)
— OO
where the integration is actually taken over a small vicinity of the point u = 0 (k = km) and the following notation is introduced
Em = exp [i(kmx + u}mt)\,
B = (x + Vt)/t, (27)
f(u) = iBu — au2.
The integral in (26) can then be evaluated by the saddle point method [e.g., Copson, 1965]. The stationary point «0 is found from the condition
/>o) = 0. (28)
Equations (27) and (28) give
uq = iB/2a. (29)
In the vicinity of point uo, function f(u) is represented by the series
f(u) = f(u0) + f(u0)(u - u0)/1!
+ f"{uo){u-u0)2/2\ +... . (30)
Since f'(uo) = 0 and f"(uo) = —2a, (30) takes the form
/(«) = f(uo) ~ a(u - u0)2, (31)
where only two first terms of expansion are retained. Substituting (31) into (26),
9(x,t) (32)
= <!>(km + uo)Em exp[f(u0)t] J exp[-a(u - u0)2t\du.
c
The point uo is the saddle point for the complex function f(u). According to the saddle point method, the path of integration on the complex plane is chosen in such a way that it passes through point u and, in the small vicinity of this point, the path is a segment of the staight line on which the function f(u) — f(uo) is real and negative. In the vicinity of uo,
f(u) - f(u0)
= —a(u — u0)2 = — \a\r2 exp[i(a + 2/3)], (33)
where
a = |a| exp(ia), u — uo = rexp(i/3).
It is clear from (33) that the straight line for which ¡3 = —a/2 must be chosen as the path of integration. On this straight line, the function f(u) —/(«o) takes on real negative values and the integral in (32) reduces to the simple expression
J exp[—a(u — uo)2t]du c
OO
= exp(—¿a/2) / exp(—|a|r2i)(ir
— OO
= exp(—ia/2)(Tr/\a\t)1/2 (34)
= (n/atf2 = [-2 V/"(«o)*]1/2.
In transformation (34), the integral is reduced to the real error integral, and the transition to the infinite limits of integration can be made under the condition \a\t 1. Since |a| Pd 10, the result given by (34) is holds true for t ^¡> 1/10.
Thus, the solution in the case of an arbitrary initial perturbation has the form
9(x,t) = (tt/at)1l2^{km + uo)
x exp + kmx) - ■ (35)
This is a wave packet moving to the left. To obtain the total solution, a term corresponding to the wavenumber k = —km must be added to the right-hand side of (35). This term, which is readily obtained by changing the sign of km and V in the right-hand side of (35), describes a wave packet moving to the right.
The quantity £ = x + Vt can be interpreted as a coordinate in the reference frame that moves together with the wave packet at the group velocity V. The right-hand side of (35) goes to zero at ^2/A\a\t 1. The width of the wave packet can be considered to be of the order of 2(2|a|i)1/2; i.e., the wave packet exists only for |£| < (2|a|i)1/2. Under this condition, it follows from (29) that
l«o| = |£|/2|a|i < (2|a|i)1/2.
Since \a\t 1, we find that |«o| km. Consequently, in the expression <&(km + Mo) in (35) can be replaced with $(fcm).
Figure 1. Propagation of thermoconvective wave packets.
In the case when the initial perturbation of temper- This solution represents two wave packets running in
ature takes place in a fixed point x = 0, function 9o(x) and its Fourior transform are written as
0o(x) = QoS(x), <&{k) = 0o/27r, (36)
where S(x) is the delta-function satisfying the relationship
the opposite directions from the point x = 0, where the initial perturbation takes place. The distributions of temperature at fixed moments of time are shown in Figure 1.
Let the initial perturbation be represented by the sum of two delta functions
OO
2ttS(x) = J exp(ikx)dk.
— OO
Note that (36) retains its form in the dimensional variables since the delta-function has the dimension of reciprocal length and the quantity 0o has the dimension of temperature multiplied by length.
As follows from (35), the asymptotic (\a\t 1) solu-
tion under the initial condition (36) takes the form
6(x, t)
= (0o/27r)(7r/crf)1/2exp ^+ kmx) -+ (0o/27r)(7r/crf)1/2exp ^i(umt - kmx) - ^ 4q^ ^ •
0o(x) = 0o[i(* + /) + S(x — /)].
(37)
The Fourier transform of (37) is
$(fc) = (0O/7r) cos kl = (©o/27r)[exp(ikl) + exp(—ikl)],
but it cannot be used in the problem with the initial condition (37) because its solution is a simple superposition of two solutions obtained for the initial perturbations specified at points x = I and x = — I. Four wave packets move away from these points: two of them move From x = I toward the left and the other two move from x = —I toward the right. Their velocities are iV, respectively. At the moment to = l/V, the centers of the packets, moving in the opposite directions, meet each other at point x = 0. The perturbation of temperature
at — I < x < I is the superposition of two wave packets 0(x,t) = (Qo/2k)(k / at)1!2
x ^exp[ikm(x — I) + iujmt — (x — I + Vt)2/4at\ (38)
+ exp[—ikm(x + /) + iujmt — (x + I — Vt)2/4at]J .
Denoting t\ = t — to (ti is positive after the meeting and negative before it), substituting t = to -Mi into the right-hand side of (38) and observing that
x — I + Vt = x + Vti, x + I — Vt = x — Vti, after algebraic transformations, we obtain
0(x,t) = (Qo/2tt)(tt/at)1!2 exp[—ikml
— (x2 + V2t\)/4ai] ^2 exp(iu}mt) cos kmx + [exp(Vrii*/2ai) — 1] exp(—ikmx + iuimt) (39)
+ [exp(—Vt\xj‘lat) — 1] exp(ikmx + iuimtŸj ,
—l<x<l.
Thus, (39) represents the perturbation of temperature as the superposition of a standing wave and two running waves moving in the opposite directions. The amplitudes of the running waves is zero for small x and t\, and in the vicinity of x = 0, the solution (39) reduces to
6(x, t)
= (©o/^^/Grf)1/2 exp(—ikml) exp(iojmt) cos kmx (40)
Equation (40) holds for \x\ < A* and |ii| < At, where 2Ax is the width of the standing wave zone, 2At is the lifetime of the standing wave. These quantities are estimated as
envelops two wavelengths, and two periods of oscillation occur for the lifetime of the standing wave.
Consider another example of the initial perturbation
90(x) = 0o if |*| < I,
90(x) = 0 if |æ| > I. The Fourier transform for function (42) is
(42)
$(fc) = {do/Trk)suikl
= (6o/2Trki)[exp(ikl) — exp(—ikl)]. (43)
With transform (43), the solution is
6(x, t)
(OO
J (1/k) exp[ik(x + I) + uit\dk (44)
— / (1/k) exp[ik(x — I) + u)t\dk
Since the dominant contribution to the integral in (44) comes from small vicinities of points k = km and k = —km, (44) can be rewritten in the form
= \ao/2nkmi) / exp[ik(x + I) + iuit\dk Voo
OO
— J exp[—ik(x + /) + iuit\dk
— OO OO
— J exp[ik(x — I) + iuit]dk
(45)
2A*
(2|a|i0)1/2, 2Atn(2\a\t0)1/2/V, (41)
exp [—ik(x — I) + iu)t\dk
where |a| Pd 10. The typical dimension of a craton is on the order of 2000 km, which corresponds to I Pd 10.
Then, to = l/V & 1, the width of standing wave zone where u> is a function of k and ui(km) = u>m. Each of the
is 2A* Pd 5, and its lifetime is 2At Pd 1/2. Since the four integrals in (45) describes a wave packet associated
wavelength is 2ir/km Pd 2 and the period of convective with the initial perturbation in the form of ¿-function,
oscillation is 2ir/uim Pri 1/5, the zone of standing wave The first packet moves from point x = —I toward the
Ф(Рг7
Jm I /
-00 — 00
Ят + ^) ехр[/(и, v)t\du dv, (50)
left, the second one moves from x = —I toward the right, the third one moves from x = I toward the left, and the =
fourth one moves from x = I toward the right. The first and fourth packets move outward from the initial wjlere
region of temperature perturbation, and the second and
third packets move inward this initial region. When M ^ ^ v _ ^ ^
they meet, the standing wave is formed in the vicinity
of x = 0, like in the problem with the initial condition Em = exp(ipmx + iqmy + iu>mt),
given by the sum of two delta functions. ti \ ■ 2 0 2
° J j[u,v) = iBiu + iB2v — a\u — Za2uv — a3v ,
5. Two-dimensional problems with = (* + ^)A, B2 = (y + V2t)/t.
initial conditions The stationary point (wo,i>o) °f the function f(u,v) is
defined by the condition
Considering two-dimensional perturbations, we introduce, for brevity, the notation kx = p and ky = q. Then, df/du = df/dv = 0,
9 9,9, r~ , » , which allows us to find
Pm + <lm = k2m, km = V37T/2. (46)
The complex frequency u> is expanded in the power series «0 = i(a3B1 - a2B2)/2(a1a3 - a|),
where
iui = iuim + iVi(p - Pm) + iV2(q - qm)
- a1(p-pm)2 - 2a2(p-pm)(q - qm)
- a3(? - qm)2 + • • •,
Vi = 2лДрт, V2 = 2л/3 q„
(47)
Since
ai = —л/bi + (32/217Г2) (9 + i\f?j)p\ a2 = (32/217Г2) (9 + iVd)pmqm,
a3 = -уД i + (32/217Г2) (9 + iVS)q„.
v0 = i(a\B2 - a2B±)/2{a±a3 - a2). (51)
d2f/du2 = —2ai, d2f/dudv = -2 a2,
d2 f/dv2 = —2a3,
the function in the neighborhood of point («о, ^o) is represented as
f(u, v) = f(u0,vo) - ai(w - Uq)2
(48) - 2a2(u - u0)(v - v0) - a3(v - v0)2, (52)
where
The solution, satisfying the initial conditions, is writ- f(u0,v0) = ~(a3Bf - 2a2B1B2 + axB2)/4(aia3 -
ten as
6(x,y,t)
OO OO
J J ${p,q)exp[i(px + qy + ut)\dp dq, (49)
Substituting expansion (52) into (50),
9{x,y,t)
= Em$(p m Uq , qm + v0) M exp[/(u0, V0)t], (53)
where <&(p, q) is the Fourier transform of the initial tem- where perature distribution $o(x,y)
M = I / exp[—tai (u — Mo)2
OO OO —00—00
®{p,q) = { 1/4tt2) J J 9o(x,y)exp[-i(px + qy)\dx dy. -2ta2(u - u0)(v - v0) - ta3(v - v0)2]du dv.
Substituting (47) for и) into (49), we get 6(x,y,t)
Using the substitution
x = v - v0, у = и - u0 + (a2/ai)(v - v0),
a
M is represented in the form
M = M1M2,
OO
Mi = / exp
= (0o/47r2)(7r/ia4) exp(iojmt)I(x, y,t), I{x,y,t)
(58)
-t[a3-^-)x2
OO
M2 = j exp(-ta1y2)dy.
dx,
27r
= J exp yipmx + iqmy - [a3(x + Vit)
0
- 2a2(x + Vit)(y + V2t) + ai(y + V2t)2] / dip,
Evaluating the integrals Mi and M2 by the saddle point where the dependence of integrand on Lp is defined by
equations (48) and (55). Since it is clear that 9 depends on r rather than on ip, for the case of the ini-
method, we reduce (53) to
9{x,y,t)
= Em$(p
m + ^07 Qm
+ vo)(Tr/a4t)exp[f(uo,vo)t\, (54)
where
a4 = (aia3 — a^)1!2, a2 = 3(1 — ¿24 31/2)/7.
The wave vector components pm and qm, satisfying condition (46), can be represented as
tial point-concentrated perturbation, we can take tp = 0 (and hence x = r, y = 0) in the integrand. Then,
I(x,y,t)
= I(r,t) = exp[(iV3r2 + i9y/3TT2t2)/4a1t\ Ii(r,t),
Pm = km COSip, qm = km sin (p
(55)
where Lp varies from 0 to 2ir. Since all the directions of the wave vector (pm, qm) are equivalent, the integration over Lp is implied in the right-hand sides of (53) and (54)'
If the initial perturbation takes place at point (0, 0), then
90(x, y) = Q08{x)8{y),
and we must substitute in (54) the Fourier transform of this function
Z7T
I1(r,t) = J exp[fi(<p)r]d<p,
0
hif) = + (3V37r/2a|)] cos ip
— (rbi/Aa\t) sin2 Lp = b2[—37rcos Lp + (r/2t) sin2 Lp\,
(59)
where
$(p, q) = 0o/4tt2.
(56)
b1 = (9 + *31/2)8/7, b2 = (4/3)(9 - ¿31 31/2)/247
The integral Ii (r, t) at large r is calculated using the It is convenient to introduce the polar coordinates r and gajjjg p0int method
i)
h(r,t) = [—2tt//"(Lp0)r]1/2 exp[/i(^0)r], (60)
x = r cosip, y = rsin^, 0 < ip < 2ir. (57)
As follows from (55) and (57), the factor Em in (54) Lp0 = tt, fi{<po) = — i[fcm + (3V37r/2a2)] = 3Trb2, (61) becomes
Em = exp (ipmX + iqmy +
= exp[ikmr cos(Lp — ip) + iu>mt\.
fi(ipo) = i[km + {3V3TT/2a24)]-rb1/2a24t = b2[(r/t)-3ir].
After simple algebraic transformations, we find an Thus, the temperature perturbation 9 induced by the asymptotic solution (valid only for sufficiently large t initial point-concentrated perturbation is given by the and r) of the problem with the initial condition speci-equation fied at the given point:
9(x,y,t)
9(x,y,t)
= (0o/47ra4)[-27r/&2^f’]1/2 exp[i(ojmt - kmr) (62)
- (i2/4at)],
where
( = r-3irt, a = (72 + *'31/2)/7.
Note that parameter a has already appeared in the onedimensional problem. Solution (62) describes a cylindrical wave because the dependence on coordinates reduces to the dependence only on r = (x2 + y2)1!2.
The function fi(ip) introduced in (59)—(61), in addition to the saddle point ipo = tt, has the second saddle point ipo = 0. The use of the latter would lead to an additional term in the right-hand side of (60). This term describes a wave running to the point r = 0 from outside and includes the factor
exp[i(umt + kmr)\ exp[-(r + 37rf)/4crf],
which becomes zero for sufficiently large positive r and t.
In the center of the wave packet, i.e., at £ = 0(r = 37rf), equation (62) violates. When £ = 0, the second derivative f”(ipo) is zero (the third derivative is also zero). In this case, the saddle point method leads to
h(r,t) = (l/2)Y(llA)l[-fIV(Vo)r/2A]1l\ (63)
instead of (60). Note that T(l/4) Pd 4. The derivation of the asymptotic relation (63) that holds for large r, is analogous to that of (30)—(34), with the exception that the Taylor series is trancated at the term including the fourth derivative, and furthermore, the integral
OO
J exp(—xi)dx = (l/4)r,l/4)
0
is used instead of the error integral.
Substituting the value of the fourth derivative
fIV{(Po) = -9tt&2
into (63), we find the solution valid in the small vicinity of the wave packet center, i.e., at r Pd 37rf,
9(x,y,t)
(64)
= (0o/87ra4i)r(l/4) (8/37r&2f’)1/4exp[i(wmi - kmr)].
6. Simple and composite sources of thermoconvective waves
A simplified approach of this section to the problems under consideration is not to expand u> in the power series but to set uj = u!m for values of k = (p2 + q2)1!2 close to km in equation (49). Then,
9(x,y,t)
<&(p, q) exp[i(px + qy + ujt)\dp dq
Pd km Ak exp(iuimt) (65)
27r
X J$(km cos (p, km sin Lp) exp[f(<p, i>)r\d<p,
0
where
f(<p, Ip) = ikm cos [ip - -Ip).
Here the polar coordinates are used and the integral is taken over a ring on the plane p, q. The radius of the ring is equal to km, but its thickness Ak is indeterminate in the framework of the simplified approach used now.
In the interval 0 < ip < 2n, function f(tp, ip) has the
following stationary points
fo = TT + V’, if 0 < Ip < 7T,
ipo = —7r + Ip, if 7T < Ip < 2lr, (66)
ipo = ip.
For the initial perturbation in the form of delta function (an elementary radiator of thermoconvective waves), substituting $ = ©o/47r2 into (65), using the saddle point method, and omitting the solution corresponding to the saddle point ipo = ip (a wave arriving at point r = 0 from outside), we find the solution for large r
9{r,t)
(67)
= (Qo/4:TT2)kmAk(2TTi/kmr)1!2 exp[i(uimt - kmr)].
In contrast to solution (62), which takes into account the wavenumber dependence of u> (i.e., dispersion), equation (67) describes a wave with nonmodulated amplitude. It is noteworthy that the simplified approach gives the factor r-1/2 in (67), which determines attenuation of
any cylindrical wave [e.g., Whitham, 1974]. Comparison of solution (67) with solutions (62) and (64) shows that, in order to take the dispersion into account, it is enough to substitute into (67) the following expressions for Ak
Ak = (TT/kma4)(ikm/b2i,t)1/2 exp(-(2/Aat), £ ^ 0,
(68)
Ak = (TT/2krna4t)T(l/A)(-2kllr/?^3b2)1l\ £ = 0.
We then use the same simplified approach for the case when the initial perturbation has the form
9(x, y) = 60, if — I < x < I, —l<y<l,
9(x,y) = 0, if |*| > I, \y\ > I,
that is, the uniform temperature perturbation 90 at the initial moment is given within a square with a side 21. For such initial condition,
$(p,q) = (90/it2) (sinpi sin ql)/pq. (69)
It is convenient to rewrite (69) in the form
®(p,l) = -(0o/47r2)[exp(i>/)
— exp( — ipl)][exp(iql) — exp( — iql)\/pq. (70) Substituting (70) into (65),
9{x,y,t)
27r
= ~(90/4:TT2)kmAkexp(iLjmt) j (exp[ipm(x + I)
0
+ iqm(y + I)] + exp[ipm(x - I) (71)
+iqm(y ~ /)] - exp[ipm{x - I) + iqm(y + /)]
- exp[ipm(x + l) + iq (y-l)f)d<p/p
mim •
Consider separately the integral corresponding to the first term. The vector (x + I, y + /) connects the apex
(—/, —I) of the square with the point (x, y). Introducing
a polar coordinate system whose origin is located at this apex,
x + I = r\ cos tpi, y + I = r\ sin ip 1. (72)
Then,
27T
Ji = J exp[ipm{x + l)+iqm(y + l)\d(p/pmqm.
0
27T
= (2/k^)J exp[ikmri cos(ip — tpi)]dip/ sin2cp (73)
0
Using the saddle point method and omitting the term that describes the wave running to the apex from outside, we obtain the expression for -J\ at sufficiently large
r 1
Ji = (2/fcj^ sin 2ipi)(2'Ki/kmri)1l2 exp(—¿fcmri). (74)
Equation (74) is valid when ip 1 is outside of a small vicinity of points 0, 7r/2, 7r, and 37r/2. Substituting these values of ip 1 into (73), we readily verify that the integral Ji goes to zero for these values of ip\.
The total solution is obtained as the sum of four integrals
9{x,y,t)
4
= -(90/47r2)kmAk^exp[i(uimt (75)
3=1
- kmrj){2/k2m sin2ipj)(2TTi/kmrj)1/2.
Here, the angles ips and ip4 corresponding to the apexes (/,—/) and (—1,1) are measured clockwise, whereas ip\ and ip2 are measured counter-clockwise. Substituting
rj = [(x±l)2 + (yil)2]1/2, ipj = arctan [(y±l)/(x±l)],
l/sin2^j = (1 + tan2^j)/2tan ipj =
= [(* ± I)2 + (y ± /)2]/2(* ±l)(y ±1), into (67), we write the solution in the form
9(x,y,t)
= — (9o/4n2)(Ak/km)(2ni/km)1 /2F(x, y) exp(iwmt), where
F(x,y)
= exp (-ikm[(x + I)2 + (y + /)2]1/2^) [(* + I)2 + (y + l)2]3iy(x + l)(y + l)
+ exp (^-ikm[(x - I)2 + (y - I)2]1/2) [(* - I)2
+ iv - 02]3/4/(* - 0 (y- 0 (76)
G
1=1
0 (
00 1.57 3.14 4.71 6.28 1 = 2
) 0 ( \ V
00 1.57 3.14 4.71 6.28 \"! A. A. A. /
0 c v/\y \/\y Vv V/ w V
00 1.57 3.14 4.71 6.28 \:a a a.j
yy\/\/SJ V/A/\/^ yj\Z\r^ w 1 1 1 1
0.00 1.57 3.14 4.71 6.28
Figure 2. G versus ip for various values of I.
-exp(^-ikm[(x - I)2 + (y + l)2]1!2^ [(x - I)2 + (y + l)2PV(x-l)(y + l)
- exp(-ikm[(x + I)2 + (y - /)2]1/2^ [(x + I)2
+ (y-l)2]3l4/(x + l)(y-l).
Equation (76) is valid for all x and y except the points located on the lines x = ±/, y = ±/ and in the neighborhood of these lines having the width of the order of the wavelength 2ir/km.
Introducing the polar coordinates x = r cos ip, y = r sin ip and using the simple relationship
[(x ± I)2 + (y ± I)2]1!2 = r + Z(± cos ip + sin ip), r^>l,
we obtain that, for r I (the great distance from the initial square), solution (76) reduces to
6(x,y,t) = -(60/4:Tr2)(Ak/km)(2Tri/kmr)1/2x
x G{ip) exp[i(ui r)], (77)
where
G(ip) = [sin(fcmZ cos ip) sin(kml sin ip)\/ cos ip sin ip
Equation (77) describes the wave running outside from the initial square. Note that G and hence the right-hand side of (77) does not depend on ip only under the condition kml <C 1 (the length of the side of the initial square is small in comparison with the thermoconvective wavelength). Under this condition, the initial perturbation given in the square acts like an initial pointwise perturbation.
Function G(ip) is plotted for various values of I in Figure 2. The directivity pattern of the thermoconvective radiation is defined by the amplitude factor |C?(-0) |. As follows from the plots in Figure 2, the directions of the most intense radiation are ip = 0, ip = 7r/2, ip = tt, ip = 37r/2. When I > 2, the real function G(ip) periodically changes its sign, and hence, argG(ip) takes on the values 0 or 7r. Thus, the phase of the wave described by (77), as well as the amplitude, is direction-dependent.
In the theory of radiation, a radiator whose dimensions are much smaller than the radiated wavelength is called simple. The simple radiator has an omnidirectional radiation. A composite radiator is one whose dimensions are not small compared to the wavelength and which emits a directional radiation. Thus, the square where the temperature is initially perturbed is a simple radiator of thermoconvective waves if kml <C 1, but this square is a composite radiator if kml > 1.
In order to account for the dispersion that leads to the amplitude and phase modulation, it is enough to substitute the expressions (68) for Ak into (77).
For r <C I (in the neighborhood of the center of the initial square), the following relationship takes place
[(x±l)2-[-(yztl)2]1!2 = 21/,2/ + (21/,2/2)r(± cos ^±sin ip),
where 21!2l is the distance between the apex and the center of the square. By using this approximate equation, we reduce (76) for r <C I to the form
6{x,y,t)
(78)
= H(l) exp(iujmt) cos[(\f2/2)kmx] cos[(V2/2)kmy\, where
H(l)
= —(Go/ir2)(Ak/km)2314(27n/kml)1 ^2 exp(—i21/2kml)-
As is seen from (78), a standing wave with square cells is formed in the central part (r <C I) of the initial square.
The sides of the cells are parallel to the sides of the initial square and the length of the cell side is 21!2tt / km-The same standing wave is formed in the neighborhood of the point (0,0) in the case when the initial temperature perturbation is given at four points: (1,1), (/,—/), (—1,1), and (—/,—/) i.e. when the initial perturbation square is replaced by the pointwise perturbation at its apexes.
Note that in order to obtain the solution (77), valid for r > i, we could use the expression (69) for $. We represents $ in the form (70) to find the solution for r <C I i.e. in the neighborhood of the center of the initial square. The result (78), obtained by the saddle point method, is holds true under the condition kml 1. (length of the side of the initial square is much greater then thermoconvective wave length). This condition is satisfied when the initial temperature perturbation covers a craton as a whole, and therefore, I Pd 10. for km =2.7.
Discussion
The distribution of temperature in the lithosphere is represented in the form
T(x,y, z,t) = 1 — z + 0(x, y, t) sin 7TZ, 0 < z < 1,
where 6(x, y, t) can be interpreted as a temperature perturbation in the middle of the lithosphere (z = 1/2). The temperature perturbation 9(x,y,t) is found under a few simple initial conditions: the initial perturbations 8o(x,y) are given on the straight line (x = 0), in the strip (—I < x <l), at the point (x = 0, y = 0), at four points, or within a square.
For a plane thermoconvective wave with the wave vector (kx, ky), the displacement of the upper surface of the layer (lithosphere) is related to the temperature perturbation by (12). This relation is easily generalized to the case of a packet of thermoconvective waves with the wave numbers close to km
uz (x, y, t) = aAT[7r(7r2 + 3^4)/(tt2 + k2m)2]6(x, y, t).
When more real boundary conditions on the upper and lower surfaces of the lithosphere are considered [Birger, 1988], the dependence of temperature perturbation on the vertical coordinate z is not determined by the function sin irz and the wave number km is not equal to 2.7. However, the displacement uz(x,y,t) remains equal to the temperature perturbation 9(x,y,t) multiplied by a constant whose value depends on the boundary conditions. This is also valid for the case when the depth dependence of physical parameters (in particular, the Andrade rheological parameter) of the lithosphere is taken into account. Thus, functions 9(x,y,t), found
above under various initial conditions, describe the vertical displacements of upper surface of the lithosphere, with accuracy to a constant factor. These displacements determine sedimentary processes.
The initial temperature point-concentrated perturbation can be treated as a source of thermoconvective waves in the lithosphere. When the initial perturbation occupies a finite area, thermoconvective waves propagate outward from this area and thermoconvective oscillations (standing waves) are settled inside the area. The thermoconvective oscillations create a system of convective cells in the lithosphere. Over the convective cells, the surface of the lithosphere subsides (or rises) forming sedimentary basins. Since the Rayleigh number for the lithosphere is not greater than Ram, the shape of cells and hence the shape of basins is determined by initial perturbations. The initial perturbation in the square considered above leads to the appearance of square cells. However, other initial conditions may lead to rectangular, hexagonal, and other shapes of cells and basins.
Thermoconvective oscillations may be considered as a mechanism inducing the formation and evolution of sedimentary basins on continental cartons [Birger, 1998]. In this paper, we have studied the excitation of thermoconvective waves and oscillations by initial perturbations of temperature. Note that these waves are also generated from initial vertical displacements of the Earth’s surface (relief perturbations).
Packets of thermoconvective waves propagate in the lithosphere with the velocity 37m/d Pd 0.15 cm/year and create sediment-filled depressions on the upper surface of the lithosphere. The age of sediments increases proportionally to the distance from the front of the wave to its source. Thermoconvective waves may be related to the development of peripheral depression on a craton. In this process known in geology [Beloussov, 1978; Kham, 1973], the initial depression occurring in a geologically active area (“geosyncline”) adjacent the craton slowly develops within the craton.
References
Beloussov, V. V., Endogenous Regimes of Continents, Ne-dra, Moscow, 1978 (in Russian).
Birger, B. I., Linear and weakly nonlinear problems of the theory of thermal convection in the earth’s mantle, Phys. Earth Planet. Inter., 50, 92-98, 1988.
Birger, B. I., On the thermoconvective mechanism for oscillatory vertical crustal movements, Phys. Earth Planet. Inter., 92, 279-291, 1995.
Birger, B. I., Rheology of the Earth and thermoconvective mechanism for sedimentary basins formation, Geophys. I. Inter., 134, 1-12, 1998.
Copson, E. T., Asymptotic Expansions, University Press, 170 pp., Cambridge, 1965.
Fleitout, L., and Yuen, D. A., Steady state, secondary convection beneath lithosphere plates with temperature- and pressure-dependent viscosity, J. Geophys. Res., 89, 92279244, 1984.
Khain, V. E., General Geotectonics, 512 pp., Nedra, Moscow, 1973 (in Russian).
Von Doetch, G., Anleitung zum praktischen Gebrauch der
Laplace-Transformation und der Z-Transformation, Springer, 365 pp., Munich, 1967.
Whitham, G. B., Linear and Nonlinear waves, John Wiley and sons, New York, 1974.
(Received November 15 1998.)
