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

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

(1)

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

(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 t^1/3 (the Andrade law).

The linearized equations of thermal convection in a horizontal layer heated from below are written as

(3)

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 d²/k, where k is the thermal diffusivity, and the velocity scale is k/d. For a Newtonian fluid, the stress (and pressure) scale khd² is usually taken, and the Rayleigh number is Ra = rgaDTd³/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

Then, the Rayleigh number is defined as

(4)

and the stress scale is kh_A d² = A(d²/k)^-1/3.

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

(5)

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

(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

(7)

where C is an arbitrary complex factor (the amplitude of temperature), k_x and k_y 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

(8)

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

(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 = k_m and frequency w = w_m does not attenuate. For the Andrade rheological model, the complex viscosity is

(10)

where G(x) is the gamma-function, and Ra _m, k_m and w_m take on the following values:

(11)

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 > Ra_m, 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 Ra_m, 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 u_z of the upper surface of the layer is equal to zero. In the the same approximation in aDT, we find

(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

(13)

the factor C in (7) and (12) is defined as C = q₀. However, the initial condition (13) does not completely determine the evolution of the perturbations. Since the wavenumber k, rather than the components k_x and k_y of the wave vector, enters the dispersion relation (9), in addition to the thermoconvective wave (7), there is a solution with the wave vector (-k_x, -k_y) 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

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

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 q₀(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

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 = k_m 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)

(14)

The viscosity analog is now the function

(15)

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

(16)

where s₁ = s/(k²_m + p²).

Since

[(3^1/6/2) (3^1/2 i)]³ = i3^1/2, (-3^-1/3)³ = -1/3,

we may assume that the expression (16) is zero when s₁ = i3^1/2, s₁ = -i3^1/2 and s₁ = -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 < 2p, must be used. Then,

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

	(17)

This expression has two singularities: at the branch point s = 0 and at the pole s = i3^1/2(k²_m + p²). In the vicinity of point s=0,

q (s)

(18)

= [q₀ / (k²_m + p²)] [1 + 2 3^-1/3 (k²_m + p²)^-2/3 s^2/3 + ],

and in the vicinity of point s = i3^1/2(k²_m + p²),

(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

(20)

where w_m = 3^1/2(k²_m + p²), 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/w_m 1/5. Thus, for t 2p/w_m, the solution takes the form

(21)

Since |C₁| 2, the amplitude of steady-state thermoconvective wave is two times greater than the initial amplitude of temperature perturbation q₀. The argument of complex number C₁ defines a phase difference. In the next sections of the paper, the factor C₁, which appears in the transient process, is omitted for brevity.

4. One-dimensional problems with initial conditions

Taking k_y = 0 and k_x = 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 - k_m,

(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

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

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

(23)

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

(24)

The solution satisfying the initial condition is represented as

(25)

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

(26)

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

(27)

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

(28)

Equations (27) and (28) give

(29)

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

(30)

Since f(u₀) = 0 and f(u₀) = -2a, (30) takes the form

(31)

where only two first terms of expansion are retained.

Substituting (31) into (26),

(32)

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

(33)

where

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(u₀) takes on real negative values and the integral in (32) reduces to the simple expression

(34)

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

(35)

This is a wave packet moving to the left. To obtain the total solution, a term corresponding to the wavenumber k = -k_m must be added to the right-hand side of (35). This term, which is readily obtained by changing the sign of k_m 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 x²/4|a|t 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| (2|a|t)^1/2. Under this condition, it follows from (29) that

|u₀| = |x|/2|a|t (2|a|t)^1/2.

Since |a|t 1, we find that |u₀| k_m. Consequently, in the expression F(k_m + u₀) in (35) can be replaced with F(k_m).

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

(36)

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

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

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

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

(37)

The Fourier transform of (37) is

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 V, respectively. At the moment t₀ = l/V, the centers of the packets, moving in the opposite directions, meet each other at point x=0. The perturbation of temperature at -l x l is the superposition of two wave packets

(38)

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

after algebraic transformations, we obtain

(39)

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

(40)

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

(41)

where |a| 10. The typical dimension of a craton is on the order of 2000 km, which corresponds to l 10. Then, t₀ = l/V 1, the width of standing wave zone is 2Dx 5, and its lifetime is 2Dt 1/2. Since the wavelength is 2p/k_m 2 and the period of convective oscillation is 2p/w_m 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

(42)

The Fourier transform for function (42) is

(43)

With transform (43), the solution is

(44)

Since the dominant contribution to the integral in (44) comes from small vicinities of points k = k_m and k = -k_m, (44) can be rewritten in the form

(45)

where w is a function of k and w(k_m) = w_m. 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 k_x = p and k_y = q. Then,

(46)

The complex frequency w is expanded in the power series

(47)

where

(48)

The solution, satisfying the initial conditions, is written as

(49)

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

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

(50)

where

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

f/ u = f/ v = 0,

which allows us to find

(51)

Since

² f/ u² = -2a₁, ² f/ u v = -2a₂,

² f/ v² = -2a₃,

the function in the neighborhood of point (u₀, v₀) is represented as

(52)

where

Substituting expansion (52) into (50),

(53)

where

Using the substitution

M is represented in the form

Evaluating the integrals M₁ and M₂ by the saddle point method, we reduce (53) to

(54)

where

The wave vector components p_m and q_m, satisfying condition (46), can be represented as

(55)

where j varies from 0 to 2p. Since all the directions of the wave vector (p_m, q_m) 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

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

(56)

It is convenient to introduce the polar coordinates r and y

(57)

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

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

(58)

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,

(59)

where

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

(60)

(61)

f₁ (j₀) = i[k_m + (33p/2a²₄)]-rb₁ /2a²₄ t = b₂ [(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:

(62)

where

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 = (x² + y²)^1/2.

The function f₁(j) introduced in (59)-(61), in addition to the saddle point j₀ = p, has the second saddle point j₀ = 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

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 f₁(j₀) is zero (the third derivative is also zero). In this case, the saddle point method leads to

(63)

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

is used instead of the error integral.

Substituting the value of the fourth derivative

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

(64)

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 = w_m for values of k = (p² + q²)^1/2 close to k_m in equation (49). Then,

(65)

where

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

(66)

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

(67)

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

(68)

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

q(x,y) = q₀, if -l x l, -l y l,

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

(69)

It is convenient to rewrite (69) in the form

(70)

Substituting (70) into (65),

(71)

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,

(72)

Then,

(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 J₁ at sufficiently large r₁

(74)

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

The total solution is obtained as the sum of four integrals

(75)

Here, the angles y₃ and y₄ corresponding to the apexes (l, -l) and (-l, l) are measured clockwise, whereas y₁ and y₂ are measured counter-clockwise. Substituting

r_j = [(x l)² + (y l)²]^1/2, y_j = arctan [(y l)/(x l)],

= [(x l)² + (y l)²]/2(x l)(y l),

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

where

(76)

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

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

[(x l)² + (y l)²]^1/2 = r + l( cosy + siny), r l,

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

q (x,y,t) = - (q₀ /4p²)(Dk/k_m) (2pi/k_m r)^1/2

(77)

where

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 k_ml 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.

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 l 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 k_ml 1, but this square is a composite radiator if k_ml > 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 r l (in the neighborhood of the center of the initial square), the following relationship takes place

[(x l)² + (y l)²]^1/2 = 2^1/2l + (2^1/2/2)r ( cosy siny),

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

(78)

where

As is seen from (78), a standing wave with square cells is formed in the central part (r 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 2^1/2p/k_m. 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 r l, we could use the expression (69) for F. We represents F in the form (70) to find the solution for r 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 k_ml 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, l 10. for k_m = 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, 0 z 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 q₀(x, y) are given on the straight line (x = 0), in the strip (-l 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 (k_x, k_y), 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 k_m

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 k_m is not equal to 2.7. However, the displacement u_z(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/d 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.