Russian Journal of Earth Sciences
Vol 1, No. 1, July 1998
Translated 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


Contents


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 period 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

eqn001.gif(1)

where eij and tij are the deviatoric strain and stress tensors, respectively, t is time, and K(t) is the creep kernel given by

eqn002.gif(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

eqn003.gif

eqn004.gif

eqn005.gif(3)

eqn006.gif

eqn007.gif

where z is the vertical coordinate, x and y are the horizontal coordinates, v is the velocity, q 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 DT between the hot lower and cold upper surfaces of the layer (both surfaces are supposed to be isothermal). The time scale is d2/k, where k is the thermal diffusivity, and the velocity scale is k/d. For a Newtonian fluid, the stress (and pressure) scale khd2 is usually taken, and the Rayleigh number is Ra = rgaDTd3/hk, where r is the density, a is the thermal expansion coefficient, g is the gravitational acceleration, and h is the Newtonian viscosity having the dimension of Pa s. For the Andrade rheological medium (the Andrade parameter A has the dimension of Pa s 1/3 ), we introduce a reference viscosity

hA = A(d2/k)2/3.

Then, the Rayleigh number is defined as

eqn008.gif(4)

and the stress scale is khA d2 = A(d2/k)-1/3.

The lithosphere is characterized by the following depth-averaged values [Birger, 1995]:

eqn009.gif

eqn010.gif(5)

eqn011.gif

Equations (3) were written in the Boussinesq approximation, which is valid if several dimensionless parameters are small; one of these parameters is aDT. This parameter is estimated for the lithosphere as aDTapprox 0.04. As the zero-order approximation in the small parameter aDT, 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

eqn012.gif(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 (1)-(6) has a solution in the form of a thermoconvective wave

eqn013.gif

eqn014.gif

eqn015.gif(7)

eqn016.gif

eqn017.gif

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, w is the complex frequency (its imaginary part describes the wave attenuation), and F(w) is the complex viscosity defined as

eqn018.gif(8)

Kast(iw) being the Laplace transform of the creep kernel. The complex viscosity is related to the wave number k by the dispersion relation

eqn019.gif(9)

This allows to find such a value of the Rayleigh number Ra m (called the minimal critical Rayleigh number) that only a wave with k = km and frequency w = wm does not attenuate. For the Andrade rheological model, the complex viscosity is

eqn020.gif(10)

where G(x) is the gamma-function, and Ra m, km and wm take on the following values:

eqn021.gif

eqn022.gif(11)

eqn023.gif

According to the estimated parameters the lithosphere (5), the Rayleigh number Ra for the lithosphere is on the same order of magnitude as Ra m. 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 Rale 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 aDT. In this approximation, the vertical displacement uz of the upper surface of the layer is equal to zero. In the the same approximation in aDT, we find

eqn024.gif(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

eqn025.gif(13)

the factor C in (7) and (12) is defined as C = q0. 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

q0(x, y, z) = q0(x, y) sin pz,

where q0(x, y) is not zero only in a limited area on the plane xy, and the solution is sought for in the form

q(x, y, z, t) = q(x, y, t) sinpz.

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 q0(x), the problem of two- dimensional temperature distribution is reduced to a one-dimensional 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

q(x, y, z, t) = q(t) exp[i(kxx + kyy)] sinpz.

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 = Ra m, 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 q(t)

eqn026.gif

eqn027.gif(14)

The viscosity analog is now the function

eqn028.gif(15)

Rearranging the denominator of the expression in the right-hand side of (14),

eqn029.gif

eqn030.gif(16)

eqn031.gif

where s1 = s/(k2m + p2).

Since

[(31/6/2) (31/2pm i)]3 = pm i31/2, (-3-1/3)3 = -1/3,

we may assume that the expression (16) is zero when s1 = i31/2, s1 = -i31/2 and s1 = -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 0le arg s < 2p, must be used. Then,

(i31/2)1/3 = (31/6/2) (31/2 + i), (-i31/2)1/3 = 31/6i,

(-1/3)1/3 = 3-1/3 (1+i31/2)/2,

and hence, (16) is zero for s1 = i31/2 but not for s1 = -i31/2 and s1 = -1/3. The Laplace transform (14) can be rewritten as

eqn032.gif
 (17)
eqn034.gif

This expression has two singularities: at the branch point s = 0 and at the pole s = i31/2(k2m + p2). In the vicinity of point s=0,

qast (s)

 (18)

= [q0 / (k2m + p2)] [1 + 2 3-1/3 (k2m + p2)-2/3 s2/3 + ldots],

and in the vicinity of point s = i31/2(k2m + p2),

eqn036.gif

eqn037.gif(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

eqn038.gif

eqn039.gif(20)

where wm = 31/2(k2m + p2), 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 2p/wmapprox 1/5. Thus, for tge 2p/wm, the solution takes the form

eqn040.gif(21)

Since |C1|approx 2, the amplitude of steady-state thermoconvective wave is two times greater than the initial amplitude of temperature perturbation q0. The argument of complex number C1 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,

eqn041.gif(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

V = 3p, a = (72+i31/2)/7.

The group velocity V of thermoconvective waves in a medium with the Andrade rheology is slightly lower than the phase velocity wm/km = 7p/2.

The initial temperature perturbation q0(x) is represented in the form of the Fourier integral

eqn042.gif(23)

where F(k) is the Fourier transform of the initial temperature

eqn043.gif(24)

The solution satisfying the initial condition is represented as

eqn044.gif(25)

Substituting the expansion of frequency (22) into (25) and denoting k - km = u, we obtain

eqn045.gif(26)

where the integration is actually taken over a small vicinity of the point u = 0 (k = km) and the following notation is introduced

eqn046.gif

eqn047.gif(27)

eqn048.gif

The integral in (26) can then be evaluated by the saddle point method [e.g., Copson, 1965]. The stationary point u0 is found from the condition

eqn049.gif(28)

Equations (27) and (28) give

eqn050.gif(29)

In the vicinity of point u0, function f(u) is represented by the series

eqn051.gif

eqn052.gif(30)

Since fprime(u0) = 0 and fprimeprime(u0) = -2a, (30) takes the form

eqn053.gif(31)

where only two first terms of expansion are retained.

Substituting (31) into (26),

eqn054.gif(32)

eqn055.gif

The point u0 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(u0) is real and negative. In the vicinity of u0,

eqn056.gif

eqn057.gif(33)

where

a = |a| exp(ia), u - u0 = r exp(ib).

It is clear from (33) that the straight line for which b = -a/2 must be chosen as the path of integration. On this straight line, the function f(u)-f(u0) takes on real negative values and the integral in (32) reduces to the simple expression

eqn058.gif

eqn059.gif

eqn060.gif(34)

eqn061.gif

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|tgg 1. Since |a|approx 10, the result given by (34) is holds true for tgg 1/10.

Thus, the solution in the case of an arbitrary initial perturbation has the form

eqn062.gif

eqn063.gif(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 = 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 x2/4|a|tgg 1. The width of the wave packet can be considered to be of the order of 2(2|a|t)1/2 ; i.e., the wave packet exists only for |x|le (2|a|t)1/2. Under this condition, it follows from (29) that

|u0| = |x|/2|a|tle (2|a|t)1/2.

Since |a|tgg 1, we find that |u0|ll km. Consequently, in the expression F(km + u0) in (35) can be replaced with F(km).

In the case when the initial perturbation of temperature takes place in a fixed point x = 0, function q0(x) and its Fourior transform are written as

eqn064.gif(36)

where d(x) is the delta-function satisfying the relationship

eqn065.gif

Note that (36) retains its form in the dimensional variables since the delta-function has the dimension of reciprocal length and the quantity Q0 has the dimension of temperature multiplied by length.

As follows from (35), the asymptotic (|a|tgg 1) solution under the initial condition (36) takes the form

fig01

q (x,t)

eqn066.gif

eqn067.gif

This solution represents two wave packets running in 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

eqn068.gif(37)

The Fourier transform of (37) is

F(k) = (Q0/p) cos kl = (Q0/2p) [ 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=l and x=-l. Four wave packets move away from these points: two of them move From x = l toward the left and the other two move from x = -l toward the right. Their velocities are pm V, respectively. At the moment t0 = l/V, the centers of the packets, moving in the opposite directions, meet each other at point x=0. The perturbation of temperature at -lle xle l is the superposition of two wave packets

eqn069.gif

eqn070.gif(38)

eqn071.gif

Denoting t1 = t - t0t1 is positive after the meeting and negative before it), substituting t = t0 + t1 into the right-hand side of (38) and observing that

x - l + Vt = x + Vt1, x + l - Vt = x - Vt1,

after algebraic transformations, we obtain

eqn072.gif

eqn073.gif

eqn074.gif(39)

eqn075.gif

eqn076.gif

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 t1, and in the vicinity of x = 0, the solution (39) reduces to

eqn077.gif

eqn078.gif(40)

Equation (40) holds for |x|le Dx and |t1|le Dt, where 2Dx is the width of the standing wave zone, 2Dt is the lifetime of the standing wave. These quantities are estimated as

eqn079.gif(41)

where |a|approx 10. The typical dimension of a craton is on the order of 2000 km, which corresponds to lapprox 10. Then, t0 = l/Vapprox 1, the width of standing wave zone is 2Dxapprox 5, and its lifetime is 2Dtapprox 1/2. Since the wavelength is 2p/kmapprox 2 and the period of convective oscillation is 2p/wmapprox 1/5, the zone of standing wave envelops two wavelengths, and two periods of oscillation occur for the lifetime of the standing wave.

Consider another example of the initial perturbation

eqn080.gif

eqn081.gif(42)

The Fourier transform for function (42) is

eqn082.gif

eqn083.gif(43)

With transform (43), the solution is

eqn084.gif

eqn085.gif(44)

eqn086.gif

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

eqn087.gif

eqn088.gif

eqn089.gif(45)

eqn090.gif

eqn091.gif

where w is a function of k and w(km) = wm. Each of the four integrals in (45) describes a wave packet associated with the initial perturbation in the form of d -function. The first packet moves from point x = -l toward the left, the second one moves from x = -l toward the right, the third one moves from x = l toward the left, and the fourth one moves from x = l toward the right. The first and fourth packets move outward from the initial region of temperature perturbation, and the second and third packets move inward this initial region. When they meet, the standing wave is formed in the vicinity of x = 0, like in the problem with the initial condition given by the sum of two delta functions.


5. Two-dimensional problems with initial conditions

Considering two-dimensional perturbations, we introduce, for brevity, the notation kx = p and ky = q. Then,

eqn092.gif(46)

The complex frequency w is expanded in the power series

eqn093.gif

eqn094.gif(47)

eqn095.gif

where

eqn096.gif

eqn097.gif

eqn098.gif(48)

eqn099.gif

The solution, satisfying the initial conditions, is written as

eqn100.gif

eqn101.gif(49)

where F(p, q) is the Fourier transform of the initial temperature distribution q0(x, y)

eqn102.gif

Substituting (47) for w into (49), we get

eqn103.gif

eqn104.gif(50)

where

u = p-pm, v = q-qm,

Em = exp(ipm x + iqm y + iwm t),

f(u,v) = iB1 u + iB2 v - a1 u2 - 2a2 uv - a3 v2,

B1 = (x + V1 t)/t, B2 = (y + V2 t)/t.

The stationary point (u0, v0) of the function f(u, v) is defined by the condition

partial f/partial u = partial f/partial v = 0,

which allows us to find

eqn105.gif

eqn106.gif(51)

Since

partial2 f/partial u2 = -2a1, partial2 f/partial upartial v = -2a2,

partial2 f/partial v2 = -2a3,

the function in the neighborhood of point (u0, v0) is represented as

eqn107.gif

eqn108.gif(52)

where

f(u0, v0) = - (a3 B21 - 2a2 B1 B2 + a1 B22) /4(a1 a3 - a22)

Substituting expansion (52) into (50),

eqn109.gif

eqn110.gif(53)

where

eqn111.gif

- 2ta2 (u-u0)(v-v0) - ta3 (v-v0)2]du dv.

Using the substitution

x = v-v0, y = u-u0 + (a2 /a1)(v-v0),

M is represented in the form

M = M1 M2,

eqn112.gif

eqn113.gif

Evaluating the integrals M1 and M2 by the saddle point method, we reduce (53) to

eqn114.gif

eqn115.gif(54)

where

a4 = (a1 a3 - a22)1/2, a24 = 3(1-i24 31/2)/7.

The wave vector components pm and qm, satisfying condition (46), can be represented as

eqn116.gif(55)

where j varies from 0 to 2p. Since all the directions of the wave vector (pm, qm) are equivalent, the integration over j is implied in the right-hand sides of (53) and (54).

If the initial perturbation takes place at point (0, 0), then

q0(x,y) = Q0d(x) d(y),

and we must substitute in (54) the Fourier transform of this function

eqn117.gif(56)

It is convenient to introduce the polar coordinates r and y

eqn118.gif(57)

As follows from (55) and (57), the factor Em in (54) becomes

Em = exp(ipm x + iqm y + iwm t)

= exp[ikm r cos(j - y) + iwm t].

Thus, the temperature perturbation q induced by the initial point-concentrated perturbation is given by the equation

eqn119.gif

eqn120.gif

eqn121.gif(58)

eqn122.gif

eqn123.gif

where the dependence of integrand on j is defined by equations (48) and (55). Since it is clear that q depends on r rather than on y, for the case of the initial point-concentrated perturbation, we can take y = 0 (and hence x = r, y = 0 ) in the integrand. Then,

eqn124.gif

eqn125.gif

eqn126.gif(59)

eqn127.gif

eqn128.gif

eqn129.gif

where

b1 = (9 + i31/2)8/7, b2 = (4/3)(9-i31 31/2)/247

The integral I1(r, t) at large r is calculated using the saddle point method

eqn130.gif(60)

eqn131.gif(61)

fprimeprime1 (j0) = i[km + (33p/2a24)]-rb1 /2a24 t = b2 [(r/t)-3p].

After simple algebraic transformations, we find an asymptotic solution (valid only for sufficiently large t and r ) of the problem with the initial condition specified at the given point:

eqn132.gif

eqn133.gif(62)

eqn134.gif

where

x = r-3pt, a = (72 + i31/2)/7.

Note that parameter a has already appeared in the one-dimensional 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 f1(j) introduced in (59)-(61), in addition to the saddle point j0 = p, has the second saddle point j0 = 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(wm t + km r)] exp[-(r + 3pt)/4at],

which becomes zero for sufficiently large positive r and t.

In the center of the wave packet, i.e., at x = 0 (r = 3pt), equation (62) violates. When x = 0, the second derivative fprimeprime1(j0) is zero (the third derivative is also zero). In this case, the saddle point method leads to

eqn135.gif(63)

instead of (60). Note that G(1/4)approx 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

eqn136.gif

is used instead of the error integral.

Substituting the value of the fourth derivative

fIV(j0) = -9pb2

into (63), we find the solution valid in the small vicinity of the wave packet center, i.e., at rapprox 3pt,

q (x,y,t)

 (64)

= (Q0 /8pa4 t) G (1/4) (8/3pb2 r)1/4 exp[i(wm t - km r)].


6. Simple and composite sources of thermoconvective waves

A simplified approach of this section to the problems under consideration is not to expand w in the power series but to set w = wm for values of k = (p2 + q2)1/2 close to km in equation (49). Then,

eqn138.gif

eqn139.gif

eqn140.gif(65)

eqn141.gif

where

f(j, y) = ikm cos(j - y).

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 Dk is indeterminate in the framework of the simplified approach used now. In the interval 0le j< 2p, function f(j, y) has the following stationary points

eqn142.gif

eqn143.gif(66)

j0 = y.

For the initial perturbation in the form of delta function (an elementary radiator of thermoconvective waves), substituting F = Q0/4p2 into (65), using the saddle point method, and omitting the solution corresponding to the saddle point j0 = y (a wave arriving at point r = 0 from outside), we find the solution for large r

q (r,t)

 (67)

= (Q0 /4p2) km Dk (2pi/km r)1/2 exp[i(wm t - km r)].

In contrast to solution (62), which takes into account the wavenumber dependence of w (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 Dk

Dk = (p/km a4)(ikm /b2 xt)1/2 exp(- x2 /4at), x 0,

 (68)

Dk = (p/2km a4 t)G(1/4)(-2k2m r/3p3 b2)1/4, x = 0.

We then use the same simplified approach for the case when the initial perturbation has the form

q(x,y) = q0, if -lle x le l, -lle y le l,

q(x,y) = 0, if |x| > l, |y| > l,

that is, the uniform temperature perturbation q0 at the initial moment is given within a square with a side 2l. For such initial condition,

eqn146.gif(69)

It is convenient to rewrite (69) in the form

eqn147.gif

eqn148.gif(70)

Substituting (70) into (65),

q (x,y,t)

eqn149.gif

eqn150.gif(71)

eqn151.gif

eqn152.gif

Consider separately the integral corresponding to the first term. The vector ( x + l, y + l ) connects the apex (-l, -l) of the square with the point (x, y). Introducing a polar coordinate system whose origin is located at this apex,

eqn153.gif(72)

Then,

eqn154.gif

eqn155.gif(73)

Using the saddle point method and omitting the term that describes the wave running to the apex from outside, we obtain the expression for J1 at sufficiently large r1

eqn156.gif(74)

Equation (74) is valid when y1 is outside of a small vicinity of points 0, p/2, p, and 3p/2. Substituting these values of y1 into (73), we readily verify that the integral J1 goes to zero for these values of y1.

The total solution is obtained as the sum of four integrals

eqn157.gif

eqn158.gif(75)

eqn159.gif

Here, the angles y3 and y4 corresponding to the apexes (l, -l) and (-l, l) are measured clockwise, whereas y1 and y2 are measured counter-clockwise. Substituting

rj = [(xpm l)2 + (ypm l)2]1/2, yj = arctan [(ypm l)/(xpm l)],

1/ sin2yj = (1+ tan2 yj)/2 tan yj=

= [(xpm l)2 + (ypm l)2]/2(xpm l)(ypm l),

into (67), we write the solution in the form

q (x,y,t)

= - (q0 /4p2)(Dk/km)(2pi/km)1/2 F (x,y) exp(iwm t),

where

eqn160.gif

eqn161.gif

eqn162.gif

eqn163.gif

eqn164.gif(76)

eqn165.gif

eqn166.gif

eqn167.gif

eqn168.gif

Equation (76) is valid for all x and y except the points located on the lines x = pm l, y = pm l and in the neighborhood of these lines having the width of the order of the wavelength 2p/km.

Introducing the polar coordinates x = r cosy, y = r siny and using the simple relationship

[(xpm l)2 + (ypm l)2]1/2 = r + l(pm cosy + siny), rgg l,

we obtain that, for rgg l (the great distance from the initial square), solution (76) reduces to

q (x,y,t) = - (q0 /4p2)(Dk/km) (2pi/km r)1/2 times

eqn169.gif(77)

where

G(y) = [ sin(km l cosy) sin(km l siny)]/ cosy siny

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 y only under the condition kmlll 1 (the length of the side of the initial square is small in comparison with the thermoconvective wavelength). Under this condition, the initial pertur-bation given in the square acts like an initial pointwise perturbation.

fig02 Function G(y) is plotted for various values of l in Figure 2. The directivity pattern of the thermoconvective radiation is defined by the amplitude factor |G(y)|. As follows from the plots in Figure 2, the directions of the most intense radiation are y = 0, y = p/2, y = p, y = 3p/2. When lge 2, the real function G(y) periodically changes its sign, and hence, arg G(y) takes on the values 0 or p. 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 kmlll 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 Dk into (77).

For rll l (in the neighborhood of the center of the initial square), the following relationship takes place

[(xpm l)2 + (ypm l)2]1/2 = 21/2l + (21/2/2)r (pm cosypm siny),

where 21/2l is the distance between the apex and the center of the square. By using this approximate equation, we reduce (76) for rll l to the form

q (x,y,t)

 (78)

= H (l) exp(iwm t) cos[(2/2)km x] cos[(2/2)km y],

where

H(l)

= -(q0 /p2)(Dk/km)23/4 (2pi/km l)1/2 exp(- i21/2km l).

As is seen from (78), a standing wave with square cells is formed in the central part (rll l) 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/2p/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: (l, l), (l, -l), (-l, l), and (-l, -l) 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 rgg l, we could use the expression (69) for F. We represents F in the form (70) to find the solution for rll l 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 kmlgg 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, lapprox 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 + q(x,y,t) sin pz, 0le zle 1,

where q(x, y, t) can be interpreted as a temperature perturbation in the middle of the lithosphere (z = 1/2). The temperature perturbation q(x, y, t) is found under a few simple initial conditions: the initial perturbations q0(x, y) are given on the straight line (x = 0), in the strip (-lle xle 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) = aDT [p(p2 + 3k2m)/(p2 + k2m)2] q(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 pz and the wave number km is not equal to 2.7. However, the displacement uz(x,y,t) remains equal to the temperature perturbation q(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 q(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 Ra m, 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 3pk/dapprox 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; Khain, 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, Nedra, 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. J. 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, 9227-9244, 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.


 Load files for print and local use.


 
This document was generated by TeXWeb (Win32, v.1.0) on March 28, 1999.